Skip to content

GitLab

Commit 46fe7aa1 authored by Matthieu Constant's avatar Matthieu Constant
Browse files

lots of bugs + testcase

parent d788dfc7
Show whitespace changes
Inline Side-by-side
.gitignore 0 → 100644
View file @ 46fe7aa1
output
*.so
*.core
*.msh
*.pos
__pycache__
*.pyc
*.aux
*.log
*.pdf
*.vtp
part-*
*.vtu
*.vtm
OUTBOX
DISPLAY
*.pvsm
DATBOX
*.png
.vscode/c_cpp_properties.json 0 → 100644
View file @ 46fe7aa1
{
"env" : {
"defaultIncludePath": [
"${workspaceFolder}",
"${workspaceFolder}/hxt",
"${workspaceFolder}/src",
"${workspaceFolder}/scontact"
],
"myCompilerPath": "/usr/bin/clang"
},
"configurations": [
{
"name": "Linux",
"intelliSenseMode": "clang-x64",
"includePath": [ "${defaultIncludePath}"],
"defines": [ "HXT_HAVE_PETSC","DIMENSION=2" ],
"compilerPath": "/usr/bin/clang",
"cStandard": "c11",
"cppStandard": "c++17",
"browse": {
"path": [ "${workspaceFolder}" ],
"limitSymbolsToIncludedHeaders": true,
"databaseFilename": ""
}
}
],
"version": 4
}
python/marblesbag/fluid.py
View file @ 46fe7aa1
......@@ -36,7 +36,7 @@ class StrongBoundary(Structure):
class fluid_problem :
def __init__(self, mesh_file_name, g, mu, rho, epsilon, strong_boundaries):
def __init__(self, mesh_file_name, g, mu, rho, epsilon, n_fluids, strong_boundaries):
nsb = len(strong_boundaries)
class Bnd :
def __init__(self, b) :
......@@ -48,11 +48,16 @@ class fluid_problem :
v[:] = self._b(x)
else :
v[:] = self._b
def np2c(a) :
tmp = numpy.require(a,"float","C")
r = c_void_p(tmp.ctypes.data)
r.tmp = tmp
return r
asb = (StrongBoundary*nsb)(*[StrongBoundary(i[0].encode(),i[1],BNDCB(Bnd(i[2]).apply)) for i in strong_boundaries])
self.asb = asb
lib.fluid_problem_new.restype = c_void_p
self._fp = c_void_p(lib.fluid_problem_new(mesh_file_name.encode(), c_double(g), c_double(mu), c_double(rho), c_double(epsilon), nsb, asb))
self._fp = c_void_p(lib.fluid_problem_new(mesh_file_name.encode(), c_double(g), np2c(mu), np2c(rho), c_double(epsilon), c_int(n_fluids), nsb, asb))
if self._fp == None :
raise NameError("cannot create fluid problem")
......
src/fluid_problem.c
View file @ 46fe7aa1
......@@ -195,21 +195,21 @@ static void compute_node_force_implicit(FluidProblem *problem, double dt, double
}
double du[DIMENSION], due[DIMENSION];
for (int j = 0; j < DIMENSION; ++j) {
du[j] = problem->particle_velocity[i*DIMENSION+j]-vf[j]/c;
due[j] = problem->particle_velocity[i*DIMENSION+j]-vfe[j]/c;
du[j] = problem->particle_velocity[i*DIMENSION+j]-vf[j]/(1-c);
due[j] = problem->particle_velocity[i*DIMENSION+j]-vfe[j]/(1-c);
}
double vol = problem->particle_volume[i];
double rhop = mass/vol;
double drho = rhop -problem->rho[ff];
double g = problem->g*drho/rhop;
double gamma = particle_drag_coeff(due,problem->mu[ff],problem->rho[ff],vol,c);
double gamma = particle_drag_coeff(due,problem->mu[ff],problem->rho[ff],vol,1-c);
double fcoeff = beta*mass/(mass+dt*gamma);
for (int d = 0; d < DIMENSION; ++d)
problem->particle_force[i*DIMENSION+d] = fcoeff*(-gamma*du[d]-vol*gradp[d]);
problem->particle_force[i*DIMENSION+1] += fcoeff*g*mass;
////pas bon
double *local_vector = &all_local_vector[N_SF*n_fields*problem->nFluids*iel+ff*N_SF*(DIMENSION+1)];
double *local_matrix = &all_local_matrix[N_SF*N_SF*n_fields*n_fields*problem->nFluids*problem->nFluids*iel+ff*N_SF*N_E*(DIMENSION+1)*(DIMENSION+1)];
for (int iphi = 0; iphi < N_SF; ++iphi) {
......@@ -225,7 +225,7 @@ static void compute_node_force_implicit(FluidProblem *problem, double dt, double
#define LOCAL_MATRIX(a,b) local_matrix[IDX+N2*a+b]
double f = fcoeff*phi[iphi];
for (int d = 0; d < DIMENSION; ++d){
LOCAL_MATRIX(d,d) -= -beta*f/c*phi[jphi]*gamma;
LOCAL_MATRIX(d,d) -= -beta*f/(1-c)*phi[jphi]*gamma;
LOCAL_MATRIX(d,DIMENSION) -= -beta*f*gamma*dt/mass*vol*dphi[jphi][d];
}
}
......@@ -259,7 +259,7 @@ static void fluid_problem_assemble_system(FluidProblem *problem, double *rhs, co
for (int iel=0; iel < mesh->n_elements; ++iel) {
const unsigned int *el = &mesh->elements[iel*N_N];
double dxdxi[DIMENSION][DIMENSION], dxidx[DIMENSION][DIMENSION], dphi[N_N][DIMENSION];
int P = n_fields;
int P = n_fields-1;
double *local_vector = &all_local_vector[N_SF*n_fields*n_fluids*iel];
double c=0,p=0, dp[DIMENSION]={0};
for (int i = 0; i < DIMENSION; ++i)
......@@ -270,6 +270,7 @@ static void fluid_problem_assemble_system(FluidProblem *problem, double *rhs, co
for (int iqp = 0; iqp< N_QP; ++iqp) {
const double jwl = QW[iqp]*detj;
double phi[N_SF];
shape_functions(QP[iqp],phi);
for (int i = 0; i < N_SF; ++i) {
double phii = phi[i];
p += solution[el[i]*n_fields+P]*phi[i];
......@@ -278,13 +279,12 @@ static void fluid_problem_assemble_system(FluidProblem *problem, double *rhs, co
for (int j = 0; j < DIMENSION; ++j)
dp[j] += dphi[i][j]*solution[el[i]*n_fields+P];
}
}
}printf("c=%g,p=%g\n",c,p);
double *local_matrix = &all_local_matrix[N_SF*N_SF*n_fields*n_fields*n_fluids*n_fluids*iel];
for (int ff = 0; ff < n_fluids; ++ff)
{
int Q = DIMENSION+ff*(DIMENSION+1);
double dq[DIMENSION]={0}, du[DIMENSION][DIMENSION]={{0}}, dqdq[DIMENSION]={0};
for (int i = 0; i< N_SF; ++i) {
for (int j = 0; j < DIMENSION; ++j) {
......@@ -297,16 +297,17 @@ static void fluid_problem_assemble_system(FluidProblem *problem, double *rhs, co
for (int iqp = 0; iqp< N_QP; ++iqp) {
double phi[N_SF];
shape_functions(QP[iqp],phi);
double u[DIMENSION]={0}, dudt[DIMENSION]={0}, dqdt=0,q=0;
double u[DIMENSION]={0}, dudt[DIMENSION]={0}, dqdt=0, q=0;
for (int i = 0; i < N_SF; ++i) {
for (int j = 0; j < DIMENSION; ++j){
u[j] += solution[el[i]*n_fields+j+ff*(DIMENSION+1)]*phi[i];
dudt[j] += (solution[el[i]*n_fields+j+ff*(DIMENSION+1)]-solution_old[el[i]*n_fields+j+ff*(DIMENSION+1)])/dt*phi[i];
}
q += solution[el[i]*n_fields+Q]*phi[i];
dqdt += (q-solution_old[el[i]*n_fields+Q]*phi[i])/dt;
dqdt += (solution[el[i]*n_fields+Q]-solution_old[el[i]*n_fields+Q])*phi[i]/dt;
//dcdt += (problem->porosity[el[i]]-problem->old_porosity[el[i]])/dt*phi[i];
}
printf("q=%g,dqdt=%g\n",q,dqdt);
const double jw = QW[iqp]*detj;
for (int iphi = 0; iphi < N_SF; ++iphi) {
double phii = phi[iphi];
......@@ -336,11 +337,10 @@ static void fluid_problem_assemble_system(FluidProblem *problem, double *rhs, co
);
}
local_vector[iphi+N_SF*P] += jw*phii*(q);
local_vector[iphi+N_SF*Q] += jw*phii*(-dqdt+divu)+jw*epsilon*dphidp;
local_vector[iphi+N_SF*Q] += jw*phii*(dqdt+divu)+jw*epsilon*dphidp;
for (int jphi = 0; jphi < N_SF; ++jphi) {
double phij = phi[jphi];
const double *dphij = dphi[jphi];
int P = DIMENSION;
//int N2 = N_SF*N_SF;
//int IDX = iphi*N2+jphi;
int N2 = n_fields*N_SF;
......@@ -357,14 +357,14 @@ static void fluid_problem_assemble_system(FluidProblem *problem, double *rhs, co
dphiidtau += dphii[j]*dtau[j];
}
for (int j = 0; j < DIMENSION; ++j) {
int U = j;
int U = j+ff*(DIMENSION+1);
LOCAL_MATRIX(U,P) += -jw*dphii[j]*phij;
LOCAL_MATRIX(Q,U) += jw*phii*dphij[j];
double utau =0;
double dphisigmadq = 0;
double dutaudq = 0;
for (int k = 0; k < DIMENSION; ++k) {
int V = k;
int V = k+ff*(DIMENSION+1);
utau += u[k]*tau[j][k];
dutaudq += u[k]*u[j]*(dq[k]*phij/(q*q)-dqdq[k]/q);
dphisigmadq += 2*mu[ff]*dphii[k]*0.5*(u[k]*dq[j]*phij/(q*q)+u[j]*dq[k]*phij/(q*q)-u[k]*dqdq[j]/q-u[j]*dqdq[k]/q);
......@@ -445,8 +445,10 @@ int fluid_problem_implicit_euler(FluidProblem *problem, double dt)
int n_fields = problem->nFields;
double *solution_old = malloc(mesh->n_nodes*n_fields*sizeof(double));
double *rhs = malloc(mesh->n_nodes*n_fields*sizeof(double));
for (int i=0; i<mesh->n_nodes*n_fields; ++i)
for (int i=0; i<mesh->n_nodes*n_fields; ++i){
solution_old[i] = problem->solution[i];
printf("old_solution=%g,solution=%g\n",solution_old[i],problem->solution[i]);}
problem->solution_explicit = solution_old;
double *dx = malloc(sizeof(double)*mesh->n_nodes*n_fields);
double norm0 = 0;
......@@ -457,6 +459,8 @@ int fluid_problem_implicit_euler(FluidProblem *problem, double dt)
for (int i=0; i<mesh->n_nodes*n_fields; ++i)
rhs[i] = 0;
fluid_problem_assemble_system(problem, rhs, solution_old, dt);
//exit(0);
if (i==1) break;
hxtLinearSystemGetRhsNorm(problem->linear_system,rhs,&norm);
printf("iter %i norm = %g\n", i, norm);
if (norm < newton_atol)
......@@ -472,6 +476,10 @@ int fluid_problem_implicit_euler(FluidProblem *problem, double dt)
for (int j=0; j<mesh->n_nodes*n_fields; ++j) {
problem->solution[j] -= dx[j];
}
for (int j=0; j<mesh->n_nodes*n_fields; ++j){
printf("sol=%g\n",dx[j]);
}
exit(0);
}
printf("total solve time : %g\n", totaltime);
free(dx);
......@@ -518,6 +526,9 @@ static void fluid_problem_compute_porosity(FluidProblem *problem)
for (int j=0; j<N_N;++j)
problem->compacity[el[j]] += sf[j]*volume[i];
}
for (int i = 0; i < mesh->n_nodes; ++i){
problem->compacity[i] = 0;
}
}
FluidProblem *fluid_problem_new(const char *mesh_file_name, double g, double *mu, double *rho, double epsilon, int n_fluids, int n_strong_boundaries, StrongBoundary *strong_boundaries)
......@@ -539,6 +550,7 @@ FluidProblem *fluid_problem_new(const char *mesh_file_name, double g, double *mu
problem->mu = malloc(n_fluids*sizeof(double));
for (int i = 0; i < n_fluids; ++i)
{
printf("mu=%g,rho=%g\n",mu[i],rho[i]);
problem->rho[i] = rho[i];
problem->mu[i] = mu[i];
}
......@@ -917,9 +929,13 @@ void fluid_problem_set_particles(FluidProblem *problem, int n, double *mass, dou
problem->old_compacity[i] = problem->compacity[i];
}
if (init) {
int Q = DIMENSION+1;
for (int i = 0; i < problem->mesh->n_nodes; ++i)
//initial fluid_0+grains
int Q = DIMENSION+0*(DIMENSION+1);
for (int i = 0; i < problem->mesh->n_nodes; ++i){
problem->solution[i*n_fields+Q] = 1-problem->compacity[i];
// printf("compacity=%g\n",problem->compacity[i]);
}
}
fluid_problem_compute_porosity(problem);
}
......
src/fluid_problem_io.c
View file @ 46fe7aa1
......@@ -38,7 +38,7 @@ int fluid_problem_export(const FluidProblem *problem, const char *output_dir, do
}
if(mesh_write_msh(problem->mesh, f))
return -1;
if(mesh_write_msh_scalar(problem->mesh, f, "porosity", t, iter, problem->porosity, 1, 0))
if(mesh_write_msh_scalar(problem->mesh, f, "compacity", t, iter, problem->compacity, 1, 0))
return -1;
int comp[]={0,1,2};
if(mesh_write_msh_vector(problem->mesh, f, "uv", t, iter, problem->solution, n_fields, comp))
......@@ -104,8 +104,8 @@ int fluid_problem_import_vtk(FluidProblem *problem, const char *filename)
problem->solution = malloc(sizeof(double)*(DIMENSION+1)*m->n_nodes);
while(mb_read_vtk_data_array(r,&a) != 0){
if(!strcasecmp(a.name,"Porosity")) {
free(problem->porosity);
problem->porosity = mb_vtk_data_array_to_double(&a,m->n_nodes,1);
free(problem->compacity);
problem->compacity = mb_vtk_data_array_to_double(&a,m->n_nodes,1);
}
else if(!strcasecmp(a.name,"Pressure")) {
double *p = mb_vtk_data_array_to_double(&a,m->n_nodes,1);
......@@ -220,21 +220,21 @@ int fluid_problem_export_vtk(const FluidProblem *problem, const char *output_dir
#endif
fprintf(f,"\n</DataArray>\n");
fprintf(f,"</Cells>\n");
fprintf(f,"<PointData Scalars=\"Porosity Pressure\" Vectors=\"Velocity\">\n");
fprintf(f,"<DataArray Name=\"Porosity\" NumberOfComponents = \"1\" type=\"Float64\" format=\"ascii\">\n");
fprintf(f,"<PointData Scalars=\"Compacity Pressure\" Vectors=\"Velocity\">\n");
fprintf(f,"<DataArray Name=\"Compacity\" NumberOfComponents = \"1\" type=\"Float64\" format=\"ascii\">\n");
for (int i = 0; i < mesh->n_nodes; ++i) {
fprintf(f, "%g ", problem->porosity[i]);
fprintf(f, "%g ", problem->compacity[i]);
}
fprintf(f,"\n</DataArray>\n");
fprintf(f,"<DataArray Name=\"Pressure\" NumberOfComponents = \"1\" type=\"Float64\" format=\"ascii\">\n");
for (int i = 0; i < mesh->n_nodes; ++i) {
fprintf(f, "%g ", problem->solution[i*(DIMENSION+2)+DIMENSION]);
fprintf(f, "%g ", problem->solution[i*(problem->nFields)+problem->nFields-1]);
}
fprintf(f,"\n</DataArray>\n");
fprintf(f,"<DataArray Name=\"Velocity\" NumberOfComponents = \"3\" type=\"Float64\" format=\"ascii\">\n");
for (int i = 0; i < mesh->n_nodes; ++i) {
for (int j = 0; j < DIMENSION; ++j)
fprintf(f, "%g ", problem->solution[i*(DIMENSION+2)+j]);
fprintf(f, "%g ", problem->solution[i*(problem->nFields)+j]);
if(DIMENSION==2)
fprintf(f, "0 ");
}
......
testcases/drop-2d/SimpleDrop/drop.py
View file @ 46fe7aa1
......@@ -53,6 +53,8 @@ def genInitialPosition(filename, r, rout, rhop, compacity) :
for i in range(x.shape[0]) :
#Add grains object with properties defined above
p.add_particle((x[i]+2*random.uniform(-e/2+r,e/2-r), y[i]+2*random.uniform(-e/2+r,e/2-r)), r, r**2 * np.pi * rhop);
if i==0:
break;
#Vertical shift of the grains to the top of the box
p.position()[:, 1] += 0.52
p.write_vtk(filename,0,0)
......@@ -98,8 +100,8 @@ outf = 1 # number of iterations between o
#Object fluid creation + Boundary condition of the fluid (field 0 is horizontal velocity; field 1 is vertical velocity; field 2 is pressure)
#Format: strong_boundaries = [(Boundary tag, Fluid field, Value)
strong_boundaries = [("Top",2,0.),("Top",1,0.),("Bottom",1,0.),("Lateral",0,0.)]
fluid = fluid.fluid_problem("mesh.msh",g,nu*rho,rho,epsilon,strong_boundaries)
strong_boundaries = [("Top",3,0.),("Top",1,0.),("Bottom",1,0.),("Lateral",0,0.),("Top",2,1.),("Bottom",2,1.),("Lateral",2,1.)]
fluid = fluid.fluid_problem("mesh.msh",g,[nu*rho],[rho],epsilon,1,strong_boundaries)
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity(),True)
ii = 0
......
testcases/drop-2d/SimpleDrop/mesh.geo
View file @ 46fe7aa1
......@@ -28,8 +28,8 @@ Field[1].NodesList = {5};
Field[2] = Threshold;
Field[2].DistMax = 0.08;
Field[2].DistMin = 0.02;
Field[2].LcMax = 0.02;
Field[2].LcMin = 0.001;
Field[2].LcMax = 0.08;
Field[2].LcMin = 0.004;
Field[2].IField = 1;
Background Field = 2;
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment