Commit f0d9da13 authored by Hugo Ginestet's avatar Hugo Ginestet
Browse files

ale

parent 8a4fed08
......@@ -16,7 +16,7 @@
* 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
* 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/>.
*/
......@@ -174,7 +174,7 @@ static void node_force_volume(FluidProblem *problem, const double *solution_old,
free(sold);
}
static void f_boundary(WeakBoundary *wbnd, FluidProblem *problem,const double *n, double *f0,const double *s,const double *ds, const double *sold, const double *dsold, const double c, const double *dc, const double rho, const double mu, const double dt, int eid, const double *data, double *f00, double *f01)
static void f_boundary(WeakBoundary *wbnd, FluidProblem *problem,const double *n, double *f0,const double *s,const double *ds, const double *sold, const double *dsold,const double *u_mesh,const double c, const double *dc, const double rho, const double mu, const double dt, int eid, const double *data, double *f00, double *f01)
{
const int vid = wbnd->vid;
const int pid = wbnd->pid;
......@@ -317,7 +317,7 @@ void fluid_problem_compute_node_particle_force(FluidProblem *problem, double dt,
}
}
static void particle_force_f(FluidProblem *problem, double *f, double *dfdu, double *dfddp, double *sol, double *dsol, const double *solold, const double c, const double *dc, double a, double dt, int iel, int ip) {
static void particle_force_f(FluidProblem *problem, double *f, double *dfdu, double *dfddp, double *sol,double *dsol, const double *solold, const double c, const double *dc, double a, double dt, int iel, int ip) {
const double *contact = &problem->particle_contact[ip*D];
double u[D], uold[D], dp[D],G[D]={0};
for (int iD = 0; iD < D; ++iD) {
......@@ -361,15 +361,16 @@ static void particle_force_f(FluidProblem *problem, double *f, double *dfdu, dou
*dfddp = gammastar*dt/mass*vol-vol;
}
static void fluid_problem_f(const FluidProblem *problem, const double *sol, double *dsol, const double *solold, const double *dsolold, const double c, const double *dc, const double cold, const double rho, const double mu, double dt, int iel, double *f0, double *f1, double *f00, double *f10, double *f01, double *f11) {
static void fluid_problem_f(const FluidProblem *problem, const double *sol, double *dsol, const double *solold, const double *dsolold, const double c, const double *dc, const double cold, const double *u_mesh,const double *u_mesh_old, const double *du_mesh_vec, const double rho, const double mu, double dt, int iel, double *f0, double *f1, double *f00, double *f10, double *f01, double *f11) {
double p = sol[P];
double taup = problem->taup[iel];
double tauc = problem->tauc[iel];
const int n_fields = fluid_problem_n_fields(problem);
double u[D], uold[D], du[D][D], duold[D][D],dp[D];
double u[D], uold[D], du[D][D], duold[D][D],dp[D],du_mesh[D][D];
if(problem->stab_param != 0)
taup = tauc = 0;
double divu = 0;
double divu_mesh = 0;
for (int iD = 0; iD < D; ++iD) {
u[iD] = sol[U+iD];
uold[iD] = solold[U+iD];
......@@ -377,17 +378,21 @@ static void fluid_problem_f(const FluidProblem *problem, const double *sol, doub
for (int jD = 0; jD < D; ++jD) {
du[iD][jD] = dsol[(U+iD)*D+jD];
duold[iD][jD] = dsolold[(U+iD)*D+jD];
du_mesh[iD][jD] = du_mesh_vec[(U+iD)*D+jD];
}
divu += du[iD][iD];
divu_mesh -= u_mesh[iD]*dc[iD];
}
f0[P] = (c-cold)/dt;// + (sol[P]-solold[P])/dt*0.1;
//f0[P] = (c-cold)/dt ;// + (sol[P]-solold[P])/dt*0.1;
f0[P] = (c-cold)/dt + divu_mesh;
//f00[P*n_fields+P] = 1/dt*0.1;
double G[D] = {0};
G[1] = -problem->g;
double rhoreduced = (problem->reduced_gravity ? (rho-problem->rho[0]) : rho);
for (int i = 0; i < D; ++i) {
f0[U+i] =
rho*(u[i]-uold[i])/dt
rho*(u[i]-uold[i])/dt
+ dp[i]
+ G[i]*rhoreduced*c
+ problem->volume_drag*u[i]*mu;
......@@ -397,9 +402,12 @@ static void fluid_problem_f(const FluidProblem *problem, const double *sol, doub
// dj(uj ui /c) = uj/c dj(ui) + ui/c dj(uj) - uj ui /(c*c) dj(c)
// = ui/c * (dj(uj)- uj/c*dj(c)) + uj/c*dj(ui)
for (int j = 0; j < D; ++j) {
f0[U+i] += rho/c*(uold[j]/c*du[i][j] + u[i]*(duold[j][j]-uold[j]/c*dc[j]));
f0[U+i] += rho/c*((uold[j]-c*u_mesh[j])*du[i][j] + u[i]*(duold[j][j]-uold[j]/c*dc[j]));
//f0[U+i] += rho/c*(u[j]-c*u_mesh[j])*du[i][j];
f00[(U+i)*n_fields+U+i] += rho/c*(duold[j][j]-uold[j]/c*dc[j]);
f01[(U+i)*n_fields*D+(U+i)*D+j] += rho/c*uold[j]/c;
//f00[(U+i)*n_fields+U+i] += rho/c*du[i][j];
//f01[(U+i)*n_fields*D+(U+i)*D+j] += rho/c*(u[j]-c*u_mesh[j]);
f01[(U+i)*n_fields*D+(U+i)*D+j] += rho/c*(uold[j]-c*u_mesh[j]);
}
for (int j = 0; j < D; ++j) {
f1[(U+i)*D+j] = mu*(du[i][j]-u[i]*dc[j]/c+ du[j][i]-u[j]*dc[i]/c);
......@@ -408,6 +416,7 @@ static void fluid_problem_f(const FluidProblem *problem, const double *sol, doub
f11[((U+i)*D+j)*n_fields*D+(U+j)*D+i] += mu;
f11[((U+i)*D+j)*n_fields*D+(U+i)*D+j] += mu;
// SUPG
double supg = uold[j]*taup;
f1[(U+i)*D+j] += supg*f0[U+i];
f10[((U+i)*D+j)*n_fields+U+i] += supg*f00[(U+i)*n_fields+U+i];
......@@ -415,11 +424,14 @@ static void fluid_problem_f(const FluidProblem *problem, const double *sol, doub
for (int k = 0; k < D; ++k)
f11[((U+i)*D+j)*n_fields*D+(U+i)*D+k] += supg*f01[(U+i)*n_fields*D+(U+i)*D+k];
}
double lsic = tauc*rho;
f1[(U+i)*D+i] += ((c-cold)/dt+divu)*lsic;
f1[(U+i)*D+i] += ((c-cold)/dt +divu_mesh +divu)*lsic;
//f1[(U+i)*D+i] += ((c-cold)/dt+divu)*lsic;
for (int j = 0; j < D; ++j) {
f11[((U+i)*D+i)*(n_fields*D)+(U+j)*D+j] += lsic;
}
//f10[((U+i)*D+i)*(n_fields)+(U+i)*D] += dc[i]/c*lsic;
f1[P*D+i] = -u[i];
f10[(P*D+i)*n_fields+U+i] += -1;
// PSPG
......@@ -433,6 +445,7 @@ static void fluid_problem_f(const FluidProblem *problem, const double *sol, doub
}
}
static void compute_weak_boundary_conditions(FluidProblem *problem, const double *solution_old, double dt, double *all_local_vector, double *all_local_matrix)
{
const Mesh *mesh = problem->mesh;
......@@ -446,6 +459,7 @@ static void compute_weak_boundary_conditions(FluidProblem *problem, const double
double *f0 = malloc(sizeof(double)*n_fields);
double *f00 = malloc(sizeof(double)*n_fields*n_fields);
double *f01 = malloc(sizeof(double)*n_fields*n_fields*D);
double *u_mesh = malloc(sizeof(double)*D);
double i_bnd = 0;
double iv_bnd_up = 0;
double iv_bnd_down = 0;
......@@ -511,6 +525,8 @@ static void compute_weak_boundary_conditions(FluidProblem *problem, const double
sold[i] = 0;
f0[i] = 0;
}
for (int i=0; i<D;++i) u_mesh[i] = 0;
for (int i = 0; i < n_fields*n_fields; ++i) {
f00[i] = 0;
}
......@@ -530,6 +546,10 @@ static void compute_weak_boundary_conditions(FluidProblem *problem, const double
s[j] += phi[i]*dof;
sold[j] += phi[i]*dofold;
}
for (int j=0; j<D; ++j){
double dofmesh = problem->velocity_mesh[nodes[i]*D+j];
u_mesh[j] += phi[i]*dofmesh;
}
}
double rho, mu;
fluid_problem_interpolate_rho_and_nu(problem,a, &rho, &mu);
......@@ -546,7 +566,7 @@ static void compute_weak_boundary_conditions(FluidProblem *problem, const double
}
}
}
f_boundary(wbnd,problem,n,f0,s,ds,sold,dsold,c,dc,rho,mu,dt,iel,dataqp,f00,f01);
f_boundary(wbnd,problem,n,f0,s,ds,sold,dsold,u_mesh,c,dc,rho,mu,dt,iel,dataqp,f00,f01);
for (int ifield = 0; ifield < n_fields; ++ifield) {
for (int iphi = 0; iphi < D; ++iphi){
local_vector[cl[iphi]+N_SF*ifield] += phi[iphi]*f0[ifield]*jw;
......@@ -579,6 +599,7 @@ static void compute_weak_boundary_conditions(FluidProblem *problem, const double
free(dsold);
free(f0);
free(f00);
free(u_mesh);
}
double fluid_problem_a_integ_volume(FluidProblem *problem)
......@@ -731,6 +752,9 @@ static void fluid_problem_volume(FluidProblem *problem, const double *solution_o
double *sold = malloc(sizeof(double)*n_fields);
double *ds = malloc(sizeof(double)*n_fields*D);
double *dsold = malloc(sizeof(double)*n_fields*D);
double *du_mesh_vec = malloc(sizeof(double)*D*D);
double *u_mesh = malloc(sizeof(double)*D);
double *u_mesh_old = malloc(sizeof(double)*D);
double *f0 = malloc(sizeof(double)*n_fields);
double *f1 = malloc(sizeof(double)*n_fields*D);
double *f00 = malloc(sizeof(double)*n_fields*n_fields);
......@@ -740,7 +764,7 @@ static void fluid_problem_volume(FluidProblem *problem, const double *solution_o
problem->kinetic_energy = 0;
for (int iel=0; iel < mesh->n_elements; ++iel) {
const int *el = &mesh->elements[iel*N_N];
double dxidx[D][D], dphi[N_N][D];
double dxidx[D][D]={0}, dphi[N_N][D]={0};
const double detj = mesh_dxidx(mesh,iel,dxidx);
grad_shape_functions(dxidx, dphi);
double *local_vector = &all_local_vector[local_size*iel];
......@@ -751,6 +775,11 @@ static void fluid_problem_volume(FluidProblem *problem, const double *solution_o
dsold[i*D+j] = 0;
}
}
for (int i =0; i < D; ++i) {
for (int j = 0; j < D; ++j){
du_mesh_vec[i*D+j] = 0;
}
}
double dc[D] = {0}, da[D]={0};
for (int i = 0; i < N_SF; ++i) {
for (int j = 0; j < n_fields; ++j) {
......@@ -761,6 +790,13 @@ static void fluid_problem_volume(FluidProblem *problem, const double *solution_o
dsold[j*D+k] += dphi[i][k]*dofold;
}
}
for (int j = 0; j<D; ++j){
double dofmesh = problem->velocity_mesh[el[i]*D+j];
for (int k = 0; k < D; ++k) {
du_mesh_vec[j*D+k] += dphi[i][k]*dofmesh;
}
}
for (int k=0; k<D; ++k){
dc[k] += problem->porosity[el[i]]*dphi[i][k];
}
......@@ -770,6 +806,8 @@ static void fluid_problem_volume(FluidProblem *problem, const double *solution_o
}
}
}
for (int iqp = 0; iqp< N_QP; ++iqp) {
double phi[N_SF];
shape_functions(QP[iqp],phi);
......@@ -777,6 +815,11 @@ static void fluid_problem_volume(FluidProblem *problem, const double *solution_o
s[i] = 0;
sold[i] = 0;
}
for (int i=0; i<D;++i)
{
u_mesh[i]=0;
u_mesh_old[i]=0;
}
for (int i = 0; i < n_fields*n_fields; ++i) {
f00[i] = 0;
}
......@@ -790,8 +833,26 @@ static void fluid_problem_volume(FluidProblem *problem, const double *solution_o
double c = 0, a = 0;
double cold = 0, aold = 0;
for (int i = 0; i < N_SF; ++i) {
c += problem->porosity[el[i]]*phi[i];
for (int j = 0; j<D; ++j){
double dofmesh = problem->velocity_mesh[el[i]*D+j];
double dofmesh_old = problem->velocity_mesh_old[el[i]*D+j];
u_mesh[j] += dofmesh*phi[i];
u_mesh_old[j] += dofmesh_old*phi[i];
}
/*
if (problem->volume_old[el[i]]>1.e-30){
double rapport = problem->volume_new[el[i]]/problem->volume_old[el[i]];
//c += (1-rapport*(1-problem->porosity[el[i]]))*phi[i];
cold += rapport*problem->oldporosity[el[i]]*phi[i];
}
else{
cold += problem->oldporosity[el[i]]*phi[i];
//c += problem->porosity[el[i]]*phi[i];
}*/
cold += problem->oldporosity[el[i]]*phi[i];
//printf("\n %.5f \n\n\n",cold);
c += problem->porosity[el[i]]*phi[i];
for (int j = 0; j < n_fields; ++j) {
double dof = solution[el[i]*n_fields+j];
double dofold = solution_old[el[i]*n_fields+j];
......@@ -807,8 +868,7 @@ static void fluid_problem_volume(FluidProblem *problem, const double *solution_o
fluid_problem_interpolate_rho_and_nu(problem,a, &rho, &mu);
const double jw = QW[iqp]*detj;
problem->kinetic_energy += rho*sqrt(s[U]*s[U]+s[U+1]*s[U+1])*sqrt(s[U]*s[U]+s[U+1]*s[U+1])/2*jw;
fluid_problem_f(problem,s,ds,sold,dsold,c,dc,cold,rho,mu,dt,iel,f0,f1,f00,f10,f01,f11);
fluid_problem_f(problem,s,ds,sold,dsold,c,dc,cold,u_mesh,u_mesh_old,du_mesh_vec,rho,mu,dt,iel,f0,f1,f00,f10,f01,f11);
for (int ifield = 0; ifield < n_fields; ++ifield) {
for (int iphi = 0; iphi < N_SF; ++iphi){
local_vector[iphi+N_SF*ifield] += phi[iphi]*f0[ifield]*jw;
......@@ -842,6 +902,9 @@ static void fluid_problem_volume(FluidProblem *problem, const double *solution_o
free(ds);
free(dsold);
free(sold);
free(du_mesh_vec);
free(u_mesh);
free(u_mesh_old);
free(f0);
free(f1);
free(f00);
......@@ -1075,6 +1138,22 @@ static void apply_boundary_conditions(FluidProblem *problem)
free(v);
}
}
static void fluid_problem_compute_node_volume(FluidProblem *problem) {
free(problem->node_volume);
Mesh *mesh = problem->mesh;
problem->node_volume = (double*)malloc(sizeof(double)*mesh->n_nodes);
for (int i = 0; i < problem->mesh->n_nodes; ++i){
problem->node_volume[i] = 0;
}
for (int iel = 0; iel < mesh->n_elements; ++iel) {
double dxidx[D][D];
const int *el = &mesh->elements[iel*N_N];
const double detj = mesh_dxidx(mesh, iel, dxidx);
const double vol = node_volume_from_detj(detj);
for (int i = 0; i < N_N; ++i)
problem->node_volume[el[i]] += vol;
}
}
int fluid_problem_implicit_euler(FluidProblem *problem, double dt, double newton_atol, double newton_rtol, int newton_max_it)
......@@ -1085,9 +1164,26 @@ int fluid_problem_implicit_euler(FluidProblem *problem, double dt, double newton
int n_fields = fluid_problem_n_fields(problem);
apply_boundary_conditions(problem);
double *solution_old = (double*)malloc(mesh->n_nodes*n_fields*sizeof(double));
double *velocity_mesh_old = (double*)malloc(mesh->n_nodes*D*sizeof(double));
double *rhs = (double*)malloc(mesh->n_nodes*n_fields*sizeof(double));
double *volume_old = (double*)malloc(mesh->n_nodes*sizeof(double));
double *volume_new = (double*)malloc(mesh->n_nodes*sizeof(double));
for (int i=0; i<mesh->n_nodes*D; ++i)
velocity_mesh_old[i] = problem->velocity_mesh[i];
for (int i=0; i<mesh->n_nodes; ++i)
volume_old[i] = problem->node_volume[i];
//fluid_problem_compute_node_volume(problem);
for (int i=0; i<mesh->n_nodes; ++i)
volume_new[i] = problem->node_volume[i];
problem->volume_new = volume_new;
problem->volume_old = volume_old;
for (int i=0; i<mesh->n_nodes*n_fields; ++i)
solution_old[i] = problem->solution[i];
problem->velocity_mesh_old = velocity_mesh_old;
problem->solution_explicit = solution_old;
double *dx = (double*)malloc(sizeof(double)*mesh->n_nodes*n_fields);
double norm0 = 0;
......@@ -1097,6 +1193,7 @@ int fluid_problem_implicit_euler(FluidProblem *problem, double dt, double newton
double norm = 0;
hxtLinearSystemZeroMatrix(problem->linear_system);
for (int i=0; i<mesh->n_nodes*n_fields; ++i){
rhs[i] = 0;
}
......@@ -1116,34 +1213,23 @@ int fluid_problem_implicit_euler(FluidProblem *problem, double dt, double newton
for (int j=0; j<mesh->n_nodes*n_fields; ++j) {
problem->solution[j] -= dx[j];
}
node_force_volume(problem,problem->solution_explicit,dt,NULL,NULL);
break; /// equation is linear
node_force_volume(problem,problem->solution_explicit,dt,NULL,NULL);
//break; /// equation is linear
}
printf("total solve time : %g\n", totaltime);
free(dx);
free(rhs);
free(solution_old);
free(velocity_mesh_old);
free(volume_old);
free(volume_new);
return i;
}
static void fluid_problem_compute_node_volume(FluidProblem *problem) {
free(problem->node_volume);
Mesh *mesh = problem->mesh;
problem->node_volume = (double*)malloc(sizeof(double)*mesh->n_nodes);
for (int i = 0; i < problem->mesh->n_nodes; ++i){
problem->node_volume[i] = 0;
}
for (int iel = 0; iel < mesh->n_elements; ++iel) {
double dxidx[D][D];
const int *el = &mesh->elements[iel*N_N];
const double detj = mesh_dxidx(mesh, iel, dxidx);
const double vol = node_volume_from_detj(detj);
for (int i = 0; i < N_N; ++i)
problem->node_volume[el[i]] += vol;
}
}
static void fluid_problem_compute_porosity(FluidProblem *problem)
{
Mesh *mesh = problem->mesh;
......@@ -1210,6 +1296,7 @@ FluidProblem *fluid_problem_new(double g, int n_fluids, double *mu, double *rho,
problem->oldporosity = NULL;
problem->solution = NULL;
problem->concentration = NULL;
problem->velocity_mesh = NULL;
problem->linear_system = NULL;
problem->n_particles = 0;
problem->particle_uvw = NULL;
......@@ -1239,6 +1326,7 @@ void fluid_problem_free(FluidProblem *problem)
free(problem->rho);
free(problem->solution);
free(problem->concentration);
free(problem->velocity_mesh);
free(problem->porosity);
free(problem->oldporosity);
free(problem->node_volume);
......@@ -1504,12 +1592,16 @@ static void fluid_problem_set_mesh(FluidProblem *problem, Mesh *mesh) {
mesh_tree_free(problem->mesh_tree);
problem->mesh_tree = mesh_tree_create(mesh);
free(problem->porosity);
free(problem->velocity_mesh);
problem->porosity = (double*)malloc(mesh->n_nodes*sizeof(double));
problem->velocity_mesh = (double*)malloc(mesh->n_nodes*D*sizeof(double));
free(problem->oldporosity);
problem->oldporosity = (double*)malloc(mesh->n_nodes*sizeof(double));
for(int i = 0; i < mesh->n_nodes; ++i){
problem->porosity[i] = 1.;
problem->oldporosity[i] = 1.;
problem->velocity_mesh[i*D] = 0.;
problem->velocity_mesh[i*D+1] = 0.;
}
free(problem->solution);
int n_fields = fluid_problem_n_fields(problem);
......@@ -1667,7 +1759,7 @@ void fluid_problem_set_particles(FluidProblem *problem, int n, double *mass, dou
problem->particle_volume[i] = volume[i];
for (int k = 0; k < D; ++k) {
problem->particle_position[i*D+k] = position[i*D+k];
problem->particle_velocity[i*D+k] = velocity[i*D+k];
problem->particle_velocity[i*D+k] = velocity[i*D+k] ;
problem->particle_contact[i*D+k] = contact[i*D+k];
}
}
......@@ -1722,13 +1814,15 @@ void fluid_problem_set_strong_boundary(FluidProblem *problem, const char *tag, i
void fluid_problem_set_weak_boundary(FluidProblem *p, const char *tag, BoundaryType type, BoundaryCallback *cb, int vid, int pid, int cid, int aid)
{
int cpt = 0;
for (int i = 0; i < p->n_weak_boundaries; ++i) {
if (strcmp(p->weak_boundaries[i].tag, tag) == 0){
printf("Weak boundary is already defined for tag \"%s\".\n", tag);
exit(1);
//printf("Weak boundary is already defined for tag \"%s\".\n", tag);
//exit(1);
cpt++;
}
}
p->n_weak_boundaries++;
if(cpt == 0) p->n_weak_boundaries++;
p->weak_boundaries = realloc(p->weak_boundaries,sizeof(WeakBoundary)*p->n_weak_boundaries);
WeakBoundary *bnd = &p->weak_boundaries[p->n_weak_boundaries-1];
bnd->tag = strdup(tag);
......@@ -1755,6 +1849,11 @@ double *fluid_problem_porosity(const FluidProblem *p)
return p->porosity;
}
double *fluid_problem_velocity_mesh(const FluidProblem *p)
{
return p->velocity_mesh;
}
double *fluid_problem_stab_param(FluidProblem *problem)
{
return problem->taup;
......
......@@ -66,9 +66,14 @@ struct FluidProblem {
MeshBoundary *boundaries;
double *porosity;
double *oldporosity;
double *velocity_mesh;
double *concentration;
double *solution;
double *solution_explicit;
double *volume_new;
double *volume_old;
double *velocity_mesh_old;
double *velocity_old;
double *node_volume;
double *element_size;
double volume_drag;
......@@ -114,6 +119,7 @@ void fluid_problem_set_weak_boundary(FluidProblem *p, const char *tag, BoundaryT
void fluid_problem_set_strong_boundary(FluidProblem *p, const char *tag, int field, int row, BoundaryCallback *apply);
double *fluid_problem_old_porosity(const FluidProblem *p);
double *fluid_problem_porosity(const FluidProblem *p);
double *fluid_problem_velocity_mesh(const FluidProblem *p);
double *fluid_problem_concentration_dg(const FluidProblem *p);
double *fluid_problem_element_size(const FluidProblem *p);
......
......@@ -345,6 +345,10 @@ class FluidProblem :
def coordinates(self) :
"""Give access to the coordinate of the mesh nodes."""
return self._get_matrix("coordinates",self.n_nodes(), 3)
def velocity_mesh(self):
return self._get_matrix("velocity_mesh", self.n_nodes(), 2)
def elements(self):
"""Read-only access to the elements of the mesh."""
......
......@@ -20,5 +20,5 @@
# along with this program (see COPYING and COPYING.LESSER files). If not,
# see <http://www.gnu.org/licenses/>.
export PYTHONPATH=@PROJECT_BINARY_DIR@:$PYTHONPATH
export PYTHONPATH=/home/hugo/Documents/ens/mms/stage/MigFlow/python:$PYTHONPATH
python $*
......@@ -62,21 +62,21 @@ def genInitialPosition(filename, r, H, ly, lx, rhop) :
p.write_vtk(filename,0,0)
# Define output directory
outputdir = "output"
outputdir = "output2"
if not os.path.isdir(outputdir) :
os.makedirs(outputdir)
# Physical parameters
g = -9.81 # gravity
r = 1e-3 # particles radius
rhop = 1500 # particles density
rho = 1000 # fluid density
nu = 1e-6 # kinematic viscosity
g = -9.81 # gravity
r = 2.e-3 # particles radius
rhop = 2000 # particles density
rho = 1000 # fluid density
nu = 1.e-6 # kinematic viscosity
# Numerical parameters
outf = 5 # number of iterations between output files
dt = 0.25e-2 # time step
tEnd = 100 # final time
outf = 1 # number of iterations between output files
dt = 0.25e-2 # time step
tEnd = 10 # final time
# Geometrical parameters
ly = 5e-2 # particles area height
......@@ -119,7 +119,7 @@ fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity(),p.contact_f
#
while t < tEnd :
# Fluid solver
fluid.implicit_euler(dt)
fluid.implicit_euler(dt, newton_atol=5e-6, newton_rtol=1e-1, newton_max_it=35, reduced_gravity=0, stab_param=0.)
# Computation of the new velocities
forces = fluid.compute_node_force(dt)
......@@ -130,7 +130,7 @@ while t < tEnd :
print("NSUB", nsub,"VMAX",vmax, "VMAX * dt", vmax * dt, "r", min(p.r()))
# NLGS iterations
for i in range(nsub) :
tol = 1e-8 #numerical tolerance for particles intersection
tol = 1e-6 #numerical tolerance for particles intersection
p.iterate(dt/nsub, forces, tol)
# Localisation of the particles in the fluid
......
L = 0.2;
H = 0.6;
y = 0;
lc = 0.01;
lc = 0.04;
Point(1) = {-L, H, 0,lc};
Point(2) = {-L, -H, 0,lc};
......
This diff is collapsed.
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