Skip to content
GitLab
Commit cd61446c authored by Jonathan Lambrechts's avatar Jonathan Lambrechts
Browse files

force implicit 2d

parent 06d91fdd
Hide whitespace changes
Inline Side-by-side
marblesbag/src/fluid_problem.c
View file @ cd61446c
#include "tools.h" #include "tools.h"
#include <stdlib.h> #include <stdlib.h>
#include <stdio.h> #include <stdio.h>
...@@ -100,212 +99,37 @@ static double invDD(const double m[3][3], double inv[3][3]) { ...@@ -100,212 +99,37 @@ static double invDD(const double m[3][3], double inv[3][3]) {
#endif #endif
#if DIMENSION == 2 #if DIMENSION==2
static void particle_drag(double u[2], double mu, double rho, double d, double c, double drag[2], double dt, double mass, double rhop, double gradp[2]) static double particle_drag_coeff(double u[2], double mu, double rho, double vol, double c)
{ {
double d = 2*sqrt(vol/M_PI);
double normu = hypot(u[0],u[1]); double normu = hypot(u[0],u[1]);
//Reynolds/|u_p-u_f| //Reynolds/|u_p-u_f|
double Re_O_u = sqrt(d)*c/mu*rho; double Re_O_u = sqrt(d)*c/mu*rho;
double Cd_u = 0.63*sqrt(normu)+4.8/sqrt(Re_O_u); double Cd_u = 0.63*sqrt(normu)+4.8/sqrt(Re_O_u);
Cd_u = Cd_u*Cd_u; Cd_u = Cd_u*Cd_u;
double f = pow(c,-1.8); double f = pow(c,-1.8);
double r = d/2.; return = f*Cd_u*rho/2*d;
double GoU = f*Cd_u*rho/2*(2*r);
double a = 0;
double h = (1+(1-c)*rhop/(c*rho)+a)/(a+1);
double F = (mass*GoU)/(mass+h*dt*GoU);
drag[0] = (u[0]-gradp[0]*dt/rhop)*F;
drag[1] = (u[1]-9.81*dt*(1-rho/rhop)-gradp[1]*dt/rhop)*F;
} }
#else #else
static void particle_drag(double u[3], double mu, double rho, double d, double c, double drag[3], double dt, double mass, double rhop, double gradp[3]) static double particle_drag_coeff(double u[3], double mu, double rho, double vol, double c)
{ {
double normu = pow(u[0]*u[0]+u[1]*u[1]+u[2]*u[2],0.5); double d = 2*pow(3./4.*vol/M_PI,1./3.);
double normu = sqrt(u[0]*u[0]+u[1]*u[1]+u[2]*u[2]);
//Reynolds/|u_p-u_f| //Reynolds/|u_p-u_f|
double Re_O_u = d*c/mu*rho; double Re_O_u = d*c/mu*rho;
double Cd_u = 0.63*sqrt(normu)+4.8/sqrt(Re_O_u); double Cd_u = 0.63*sqrt(normu)+4.8/sqrt(Re_O_u);
Cd_u = Cd_u*Cd_u; Cd_u = Cd_u*Cd_u;
double f = pow(c,-1.8); double f = pow(c,-1.8);
double r = d/2.; double r = d/2.;
double GoU = f*Cd_u*rho/2*(M_PI*r*r); return f*Cd_u*rho/2*(M_PI*r*r);
double a = 0;
double h = (1+(1-c)*rhop/(c*rho)+a)/(a+1);
double F = (mass*GoU)/(mass+h*dt*GoU);
drag[0] = (u[0]-gradp[0]*dt/rhop)*F;
drag[1] = (u[1]-9.81*dt*(1-rho/rhop)-gradp[1]*dt/rhop)*F;
drag[2] = (u[2]-gradp[2]*dt/rhop)*F;
} }
#endif #endif
#if DIMENSION == 2
void fluid_problem_compute_node_particle_force(FluidProblem *problem, double dt, double *particle_force) {
double *porosity = problem->porosity;
Mesh *mesh = problem->mesh;
const double *solution = problem->solution;
for (int i = 0; i < problem->n_particles; ++i) {
double *xi = problem->particle_uvw + 2*i;
double phi[3] = {1-xi[0]-xi[1],xi[0],xi[1]};
int iel = problem->particle_element_id[i];
double vol = problem->particle_volume[i];
double mass = problem->particle_mass[i];
particle_force[i*2+0] = 0;
particle_force[i*2+1] = problem->g*(mass-problem->rho*vol);
if(iel < 0){
continue;
}
int *tri = problem->mesh->elements+iel*3;
double C[] = {porosity[tri[0]],porosity[tri[1]],porosity[tri[2]]};
double U[] = {solution[tri[0]*3+0],solution[tri[1]*3+0],solution[tri[2]*3+0]};
double V[] = {solution[tri[0]*3+1],solution[tri[1]*3+1],solution[tri[2]*3+1]};
double P[] = {solution[tri[0]*3+2],solution[tri[1]*3+2],solution[tri[2]*3+2]};
double vf[2] = {
phi[0]*U[0]+phi[1]*U[1]+phi[2]*U[2],
phi[0]*V[0]+phi[1]*V[1]+phi[2]*V[2]
};
double c = phi[0]*C[0]+phi[1]*C[1]+phi[2]*C[2];
const double x[3][2] = {
{mesh->x[tri[0]*3+0],mesh->x[tri[0]*3+1]},
{mesh->x[tri[1]*3+0],mesh->x[tri[1]*3+1]},
{mesh->x[tri[2]*3+0],mesh->x[tri[2]*3+1]}};
const double dxdxi[2][2] = {
{x[1][0]-x[0][0],x[2][0]-x[0][0]},
{x[1][1]-x[0][1],x[2][1]-x[0][1]}};
const double detj = dxdxi[0][0]*dxdxi[1][1]-dxdxi[0][1]*dxdxi[1][0];
const double dxidx[2][2] = {
{dxdxi[1][1]/detj,-dxdxi[0][1]/detj},
{-dxdxi[1][0]/detj,dxdxi[0][0]/detj}
};
const double dphi[][2] = {
{-dxidx[0][0]-dxidx[1][0],-dxidx[0][1]-dxidx[1][1]},
{dxidx[0][0],dxidx[0][1]},
{dxidx[1][0],dxidx[1][1]}
};
double gradp[2] = {
dphi[0][0]*P[0]+dphi[1][0]*P[1]+dphi[2][0]*P[2],
dphi[0][1]*P[0]+dphi[1][1]*P[1]+dphi[2][1]*P[2]
};
double du[2] = {
problem->particle_velocity[i*2+0]-vf[0]/c,
problem->particle_velocity[i*2+1]-vf[1]/c
};
double rhop = mass/vol;
double d = 2*sqrt(vol/M_PI);
double drag[2];
particle_drag(du,problem->mu,problem->rho,d,c,drag, dt, mass, rhop, gradp);
/*FILE *log = fopen("DragForce.txt", "at");
if (!log) log = fopen("DragForce.txt", "wt");
fprintf(log, "%d, %e, %e, %e, %e, %e, %e\n", i, du[0], du[1], drag[0], drag[1], mass, vol);
fclose(log);*/
particle_force[i*2+0] += -drag[0]-vol*gradp[0];
particle_force[i*2+1] += -drag[1]-vol*gradp[1];
//printf("%e,%e,%e,%e,%e,%e\n",particle_force[i*2+0],particle_force[i*2+1],-drag[0],-drag[1],0.,problem->g*(mass-problem->rho*vol));
for (int iphi = 0; iphi < 3; ++iphi) {
problem->node_force[tri[iphi]*2+0] += drag[0]*phi[iphi];
problem->node_force[tri[iphi]*2+1] += drag[1]*phi[iphi];
}
}
for (int i = 0; i < mesh->n_nodes; ++i) {
problem->node_force[i*2+0] /= problem->node_volume[i];
problem->node_force[i*2+1] /= problem->node_volume[i];
}
}
#else
void fluid_problem_compute_node_particle_force(FluidProblem *problem, double dt, double *particle_force) { void fluid_problem_compute_node_particle_force(FluidProblem *problem, double dt, double *particle_force) {
for (int i = 0; i < problem->n_particles*DIMENSION; ++i) { for (int i = 0; i < problem->n_particles*DIMENSION; ++i) {
particle_force[i] = problem->particle_force[i]; particle_force[i] = problem->particle_force[i];
} }
/*
double *porosity = problem->porosity;
Mesh *mesh = problem->mesh;
const double *solution = problem->solution;
for (int i = 0; i < mesh->n_nodes*3; ++i) {
problem->node_force[i] = 0.0;
}
for (int i = 0; i < problem->n_particles; ++i) {
double *xi = problem->particle_uvw + 3*i;
double phi[4] = {1-xi[0]-xi[1]-xi[2],xi[0],xi[1],xi[2]};
int iel = problem->particle_element_id[i];
double vol = problem->particle_volume[i];
double mass = problem->particle_mass[i];
particle_force[i*3+0] = 0;
particle_force[i*3+1] = problem->g*(mass-problem->rho*vol);
particle_force[i*3+2] = 0;
if(iel < 0){
continue;
}
int *tet = problem->mesh->elements+iel*4;
double C[] = {porosity[tet[0]],porosity[tet[1]],porosity[tet[2]],porosity[tet[3]]};
double U[] = {solution[tet[0]*4+0],solution[tet[1]*4+0],solution[tet[2]*4+0],solution[tet[3]*4+0]};
double V[] = {solution[tet[0]*4+1],solution[tet[1]*4+1],solution[tet[2]*4+1],solution[tet[3]*4+1]};
double W[] = {solution[tet[0]*4+2],solution[tet[1]*4+2],solution[tet[2]*4+2],solution[tet[3]*4+2]};
double P[] = {solution[tet[0]*4+3],solution[tet[1]*4+3],solution[tet[2]*4+3],solution[tet[3]*4+3]};
double vf[3] = {
phi[0]*U[0]+phi[1]*U[1]+phi[2]*U[2]+phi[3]*U[3],
phi[0]*V[0]+phi[1]*V[1]+phi[2]*V[2]+phi[3]*V[3],
phi[0]*W[0]+phi[1]*W[1]+phi[2]*W[2]+phi[3]*W[3]
};
double c = phi[0]*C[0]+phi[1]*C[1]+phi[2]*C[2]+phi[3]*C[3];
const double x[4][3] = {
{mesh->x[tet[0]*3+0],mesh->x[tet[0]*3+1],mesh->x[tet[0]*3+2]},
{mesh->x[tet[1]*3+0],mesh->x[tet[1]*3+1],mesh->x[tet[1]*3+2]},
{mesh->x[tet[2]*3+0],mesh->x[tet[2]*3+1],mesh->x[tet[2]*3+2]},
{mesh->x[tet[3]*3+0],mesh->x[tet[3]*3+1],mesh->x[tet[3]*3+2]}};
const double dxdxi[3][3] = {
{x[1][0]-x[0][0],x[2][0]-x[0][0],x[3][0]-x[0][0]},
{x[1][1]-x[0][1],x[2][1]-x[0][1],x[3][1]-x[0][1]},
{x[1][2]-x[0][2],x[2][2]-x[0][2],x[3][2]-x[0][2]}};
const double detj = dxdxi[0][0]*(dxdxi[1][1]*dxdxi[2][2]-dxdxi[2][1]*dxdxi[1][2])-dxdxi[1][0]*(dxdxi[0][1]*dxdxi[2][2]-dxdxi[2][1]*dxdxi[0][2])+dxdxi[2][0]*(dxdxi[0][1]*dxdxi[1][2]-dxdxi[1][1]*dxdxi[0][2]);
const double dxidx[3][3] = {
{(dxdxi[1][1]*dxdxi[2][2]-dxdxi[2][1]*dxdxi[1][2])/detj,-(dxdxi[0][1]*dxdxi[2][2]-dxdxi[2][1]*dxdxi[0][2])/detj, (dxdxi[0][1]*dxdxi[1][2]-dxdxi[1][1]*dxdxi[0][2])/detj},
{-(dxdxi[1][0]*dxdxi[2][2]-dxdxi[2][0]*dxdxi[1][2])/detj,(dxdxi[0][0]*dxdxi[2][2]-dxdxi[2][0]*dxdxi[0][2])/detj,-(dxdxi[0][0]*dxdxi[1][2]-dxdxi[1][0]*dxdxi[0][2])/detj },
{(dxdxi[1][0]*dxdxi[2][1]-dxdxi[2][0]*dxdxi[1][1])/detj,-(dxdxi[0][0]*dxdxi[2][1]-dxdxi[2][0]*dxdxi[0][1])/detj, (dxdxi[0][0]*dxdxi[1][1]-dxdxi[1][0]*dxdxi[0][1])/detj}
};
const double dphi[][3] = {
{-dxidx[0][0]-dxidx[1][0]-dxidx[2][0],-dxidx[0][1]-dxidx[1][1]-dxidx[2][1],-dxidx[0][2]-dxidx[1][2]-dxidx[2][2]},
{dxidx[0][0],dxidx[0][1],dxidx[0][2]},
{dxidx[1][0],dxidx[1][1],dxidx[1][2]},
{dxidx[2][0],dxidx[2][1],dxidx[2][2]}
};
double gradp[3] = {
dphi[0][0]*P[0]+dphi[1][0]*P[1]+dphi[2][0]*P[2]+dphi[3][0]*P[3],
dphi[0][1]*P[0]+dphi[1][1]*P[1]+dphi[2][1]*P[2]+dphi[3][1]*P[3],
dphi[0][2]*P[0]+dphi[1][2]*P[1]+dphi[2][2]*P[2]+dphi[3][2]*P[3]
};
double du[3] = {
problem->particle_velocity[i*3+0]-vf[0]/c,
problem->particle_velocity[i*3+1]-vf[1]/c,
problem->particle_velocity[i*3+2]-vf[2]/c
};
double rhop = mass/vol;
double d = 2*pow(3./4.*vol/M_PI,1./3.);
double drag[3];
particle_drag(du,problem->mu,problem->rho,d,c,drag, dt, mass, rhop,gradp);
particle_force[i*3+0] += -drag[0]-vol*gradp[0];
particle_force[i*3+1] += -drag[1]-vol*gradp[1];
particle_force[i*3+2] += -drag[2]-vol*gradp[2];
//printf("drag[0]=%e\n",drag[0]);
//printf("drag[1]=%e\n",drag[1]);
//printf("drag[2]=%e\n",drag[2]);
for (int iphi = 0; iphi < 4; ++iphi) {
problem->node_force[tet[iphi]*3+0] += drag[0]*phi[iphi];
problem->node_force[tet[iphi]*3+1] += drag[1]*phi[iphi];
problem->node_force[tet[iphi]*3+2] += drag[2]*phi[iphi];
}
}
for (int i = 0; i < mesh->n_nodes; ++i) {
problem->node_force[i*3+0] /= problem->node_volume[i];
problem->node_force[i*3+1] /= problem->node_volume[i];
problem->node_force[i*3+2] /= problem->node_volume[i];
}*/
} }
static void compute_node_force_implicit(FluidProblem *problem, double dt, double *all_local_vector, double *all_local_matrix) { static void compute_node_force_implicit(FluidProblem *problem, double dt, double *all_local_vector, double *all_local_matrix) {
...@@ -313,106 +137,70 @@ static void compute_node_force_implicit(FluidProblem *problem, double dt, double ...@@ -313,106 +137,70 @@ static void compute_node_force_implicit(FluidProblem *problem, double dt, double
Mesh *mesh = problem->mesh; Mesh *mesh = problem->mesh;
const double *solution = problem->solution; const double *solution = problem->solution;
for (int i = 0; i < problem->n_particles; ++i) { for (int i = 0; i < problem->n_particles; ++i) {
double *xi = problem->particle_uvw + 3*i;
double phi[4] = {1-xi[0]-xi[1]-xi[2],xi[0],xi[1],xi[2]};
int iel = problem->particle_element_id[i]; int iel = problem->particle_element_id[i];
double vol = problem->particle_volume[i]; double vol = problem->particle_volume[i];
double mass = problem->particle_mass[i]; double mass = problem->particle_mass[i];
if(iel < 0){ if(iel < 0){
continue; continue;
} }
double *xi = problem->particle_uvw + DIMENSION*i;
uint32_t *tet = problem->mesh->elements+iel*4; double phi[N_SF];
double C[] = {porosity[tet[0]],porosity[tet[1]],porosity[tet[2]],porosity[tet[3]]}; shape_functions(xi,phi);
double U[] = {solution[tet[0]*4+0],solution[tet[1]*4+0],solution[tet[2]*4+0],solution[tet[3]*4+0]}; uint32_t *el=problem->mesh->elements+iel*N_SF;
double V[] = {solution[tet[0]*4+1],solution[tet[1]*4+1],solution[tet[2]*4+1],solution[tet[3]*4+1]}; double *se = problem->solution_explicit;
double W[] = {solution[tet[0]*4+2],solution[tet[1]*4+2],solution[tet[2]*4+2],solution[tet[3]*4+2]}; double c=0, vf[DIMENSION]={0}, vfe[DIMENSION]={0}, p=0, gradp[DIMENSION]={0};
double P[] = {solution[tet[0]*4+3],solution[tet[1]*4+3],solution[tet[2]*4+3],solution[tet[3]*4+3]}; double dxdxi[DIMENSION][DIMENSION],dxidx[DIMENSION][DIMENSION];
double vf[3] = { for (int k=0; k<DIMENSION; ++k)
phi[0]*U[0]+phi[1]*U[1]+phi[2]*U[2]+phi[3]*U[3], for (int j=0; j<DIMENSION; ++j)
phi[0]*V[0]+phi[1]*V[1]+phi[2]*V[2]+phi[3]*V[3], dxdxi[k][j] = mesh->x[el[j+1]*3+k]-mesh->x[el[0]*3+k];
phi[0]*W[0]+phi[1]*W[1]+phi[2]*W[2]+phi[3]*W[3] invDD(dxdxi,dxidx);
}; double dphi[N_SF][DIMENSION];
double c = phi[0]*C[0]+phi[1]*C[1]+phi[2]*C[2]+phi[3]*C[3]; grad_shape_functions(dxidx, dphi);
const double x[4][3] = { for (int iphi=0; iphi < N_SF; ++iphi) {
{mesh->x[tet[0]*3+0],mesh->x[tet[0]*3+1],mesh->x[tet[0]*3+2]}, c += porosity[el[iphi]]*phi[iphi];
{mesh->x[tet[1]*3+0],mesh->x[tet[1]*3+1],mesh->x[tet[1]*3+2]}, for (int j=0; j < DIMENSION; ++j) {
{mesh->x[tet[2]*3+0],mesh->x[tet[2]*3+1],mesh->x[tet[2]*3+2]}, vf[j] += solution[el[iphi]*n_fields+j]*phi[iphi];
{mesh->x[tet[3]*3+0],mesh->x[tet[3]*3+1],mesh->x[tet[3]*3+2]}}; vfe[j] += se[el[iphi]*n_fields+j]*phi[iphi];
const double dxdxi[3][3] = { gradp[j] += dphi[iphi][j]*solution[el[iphi]*n_fields+DIMENSION];
{x[1][0]-x[0][0],x[2][0]-x[0][0],x[3][0]-x[0][0]}, }
{x[1][1]-x[0][1],x[2][1]-x[0][1],x[3][1]-x[0][1]}, p += solution[el[iphi]*n_fields+DIMENSION]*phi[iphi];
{x[1][2]-x[0][2],x[2][2]-x[0][2],x[3][2]-x[0][2]}}; }
const double detj = dxdxi[0][0]*(dxdxi[1][1]*dxdxi[2][2]-dxdxi[2][1]*dxdxi[1][2])-dxdxi[1][0]*(dxdxi[0][1]*dxdxi[2][2]-dxdxi[2][1]*dxdxi[0][2])+dxdxi[2][0]*(dxdxi[0][1]*dxdxi[1][2]-dxdxi[1][1]*dxdxi[0][2]); double du[DIMENSION], due[DIMENSION];
const double dxidx[3][3] = { for (int j = 0; j < DIMENSION; ++j) {
{(dxdxi[1][1]*dxdxi[2][2]-dxdxi[2][1]*dxdxi[1][2])/detj,-(dxdxi[0][1]*dxdxi[2][2]-dxdxi[2][1]*dxdxi[0][2])/detj, (dxdxi[0][1]*dxdxi[1][2]-dxdxi[1][1]*dxdxi[0][2])/detj}, du[j] = problem->particle_velocity[i*DIMENSION+j]-vf[j]/c;
{-(dxdxi[1][0]*dxdxi[2][2]-dxdxi[2][0]*dxdxi[1][2])/detj,(dxdxi[0][0]*dxdxi[2][2]-dxdxi[2][0]*dxdxi[0][2])/detj,-(dxdxi[0][0]*dxdxi[1][2]-dxdxi[1][0]*dxdxi[0][2])/detj }, due[j] = problem->particle_velocity[i*DIMENSION+j]-vfe[j]/c;
{(dxdxi[1][0]*dxdxi[2][1]-dxdxi[2][0]*dxdxi[1][1])/detj,-(dxdxi[0][0]*dxdxi[2][1]-dxdxi[2][0]*dxdxi[0][1])/detj, (dxdxi[0][0]*dxdxi[1][1]-dxdxi[1][0]*dxdxi[0][1])/detj} }
};
const double dphi[][3] = {
{-dxidx[0][0]-dxidx[1][0]-dxidx[2][0],-dxidx[0][1]-dxidx[1][1]-dxidx[2][1],-dxidx[0][2]-dxidx[1][2]-dxidx[2][2]},
{dxidx[0][0],dxidx[0][1],dxidx[0][2]},
{dxidx[1][0],dxidx[1][1],dxidx[1][2]},
{dxidx[2][0],dxidx[2][1],dxidx[2][2]}
};
double gradp[3] = {
dphi[0][0]*P[0]+dphi[1][0]*P[1]+dphi[2][0]*P[2]+dphi[3][0]*P[3],
dphi[0][1]*P[0]+dphi[1][1]*P[1]+dphi[2][1]*P[2]+dphi[3][1]*P[3],
dphi[0][2]*P[0]+dphi[1][2]*P[1]+dphi[2][2]*P[2]+dphi[3][2]*P[3]
};
double du[3] = {
problem->particle_velocity[i*3+0]-vf[0]/c,
problem->particle_velocity[i*3+1]-vf[1]/c,
problem->particle_velocity[i*3+2]-vf[2]/c
};
problem->particle_force[i*3+0] = -du[0]*mass/dt;
problem->particle_force[i*3+1] = -du[1]*mass/dt;
problem->particle_force[i*3+2] = -du[2]*mass/dt;
double rhop = mass/vol; double rhop = mass/vol;
double drho = rhop -problem->rho; double drho = rhop -problem->rho;
double drag[3] = { double gamma = particle_drag_coeff(due,problem->mu,problem->rho,vol,c);
mass*du[0]/dt - gradp[0]*vol, double g = problem->g*drho/rhop;
mass*du[1]/dt +problem->g*vol*drho - gradp[1]*vol, double fcoeff = mass/(mass+dt*gamma);
mass*du[2]/dt - gradp[2]*vol
}; 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;
double *local_vector = all_local_vector+iel*n_fields*N_SF; double *local_vector = all_local_vector+iel*n_fields*N_SF;
double *local_matrix = all_local_matrix+iel*n_fields*n_fields*N_SF*N_SF; double *local_matrix = all_local_matrix+iel*n_fields*n_fields*N_SF*N_SF;
for (int iphi = 0; iphi < 4; ++iphi) { for (int iphi = 0; iphi < N_SF; ++iphi) {
local_vector[iphi+N_SF*0] -= drag[0]*phi[iphi]; for (int d = 0; d < DIMENSION; ++d)
local_vector[iphi+N_SF*1] -= drag[1]*phi[iphi]; local_vector[iphi+N_SF*d] += fcoeff*gamma*(du[d]-dt/mass*vol*gradp[d])*phi[iphi];
local_vector[iphi+N_SF*2] -= drag[2]*phi[iphi]; local_vector[iphi+N_SF*1] += fcoeff*gamma*(-dt*g)*phi[iphi];
} }
for (int iphi = 0; iphi < 4; ++iphi) { for (int iphi = 0; iphi < N_SF; ++iphi) {
for (int jphi = 0; jphi < 4; ++jphi) { for (int jphi = 0; jphi < N_SF; ++jphi) {
int N2 = n_fields*N_SF; int N2 = n_fields*N_SF;
int IDX = (iphi*N2+jphi)*n_fields; int IDX = (iphi*N2+jphi)*n_fields;
#define LOCAL_MATRIX(a,b) local_matrix[IDX+N2*a+b] #define LOCAL_MATRIX(a,b) local_matrix[IDX+N2*a+b]
LOCAL_MATRIX(0,0) += phi[jphi]*phi[iphi]/(c*dt)*mass; double f = fcoeff*phi[iphi];
LOCAL_MATRIX(1,1) += phi[jphi]*phi[iphi]/(c*dt)*mass; for (int d = 0; d < DIMENSION; ++d){
LOCAL_MATRIX(2,2) += phi[jphi]*phi[iphi]/(c*dt)*mass; LOCAL_MATRIX(d,d) += -f/c*phi[jphi]*gamma;
LOCAL_MATRIX(d,DIMENSION) += -f*gamma*dt/mass*vol*dphi[jphi][d];
}
LOCAL_MATRIX(0,3) += phi[iphi]*vol*dphi[jphi][0];
LOCAL_MATRIX(1,3) += phi[iphi]*vol*dphi[jphi][1];
LOCAL_MATRIX(2,3) += phi[iphi]*vol*dphi[jphi][2];
} }
} }
/*
for (int iphi = 0; iphi < 4; ++iphi) {
problem->node_force[tet[iphi]*3+0] += drag[0]*phi[iphi];
problem->node_force[tet[iphi]*3+1] += drag[1]*phi[iphi];
problem->node_force[tet[iphi]*3+2] += drag[2]*phi[iphi];
}*/
} }
/*for (int i = 0; i < mesh->n_nodes; ++i) {
problem->node_force[i*3+0] /= problem->node_volume[i];
problem->node_force[i*3+1] /= problem->node_volume[i];
problem->node_force[i*3+2] /= problem->node_volume[i];
}*/
} }
#endif
static void fluid_problem_assemble_system(FluidProblem *problem, double *rhs, const double *solution_old, double dt) static void fluid_problem_assemble_system(FluidProblem *problem, double *rhs, const double *solution_old, double dt)
{ {
...@@ -587,6 +375,7 @@ int fluid_problem_implicit_euler(FluidProblem *problem, double dt) ...@@ -587,6 +375,7 @@ int fluid_problem_implicit_euler(FluidProblem *problem, double dt)
double *rhs = 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]; solution_old[i] = problem->solution[i];
problem->solution_explicit = solution_old;
double *dx = malloc(sizeof(double)*mesh->n_nodes*n_fields); double *dx = malloc(sizeof(double)*mesh->n_nodes*n_fields);
double norm0 = 0; double norm0 = 0;
int i; int i;
......
marblesbag/src/fluid_problem.h
View file @ cd61446c
...@@ -20,6 +20,7 @@ typedef struct { ...@@ -20,6 +20,7 @@ typedef struct {
double *porosity; double *porosity;
double *old_porosity; double *old_porosity;
double *solution; double *solution;
double *solution_explicit;
double *node_volume; double *node_volume;
int n_strong_boundaries; int n_strong_boundaries;
StrongBoundary *strong_boundaries; StrongBoundary *strong_boundaries;
......
testcases/drop-3d/metzger/drop.py
View file @ cd61446c
...@@ -47,19 +47,18 @@ rhop = 2450 ...@@ -47,19 +47,18 @@ rhop = 2450
nu = 1.17/1030. nu = 1.17/1030.
rho = 1030#rhop-drho/compacity rho = 1030#rhop-drho/compacity
V = 0.5 # todo : estimate V base on limit velocity V = 0.5 # todo : estimate V base on limit velocity
print('V',V)
tEnd = 1000 tEnd = 1000
#numerical parameters #numerical parameters
lcmin = 0.0004 # approx r*100 but should match the mesh size lcmin = 0.0004 # approx r*100 but should match the mesh size
dt = 5e-5*1000 dt = 5e-5*1000
alpha = 1e-3 alpha = 1e-3
epsilon = alpha*lcmin**2 /nu epsilon = alpha*lcmin**2/nu
print('epsilon',epsilon) print('epsilon',epsilon)
shutil.copy("mesh.msh", outputdir +"/mesh.msh") shutil.copy("mesh.msh", outputdir +"/mesh.msh")
#scontact2Interface.MeshLoader(p, "funnel.msh", ("Funnel", "Top", "Bottom", "Lateral")) #scontact2Interface.MeshLoader(p, "funnel.msh", ("Funnel", "Top", "Bottom", "Lateral"))
p.write(outputdir+"/part-00000") p.write_vtk(outputdir,0,0)
mu = nu*rho mu = nu*rho
genInitialPosition(r, R, rhop, compacity) genInitialPosition(r, R, rhop, compacity)
...@@ -124,13 +123,13 @@ while t < tEnd : ...@@ -124,13 +123,13 @@ while t < tEnd :
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity()) fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())
fluid.implicit_euler(dt) fluid.implicit_euler(dt)
forces = fluid.compute_node_force(dt) forces = fluid.compute_node_force(dt)
vn = p.velocity() + forces * dt / p.mass() vn = p.velocity() + forces*dt/p.mass()
vmax = np.max(np.hypot(vn[:, 0], vn[:, 1])) vmax = np.max(np.hypot(vn[:, 0], vn[:, 1]))
nsub = max(1, int(np.ceil((vmax * dt * 4)/min(p.r())))) nsub = max(1, int(np.ceil((vmax * dt * 4)/min(p.r()))))
print("NSUB", nsub,"VMAX",vmax, "VMAX * dt", vmax * dt, "r", min(p.r())) print("NSUB", nsub,"VMAX",vmax, "VMAX * dt", vmax * dt, "r", min(p.r()))
for i in range(nsub) : for i in range(nsub) :
p.iterate(dt/nsub, forces) p.iterate(dt/nsub, forces)
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity()) fluid.set_particles(p.mass(),p.volume(),p.position(),p.velocity())
t += dt t += dt
if ii %outf == 0 : if ii %outf == 0 :
ioutput = int(ii/outf) + 1 ioutput = int(ii/outf) + 1
......
Supports Markdown
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