Commit e0b488f1 authored by Jonathan Lambrechts's avatar Jonathan Lambrechts
Browse files

broken

parent edafee0e
Pipeline #4689 passed with stage
in 23 seconds
......@@ -22,6 +22,7 @@
set(SRC
fluid_problem.c
fluid_problem_concentration.c
mesh.c
mesh_find.c
hxt_linear_system.c
......
This diff is collapsed.
......@@ -43,7 +43,7 @@ typedef struct {
} StrongBoundary;
typedef struct FluidProblem FluidProblem;
typedef void f_cb(FluidProblem *problem,const double *n, double *f,const double *s,const double *ds,const double c,const double *dc,const double dt, int eid, const double *v);
typedef void f_cb(FluidProblem *problem,const double *n, double *f,const double *s,const double *ds,const double c,const double *dc,const double rho, const double mu, const double dt, int eid, const double *v);
typedef struct {
char *tag;
f_cb *cb;
......@@ -61,6 +61,8 @@ struct FluidProblem {
MeshTree *mesh_tree;
double *porosity;
double *oldporosity;
double *concentration;
double *oldconcentration;
double *solution;
double *solution_explicit;
double *node_volume;
......@@ -89,6 +91,7 @@ struct FluidProblem {
void fluid_problem_compute_node_particle_force(FluidProblem *problem, double dt, double *particle_force);
int fluid_problem_implicit_euler(FluidProblem *problem, double dt, double newton_atol, double newton_rtol, int newton_max_it, int reduced_gravity, double stab_param);
void fluid_problem_set_particles(FluidProblem *problem, int n, double *mass, double *volume, double *position, double *velocity, long int *elid, int reload);
void fluid_problem_set_concentration_cg(FluidProblem *problem, double *concentration);
void fluid_problem_free(FluidProblem *problem);
FluidProblem *fluid_problem_new(double g, int n_fluids, double *mu, double *rho, double coeffStab, const char *petsc_solver_type);
void fluid_problem_load_msh(FluidProblem *problem, const char *mesh_file_name);
......@@ -106,10 +109,12 @@ int fluid_problem_n_elements(const FluidProblem *p);
int *fluid_problem_elements(const FluidProblem *p);
double *fluid_problem_old_porosity(const FluidProblem *p);
double *fluid_problem_porosity(const FluidProblem *p);
double *fluid_problem_concentration_dg(const FluidProblem *p);
void fluid_problem_set_wall_boundary(FluidProblem *p, const char *tag, BoundaryCallback *pressurecb);
int fluid_problem_private_n_physical(const FluidProblem *p);
void fluid_problem_private_physical(const FluidProblem *p, int i, char **phys_name, int *phys_dim, int *phys_id, int *phys_n_nodes, int **phys_nodes, int *phys_n_boundaries, int **phys_boundaries);
void fluid_problem_private_set_elements(FluidProblem *p, int n_nodes, double *x, int n_elements, int *elements);
void fluid_problem_private_add_physical(FluidProblem *p, int tag, int dim, const char *name, int n_nodes, int *nodes, int n_bnd, int *bnd);
void fluid_problem_advance_concentration(FluidProblem *problem, double dt);
#endif
#include "fluid_problem.h"
#include "fluid_problem_fem.h"
#include <math.h>
void fluid_problem_advance_concentration(FluidProblem *problem, double dt) {
const Mesh *mesh = problem->mesh;
const int n_fields = fluid_problem_n_fields(problem);
double *velocity = malloc(D*problem->mesh->n_nodes*sizeof(double));
double *rhs = malloc(mesh->n_elements*N_N*sizeof(double));
double *det = malloc(mesh->n_elements*sizeof(double));
for (int i = 0; i < mesh->n_elements*N_N; ++i) {
rhs[i] = 0;
}
for (int in = 0; in < mesh->n_nodes; ++in) {
for (int id = 0; id < D; ++id) {
velocity[in*D+id] = problem->solution[in*n_fields+id]/problem->porosity[in];
}
}
/*volume term*/
for (int iel=0; iel < mesh->n_elements; ++iel) {
const unsigned int *el = &mesh->elements[iel*N_N];
double dxdxi[D][D], dxidx[D][D], dphi[N_N][D];
for (int i = 0; i < D; ++i)
for (int j = 0; j < D; ++j)
dxdxi[i][j] = mesh->x[el[j+1]*3+i] - mesh->x[el[0]*3+i];
const double detj = invDD(dxdxi, dxidx);
det[iel] = detj;
grad_shape_functions(dxidx, dphi);
double da[D] = {0};
for (int i = 0; i < N_SF; ++i) {
for (int id = 0; id < D; ++id) {
da[id] += dphi[i][id]*problem->concentration[iel*N_N+i];
}
}
for (int iqp = 0; iqp< N_QP; ++iqp) {
double phi[N_SF];
shape_functions(QP[iqp],phi);
const double jw = QW[iqp]*detj;
double u[D] = {0};
for (int i = 0; i < N_SF; ++i) {
for (int id = 0; id < D; ++id) {
u[id] += phi[i]*velocity[el[i]*D+id];
}
}
for (int i = 0; i < N_SF; ++i) {
double r = 0;
for (int id = 0; id < D; ++id) {
r -= u[id]*da[id];
}
rhs[iel*N_N+i] += r*phi[i]*jw;
}
}
}
/*interface_term*/
#if DIMENSION == 2
int elbnd[][2] = {{0,1},{1,2},{2,0},{1,0},{2,1},{0,2}};
#else
int elbnd[4][3] = {{0,1,2},{0,2,3},{0,3,1},{1,3,2}};
#endif
for (int ii = 0; ii < mesh->n_interfaces; ++ii) {
const int *interface = &mesh->interfaces[ii*4];
const int iel0 = interface[0];
const int iel1 = interface[2];
if (iel1 < 0) continue;
const int icl0 = interface[1];
const int icl1 = interface[3];
int nodes[D];
double *x[D];
const unsigned int *el0 = &mesh->elements[iel0*N_N];
const unsigned int *el1 = &mesh->elements[iel1*N_N];
for (int i = 0; i < D; ++i){
nodes[i] = el0[elbnd[icl0][i]];
x[i] = &mesh->x[nodes[i]*3];
}
double dxdxi[D][D], dxidx[D][D], dphi[N_N][D];
for (int i = 0; i < D; ++i)
for (int j = 0; j < D; ++j)
dxdxi[i][j] = mesh->x[el0[j+1]*3+i] - mesh->x[el0[0]*3+i];
invDD(dxdxi, dxidx);
grad_shape_functions(dxidx, dphi);
#if DIMENSION == 2
const double dx[2] = {x[1][0]-x[0][0], x[1][1]-x[0][1]};
const double detbnd = sqrt(dx[0]*dx[0] + dx[1]*dx[1]);
const double n[2] = {dx[1]/detbnd,-dx[0]/detbnd};
#else
const double dx[2][3] = {{x[1][0]-x[0][0], x[1][1]-x[0][1], x[1][2]-x[0][2]},{x[2][0]-x[0][0], x[2][1]-x[0][1], x[2][2]-x[0][2]}};
const double nn[3] = {
dx[0][1]*dx[1][2] - dx[0][2]*dx[1][1],
dx[0][2]*dx[1][0] - dx[0][0]*dx[1][2],
dx[0][0]*dx[1][1] - dx[0][1]*dx[1][0]};
const double detbnd = sqrt(nn[0]*nn[0] + nn[1]*nn[1] + nn[2]*nn[2]);
const double n[3] = {-nn[0]/detbnd,-nn[1]/detbnd,-nn[2]/detbnd};
#endif
for (int iqp = 0; iqp < N_LQP; ++iqp) {
double xi0[D],xi1[D];
#if DIMENSION == 2
double z = LQP[iqp];
switch(icl0) {
case 0 : xi0[0] = z; xi0[1] = 0; break;
case 1 : xi0[0] = 1-z; xi0[1] = z; break;
case 2 : xi0[0] = 0; xi0[1] = 1-z; break;
case 3 : xi0[0] = 1-z; xi0[1] = 0; break;
case 4 : xi0[0] = z; xi0[1] = 1-z; break;
case 5 : xi0[0] = 0; xi0[1] = z; break;
default:;
}
switch(icl1) {
case 0 : xi1[0] = z; xi1[1] = 0; break;
case 1 : xi1[0] = 1-z; xi1[1] = z; break;
case 2 : xi1[0] = 0; xi1[1] = 1-z; break;
case 3 : xi1[0] = 1-z; xi1[1] = 0; break;
case 4 : xi1[0] = z; xi1[1] = 1-z; break;
case 5 : xi1[0] = 0; xi1[1] = z; break;
default:;
}
#else
double zx = LQP[iqp][0];
double zy = LQP[iqp][1];
switch(icl0) {
case 0 : xi0[0] = zx; xi0[1] = zy; xi0[2] = 0; break;
case 1 : xi0[0] = 0; xi0[1] = zx; xi0[2] = zy; break;
case 2 : xi0[0] = zy; xi0[1] = 0; xi0[2] = zx; break;
case 3 : xi0[0] = 1-zx-zy; xi0[1] = zy; xi0[2] = zx; break;
default:;
}
switch(icl1) {
case 0 : xi1[0] = zx; xi1[1] = zy; xi1[2] = 0; break;
case 1 : xi1[0] = 0; xi1[1] = zx; xi1[2] = zy; break;
case 2 : xi1[0] = zy; xi1[1] = 0; xi1[2] = zx; break;
case 3 : xi1[0] = 1-zx-zy; xi1[1] = zy; xi1[2] = zx; break;
default:;
}
#endif
double phi0[N_SF], phi1[N_SF];
shape_functions(xi0,phi0);
shape_functions(xi1,phi1);
double a0 = 0, a1 = 0, un = 0;
for (int i = 0; i < N_SF; ++i) {
a0 += problem->concentration[iel0*N_N+i]*phi0[i];
a1 += problem->concentration[iel1*N_N+i]*phi1[i];
for (int id = 0; id < D; ++id) {
un -= velocity[el0[i]*D+id]*phi0[i]*n[id];
}
}
const double a = un > 0 ? a0 : a1;
const double r0 = -un*(a-a0);
const double r1 = un*(a-a1);
const double jw = LQW[iqp]*detbnd;
for (int iphi = 0; iphi < N_SF; ++iphi){
rhs[iel0*N_SF+iphi] += phi0[iphi]*r0*jw;
rhs[iel1*N_SF+iphi] += phi1[iphi]*r1*jw;
}
}
}
/*boundary_term*/
/*resolution*/
#if DIMENSION == 2
const double invmassmatrix[3][3] = {
{18,-6,-6},
{-6,18,-6},
{-6,-6,18}};
#else
const double invmassmatrix[4][4] = {
{ 96,-24,-24,-24},
{-24, 96,-24,-24},
{-24,-24, 96,-24},
{-24,-24,-24, 96}};
#endif
for (int iel = 0; iel < mesh->n_elements; ++iel) {
for (int in = 0; in < N_SF; ++in) {
for (int jn = 0; jn < N_SF; ++jn) {
problem->concentration[iel*N_SF+in] += invmassmatrix[in][jn]*rhs[iel*N_SF+jn]/det[iel]*dt;
}
}
}
free(velocity);
free(rhs);
}
#define D DIMENSION
#define deps 1e-8
#if D==2
#define N_SF 3
#define N_N 3
#define N_QP 3
#define N_LQP 2
static double LQP[] = {0.21132486540518713, 0.7886751345948129};
static double LQW[] = {0.5,0.5};
static double QP[][2] = {{1./6,1./6},{1./6,2./3},{2./3,1./6}};
static double QW[] = {1./6,1./6,1./6};
static void shape_functions(double *uvw, double *sf)
{
sf[0] = 1-uvw[0]-uvw[1];
sf[1] = uvw[0];
sf[2] = uvw[1];
}
static void grad_shape_functions(const double dxidx[2][2], double gsf[3][2])
{
for (int i = 0; i < 2; ++i) {
gsf[0][i] = -dxidx[0][i]-dxidx[1][i];
gsf[1][i] = dxidx[0][i];
gsf[2][i] = dxidx[1][i];
}
}
static double detDD(const double m[2][2])
{
return m[0][0]*m[1][1] - m[0][1]*m[1][0];
}
static double invDD(const double m[2][2], double inv[2][2])
{
double d = detDD(m);
inv[0][0] = m[1][1]/d;
inv[0][1] = -m[0][1]/d;
inv[1][0] = -m[1][0]/d;
inv[1][1] = m[0][0]/d;
return d;
}
#else
#define N_SF 4
#define N_N 4
#define N_QP 5
#define N_LQP 3
static double QP[][3] = {
{0.25, 0.25, 0.25},
{0.5, 1./6, 1./6},
{1./6, 0.5, 1./6},
{1./6, 1./6, 0.5},
{1./6, 1./6, 1./6}
};
static double LQP[][2] = {{1./6,1./6},{1./6,2./3},{2./3,1./6}};
static double QW[] = {-16./120, 9./120, 9./120, 9./120, 9./120};
static double LQW[] = {1./6,1./6,1./6};
static void shape_functions(double *uvw, double *sf)
{
sf[0] = 1-uvw[0]-uvw[1]-uvw[2];
sf[1] = uvw[0];
sf[2] = uvw[1];
sf[3] = uvw[2];
}
static void grad_shape_functions(const double dxidx[3][3], double gsf[4][3])
{
for (int i = 0; i < 3; ++i) {
gsf[0][i] = -dxidx[0][i]-dxidx[1][i]-dxidx[2][i];
gsf[1][i] = dxidx[0][i];
gsf[2][i] = dxidx[1][i];
gsf[3][i] = dxidx[2][i];
}
};
static double detDD(const double m[3][3])
{
return m[0][0]*(m[1][1]*m[2][2]-m[2][1]*m[1][2])
-m[1][0]*(m[0][1]*m[2][2]-m[2][1]*m[0][2])
+m[2][0]*(m[0][1]*m[1][2]-m[1][1]*m[0][2]);
}
static double invDD(const double m[3][3], double inv[3][3]) {
double d = detDD(m);
inv[0][0] = (m[1][1]*m[2][2]-m[2][1]*m[1][2])/d;
inv[0][1] = -(m[0][1]*m[2][2]-m[2][1]*m[0][2])/d;
inv[0][2] = (m[0][1]*m[1][2]-m[1][1]*m[0][2])/d;
inv[1][0] = -(m[1][0]*m[2][2]-m[2][0]*m[1][2])/d;
inv[1][1] = (m[0][0]*m[2][2]-m[2][0]*m[0][2])/d;
inv[1][2] = -(m[0][0]*m[1][2]-m[1][0]*m[0][2])/d;
inv[2][0] = (m[1][0]*m[2][1]-m[2][0]*m[1][1])/d;
inv[2][1] = -(m[0][0]*m[2][1]-m[2][0]*m[0][1])/d;
inv[2][2] = (m[0][0]*m[1][1]-m[1][0]*m[0][1])/d;
return d;
};
#endif
......@@ -27,6 +27,62 @@
#include <string.h>
#include "tools.h"
int edgecmp(const void *pe0, const void *pe1) {
int *e0 = (int*)pe0;
int *e1 = (int*)pe1;
if (e0[0] != e1[0]) return e0[0]-e1[0];
return e0[1]-e1[1];
}
static void mesh_gen_edges(Mesh *m) {
#if DIMENSION == 2
int n_half_edges = m->n_elements*3;
int *halfedges = malloc(sizeof(int)*n_half_edges*4);
int edgemap[3][2] = {{0,1},{1,2},{2,0}};
for (int i=0; i<m->n_elements; ++i) {
for (int j = 0; j < 3; ++j) {
int n0 = m->elements[i*3+edgemap[j][0]], n1 = m->elements[i*3+edgemap[j][1]];
if (n0 < n1) {
halfedges[(i*3+j)*4+0] = n0;
halfedges[(i*3+j)*4+1] = n1;
halfedges[(i*3+j)*4+3] = j;
}
else {
halfedges[(i*3+j)*4+0] = n1;
halfedges[(i*3+j)*4+1] = n0;
halfedges[(i*3+j)*4+3] = j+3;
}
halfedges[(i*3+j)*4+2] = i;
}
}
qsort(halfedges,n_half_edges,4*sizeof(int),edgecmp);
int n_edges = 0;
for (int i = 0; i < n_half_edges; ++i) {
if (i == 0 || edgecmp(halfedges+4*(i-1),halfedges+4*i)) n_edges++;
}
int *edges = malloc(n_edges*4*sizeof(int));
n_edges = 0;
for (int i = 0; i < n_half_edges; ++i) {
if (i == 0 || edgecmp(halfedges+4*(i-1),halfedges+4*i)){
n_edges += 1;
edges[(n_edges-1)*4+0] = halfedges[4*i+2];
edges[(n_edges-1)*4+1] = halfedges[4*i+3];
edges[(n_edges-1)*4+2] = -1;
edges[(n_edges-1)*4+3] = -1;
}
else {
edges[(n_edges-1)*4+2] = halfedges[4*i+2];
edges[(n_edges-1)*4+3] = halfedges[4*i+3];
}
}
m->n_interfaces = n_edges;
m->interfaces = edges;
free(halfedges);
#else
// todo Matthieu
#endif
}
static int check_word_in_file(FILE *f, const char *w) {
char word[256];
if(fscanf(f,"%255s", word) != 1){
......@@ -131,6 +187,7 @@ int bnd_element_cmp(const void *ps0, const void *ps1) {
Mesh *mesh_load_msh(const char *filename)
{
Mesh *m = malloc(sizeof(Mesh));
m->interfaces = NULL;
FILE *f = fopen(filename, "r");
if (!f){
printf("Cannot open file \"%s\".\n", filename);
......@@ -332,6 +389,7 @@ Mesh *mesh_load_msh(const char *filename)
}
}
free(bnd_elements);
mesh_gen_edges(m);
return m;
}
......@@ -375,6 +433,7 @@ void mesh_free(Mesh *m)
free(m->phys_dimension);
free(m->phys_id);
free(m->phys_n_nodes);
free(m->interfaces);
free(m);
}
......
......@@ -40,6 +40,8 @@ typedef struct {
int **phys_boundaries; // eid, side
int *phys_dimension;
int *phys_id;
int n_interfaces;
int *interfaces; // eid0, closure0, eid1, closure1
} Mesh;
Mesh *mesh_load_msh(const char *filename);
......
......@@ -65,7 +65,7 @@ def _is_string(s) :
class FluidProblem :
"""Create numerical representation of the fluid."""
def __init__(self, dim, g, mu, rho, coeff_stab=0.01, petsc_solver_type="-pc-type lu"):
def __init__(self, dim, g, mu, rho, coeff_stab=0.01, petsc_solver_type="-pc_type lu"):
"""Build the fluid structure.
Keyword arguments:
......@@ -193,6 +193,7 @@ class FluidProblem :
it -- computation iteration
"""
v = self.solution()[:,:self._dim]
cell_data = []
if self._dim == 2 :
v = np.insert(v,self._dim,0,axis=1)
data = [
......@@ -202,7 +203,7 @@ class FluidProblem :
("old_porosity",self.old_porosity())
]
if self._n_fluids == 2 :
data.append(("concentration",self.solution()[:,[self._dim+1]]))
cell_data.append(("concentration",self.concentration_dg()))
field_data = []
for i in range(self._n_physicals()) :
name,dim,pid,nodes,bnd = self._physical(i)
......@@ -211,7 +212,7 @@ class FluidProblem :
connectivity = self.elements()
types = np.repeat([5 if self._dim == 2 else 10],connectivity.shape[0])
offsets = np.cumsum(np.repeat([self._dim+1],connectivity.shape[0]))
VTK.write(output_dir+"/fluid",it,t,(types,connectivity,offsets),self.coordinates(),data,field_data)
VTK.write(output_dir+"/fluid",it,t,(types,connectivity,offsets),self.coordinates(),data,field_data,cell_data)
def import_vtk(self, filename) :
"""Read output file to reload computation.
......@@ -271,7 +272,9 @@ class FluidProblem :
newton_rtol -- optional argument to specify relative tolerance for nonlinear solver
newton_max_it -- optional argument to specify maximum number of iterations for nonlinear solver
"""
return self._lib.fluid_problem_implicit_euler(self._fp, c_double(dt), c_double(newton_atol), c_double(newton_rtol), c_int(newton_max_it), c_int(reduced_gravity), c_double(stab_param))
#nit = self._lib.fluid_problem_implicit_euler(self._fp, c_double(dt), c_double(newton_atol), c_double(newton_rtol), c_int(newton_max_it), c_int(reduced_gravity), c_double(stab_param))
self._lib.fluid_problem_advance_concentration(self._fp,c_double(dt))
#return nit
def set_particles(self, mass, volume, position, velocity, reload = False):
"""Set location of the grains in the mesh and compute the porosity in each cell.
......@@ -327,6 +330,20 @@ class FluidProblem :
"""Return the porosity at nodes"""
return self._get_matrix("porosity", self.n_nodes(), 1)
def set_concentration_cg(self,concentration):
"""Set the concentration at nodes"""
def np2c(a) :
tmp = np.require(a,"float","C")
r = c_void_p(tmp.ctypes.data)
r.tmp = tmp
return r
concentration = concentration.reshape((self.n_nodes(),1))
self._lib.fluid_problem_set_concentration_cg(self._fp,np2c(concentration))
def concentration_dg(self):
"""Return the concentration at discontinuous nodes"""
return self._get_matrix("concentration_dg", self.n_elements(), self._dim+1)
def old_porosity(self):
"""Return the porosity at previous time step"""
return self._get_matrix("old_porosity", self.n_nodes(), 1)
......
L = 0.05;
H = 0.6;
y = 0;
lc = 0.001;
Point(1) = {-L, H, 0};
Point(2) = {-L, -H, 0};
Point(3) = {L, -H, 0};
Point(4) = {L, H, 0};
Line(1) = {1, 2};
Line(2) = {2, 3};
Line(3) = {3, 4};
Line(4) = {4, 1};
Line Loop(1) = {1:4};
Plane Surface(1) = {1};
Physical Line("Bottom") = {2};
Physical Line("Lateral") = {1,3};
Physical Line("Top") = {4};
Physical Surface("Domain") = {1};
Physical Point("PtFix") = {1};
Mesh.Algorithm = 5;
Merge "lc.pos";
Field[1] = PostView;
Field[1].IView = 0;
Background Field = 1;
Mesh.CharacteristicLengthExtendFromBoundary = 0;
Mesh.CharacteristicLengthFromPoints = 0;
Mesh.MshFileVersion = 2;
\ No newline at end of file
# MigFlow - Copyright (C) <2010-2018>
# <Universite catholique de Louvain (UCL), Belgium
# Universite de Montpellier, France>
#
# List of the contributors to the development of MigFlow: see AUTHORS file.
# Description and complete License: see LICENSE file.
#
# This program (MigFlow) is free software:
# you can redistribute it and/or modify it under the terms of the GNU Lesser General
# Public License as published by the Free Software Foundation, either version
# 3 of the License, or (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU Lesser General Public License for more details.
#
# You should have received a copy of the GNU Lesser General Public License
# along with this program (see COPYING and COPYING.LESSER files). If not,
# see <http://www.gnu.org/licenses/>.
#!/usr/bin/env python
from migflow import fluid
from migflow import scontact
import numpy as np
import os
import time
import shutil
import random
outputdir = "output"
if not os.path.isdir(outputdir) :
os.makedirs(outputdir)
t = 0
ii = 0
#physical parameters
g = -9.81
rhof = 1030 # fluid density
rhop = 2450 # grains density
r=154e-6 # grains radii
R = 3.3e-3 # drop radius
compacity = .20 # solid volume fraction in the drop
phi = 1 - compacity
nuf = 1.17/rhof # kinematic viscosity
muf = nuf*rhof # dynamic viscosity
##mixture properties
rhom = (1-phi)*rhop+phi*rhof #mixture density