Skip to content

GitLab

Commit 3e767b64 authored by Matthieu Constant's avatar Matthieu Constant
Browse files

option stabilisation

parent 91e721f3
Hide whitespace changes
Inline Side-by-side
fluid/fluid_problem.c
View file @ 3e767b64
......@@ -495,7 +495,7 @@ static void compute_stab_parameters(FluidProblem *problem, double dt) {
}
}
static void fluid_problem_f(FluidProblem *problem, double *f0, double *f1, double *sol, double *dsol, const double *solold, const double c, const double *dc, const double cold, double dt, int iel, int reduced_gravity) {
static void fluid_problem_f(FluidProblem *problem, double *f0, double *f1, double *sol, double *dsol, const double *solold, const double c, const double *dc, const double cold, double dt, int iel, int reduced_gravity, double stab_param) {
double p = sol[P];
double taup = problem->taup[iel];
double tauc = problem->tauc[iel];
......@@ -562,10 +562,10 @@ static void fluid_problem_f(FluidProblem *problem, double *f0, double *f1, doubl
double R = dudt + u[i]*drhodt/rho + dp[i]/rho + u[i]*divu/c + utau[i]/c + uugradrho[i]/rho + (i==1 ? -(reduced_gravity ? 0 : 1)*c*problem->g : 0.);
f0[U+i] = rho*dudt + u[i]*drhodt + (i==1 ? -(reduced_gravity ? 0 : 1)*c*rho*problem->g : 0.);
for (int j = 0; j < D; ++j) {
f1[(U+i)*D+j] = -rho*u[i]*u[j]/c + mu*(tau[i][j]+tau[j][i]) + (i==j ? -p : 0) + R*u[j]*rho*taup;
f1[(U+i)*D+j] = -rho*u[i]*u[j]/c + mu*(tau[i][j]+tau[j][i]) + (i==j ? -p : 0) + (stab_param==0 ? 1 : 0)*R*u[j]*rho*taup;
}
f1[(U+i)*D+i] += divu*tauc*rho;
f1[P*D+i] = -u[i] + 0.025*dp[i];
f1[(U+i)*D+i] += (stab_param==0 ? 1 : 0)*divu*tauc*rho;
f1[P*D+i] = -u[i] + (stab_param==0 ? taup*R : stab_param*dp[i]);
if(problem->n_fluids == 2){
f1[P*D+i] += taua*(divu+(c-cold)/dt)*u[i];
f1[A*D+i] = -u[i]*a + taua*Ra*u[i] + taup*R*a + da[i]*(1e-7+extraa);
......@@ -735,7 +735,7 @@ static void compute_weak_boundary_conditions(FluidProblem *problem, double dt, d
}
static void fluid_problem_volume(FluidProblem *problem, const double *solution_old, double dt, double *all_local_vector, double *all_local_matrix, int reduced_gravity)
static void fluid_problem_volume(FluidProblem *problem, const double *solution_old, double dt, double *all_local_vector, double *all_local_matrix, int reduced_gravity, double stab_param)
{
const Mesh *mesh = problem->mesh;
const double *porosity = problem->porosity;
......@@ -791,7 +791,7 @@ static void fluid_problem_volume(FluidProblem *problem, const double *solution_o
}
}
const double jw = QW[iqp]*detj;
fluid_problem_f(problem,f0,f1,s,ds,sold,c,dc,cold,dt,iel, reduced_gravity);
fluid_problem_f(problem,f0,f1,s,ds,sold,c,dc,cold,dt,iel, reduced_gravity, stab_param);
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;
......@@ -803,12 +803,12 @@ static void fluid_problem_volume(FluidProblem *problem, const double *solution_o
for (int jfield = 0; jfield < n_fields; ++jfield) {
double store = s[jfield];
s[jfield] += deps;
fluid_problem_f(problem,g0,g1,s,ds,sold,c,dc,cold,dt,iel,reduced_gravity);
fluid_problem_f(problem,g0,g1,s,ds,sold,c,dc,cold,dt,iel,reduced_gravity,stab_param);
s[jfield] = store;
for (int jd =0; jd < D; ++jd){
store = ds[jfield*D+jd];
ds[jfield*D+jd] += deps;
fluid_problem_f(problem,h0+jd*n_fields,h1+jd*D*n_fields,s,ds,sold,c,dc,cold,dt,iel,reduced_gravity);
fluid_problem_f(problem,h0+jd*n_fields,h1+jd*D*n_fields,s,ds,sold,c,dc,cold,dt,iel,reduced_gravity,stab_param);
ds[jfield*D+jd] = store;
}
int N2 = n_fields*N_SF;
......@@ -850,7 +850,7 @@ static void fluid_problem_volume(FluidProblem *problem, const double *solution_o
free(h1);
}
static void fluid_problem_assemble_system(FluidProblem *problem, double *rhs, const double *solution_old, double dt, int reduced_gravity)
static void fluid_problem_assemble_system(FluidProblem *problem, double *rhs, const double *solution_old, double dt, int reduced_gravity, double stab_param)
{
HXTLinearSystem *lsys = problem->linear_system;
const Mesh *mesh = problem->mesh;
......@@ -884,7 +884,7 @@ static void fluid_problem_assemble_system(FluidProblem *problem, double *rhs, co
}
node_force_volume(problem, solution_old, dt, all_local_vector, all_local_matrix,reduced_gravity);
compute_weak_boundary_conditions(problem, dt, all_local_vector, all_local_matrix);
fluid_problem_volume(problem, solution_old, dt, all_local_vector, all_local_matrix,reduced_gravity);
fluid_problem_volume(problem, solution_old, dt, all_local_vector, all_local_matrix,reduced_gravity,stab_param);
for (int iel=0; iel < mesh->n_elements; ++iel) {
double *local_vector = &all_local_vector[local_size*iel];
double *local_matrix = &all_local_matrix[local_size*local_size*iel];
......@@ -938,7 +938,7 @@ static void apply_boundary_conditions(const Mesh *mesh, double *solution, int n_
}
}
int fluid_problem_implicit_euler(FluidProblem *problem, double dt, double newton_atol, double newton_rtol, int newton_max_it, int reduced_gravity)
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)
{
static double totaltime = 0;
struct timespec tic, toc;
......@@ -963,7 +963,7 @@ int fluid_problem_implicit_euler(FluidProblem *problem, double dt, double newton
//printf("solution=%g\n",problem->solution[i]);
}
//break;
fluid_problem_assemble_system(problem, rhs, solution_old, dt,reduced_gravity);
fluid_problem_assemble_system(problem, rhs, solution_old, dt,reduced_gravity,stab_param);
hxtLinearSystemGetRhsNorm(problem->linear_system,rhs,&norm);
printf("iter %i norm = %.16g\n", i, norm);
if (norm < newton_atol)
......@@ -1362,30 +1362,18 @@ void fluid_problem_adapt_mesh(FluidProblem *problem, double lcmax, double lcmin,
{
struct timespec tic, toc;
fluid_problem_adapt_gen_mesh(problem, lcmax, lcmin, n_el, cmax, cmin);
clock_get(&tic);
Mesh *new_mesh = mesh_load_msh("adapt/mesh.msh");
clock_get(&toc);
printf("Time mesh_load : %.6e \n",clock_diff(&tic,&toc));
int n_fields = fluid_problem_n_fields(problem);
double *new_solution = (double*)malloc(sizeof(double)*new_mesh->n_nodes*n_fields);
double *new_xi = (double*)malloc(sizeof(double)*new_mesh->n_nodes*D);
clock_get(&tic);
int *new_eid = (int*)malloc(sizeof(int)*new_mesh->n_nodes);
for (int i = 0; i < new_mesh->n_nodes; ++i)
new_eid[i] = -1;
clock_get(&toc);
printf("Time new_eid : %.6e \n",clock_diff(&tic,&toc));
clock_get(&tic);
double *newx = (double*)malloc(sizeof(double)*D*new_mesh->n_nodes);
for (int i = 0;i < new_mesh->n_nodes;++i)
for (int k = 0; k < D; ++k)
newx[i*D+k] = new_mesh->x[i*3+k];
clock_get(&toc);
printf("Time newx : %.6e \n",clock_diff(&tic,&toc));
clock_get(&tic);
mesh_tree_particle_to_mesh(problem->mesh_tree,new_mesh->n_nodes, newx, new_eid, new_xi);
clock_get(&toc);
printf("Time find new point in old mesh : %.6e \n",clock_diff(&tic,&toc));
for (int i = 0; i < new_mesh->n_nodes; ++i) {
unsigned int *el = problem->mesh->elements+new_eid[i]*N_N;
if(new_eid[i] < 0) {
......@@ -1473,11 +1461,7 @@ void fluid_problem_set_particles(FluidProblem *problem, int n, double *mass, dou
problem->particle_element_id[i]=elid[i];
}
}
clock_get(&tic);
mesh_tree_particle_to_mesh(problem->mesh_tree, n, position, problem->particle_element_id, problem->particle_uvw);
clock_get(&toc);
printf("Time mesh_particles_to_mesh : %.6e \n",clock_diff(&tic,&toc));
clock_get(&tic);
for (int i = 0; i < n; ++i) {
problem->particle_mass[i] = mass[i];
problem->particle_volume[i] = volume[i];
......@@ -1491,12 +1475,7 @@ void fluid_problem_set_particles(FluidProblem *problem, int n, double *mass, dou
problem->oldporosity[i] = problem->porosity[i];
}
}
clock_get(&toc);
printf("Time initialize particle_fluid : %.6e \n",clock_diff(&tic,&toc));
clock_get(&tic);
fluid_problem_compute_porosity(problem);
clock_get(&toc);
printf("Time compute_porosity : %.6e \n",clock_diff(&tic,&toc));
}
int *fluid_problem_particle_element_id(FluidProblem *problem)
......
fluid/fluid_problem.h
View file @ 3e767b64
......@@ -88,7 +88,7 @@ int fluid_problem_export_vtk(const FluidProblem *problem, const char *output_dir
int fluid_problem_import_vtk(FluidProblem *problem, const char *filename);
// complete force on the particle (including gravity)
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);
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_free(FluidProblem *problem);
FluidProblem *fluid_problem_new(double g, int n_fluids, double *mu, double *rho, double coeffStab, const char *petsc_solver_type);
......
python/fluid.py
View file @ 3e767b64
......@@ -167,7 +167,7 @@ class FluidProblem :
self._lib.fluid_problem_compute_node_particle_force(self._fp, c_double(dt), c_void_p(forces.ctypes.data))
return forces
def implicit_euler(self, dt, newton_atol=5e-6, newton_rtol=1e-5, newton_max_it=10, reduced_gravity=0) :
def implicit_euler(self, dt, newton_atol=5e-6, newton_rtol=1e-5, newton_max_it=10, reduced_gravity=0, stab_param=0.) :
"""Solve the fluid equations.
Keyword arguments:
......@@ -176,7 +176,7 @@ 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))
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))
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.
......
testcases/reduction2d/fd_depot-fluid.py
View file @ 3e767b64
......@@ -82,7 +82,7 @@ tic = time.time()
while t < tEnd :
#if ii%100==0 and ii != 0:
# fluid.adapt_mesh(0.1,10,8e-4)
fluid.implicit_euler(dt)
fluid.implicit_euler(dt, stab_param=1e-5)
forces = fluid.compute_node_force(dt)
vn = p.velocity() + forces * dt / p.mass()
vmax = np.max(np.hypot(vn[:, 0], vn[:, 1]))
......
testcases/depot-3d/smallDep3d/dep.py → testcases/smallDepot-3d/depot.py
View file @ 3e767b64
......@@ -20,6 +20,10 @@
# see <http://www.gnu.org/licenses/>.
#!/usr/bin/env python
# TESTCASE DESCRIPTION
# Deposit of 3D particles in a box full of fluid.
from migflow import fluid
from migflow import scontact
from migflow import lmgc90Interface
......@@ -31,8 +35,22 @@ import shutil
import random
def genInitialPosition(filename, rM, rm, ly, lxz, rhop) :
"""Set all the particles centre positions and create the particles objects to add in the computing structure.
For this testcase, radius are contained in a file radius.csv created with a special distribution
Keyword arguments:
filename -- name of the output file
rM -- maximum radius
rm -- minimum radius
ly -- particles area height
lxz -- particles area width
rhop -- particles density
"""
# Particles structure builder
p = scontact.ParticleProblem(3)
# Load the mesh.msh file specifying physical boundary names
p.load_msh_boundaries("mesh.msh", ["Top","Bottom","Z","X"])
# Positioning the particles on a regular grid
for x in np.arange(-lxz+rM, lxz-rM, 2*rM):
for y in np.arange(0.015+rM, 0.015+ly-rM, 2*rM):
for z in np.arange(-lxz+rM, lxz-rM, 2*rM):
......@@ -46,58 +64,58 @@ outputdir = "output"
if not os.path.isdir(outputdir) :
os.makedirs(outputdir)
#physical parameters
g = -9.81
rho = 1000
rhop = 2500
nu = 1e-6
V = 0.5 # todo : estimate V base on limit velocity
print('V',V)
tEnd = 5.
# Physical parameters
g = -9.81 # gravity
rho = 1000 # fluid density
rhop = 2500 # particle density
nu = 1e-6 # fluid viscosity
# r min
rm = 3e-4
# r max
rM = 5e-4
# H
ly = 0.005
# L | P
lxz = 0.025
# Numerical parameters
tEnd = 5. # final time
dt = 1e-4 # time step
outf = 10 # number of iterations between output files
genInitialPosition(outputdir, rM, rm, ly, lxz, rhop)
# Geometrical parameters
rm = 3e-4 # minimum radius
rM = 5e-4 # maximum radius
ly = 0.005 # particles area height
lxz = 0.025 # particles area width
#
# PARTICLE PROBLEM
#
# Particles structure creation and setting particles properties
genInitialPosition(outputdir, rM, rm, ly, lxz, rhop)
friction = 0.1
#lmgc90Interface.scontactTolmgc90(outputdir,3,0,friction,assume_box=True)
#p = lmgc90Interface.ParticleProblem(3)
p = scontact.ParticleProblem(3)
p.read_vtk(outputdir,0)
lmgc90Interface.scontactTolmgc90(outputdir,3,0,friction,assume_box=True)
p = lmgc90Interface.ParticleProblem(3)
#numerical parameters
dt = 1e-4 # time step
# Initial time and iteration
t = 0
ii = 0
#
# FLUID PROBLEM
#
fluid = fluid.FluidProblem(3,g,[nu*rho],[rho])
#Set the mesh geometry for the fluid computation
fluid.load_msh("mesh.msh")
fluid.set_strong_boundary("Top",1,0.)
fluid.set_strong_boundary("X",0,0.)
fluid.set_strong_boundary("Z",2,0.)
fluid.set_strong_boundary("Bottom",1,0.)
fluid.set_strong_boundary("Top",3,0.)
#fluid.solution()[:,3] = rho*g*(fluid.coordinates()[:,1]-0.025)
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())
fluid.export_vtk(outputdir,0,0)
tic = time.time()
#Computation loop
outf = 1
t = 0.
ii = 0
tic = time.clock()
#
# COMPUTATION LOOP
#
while t < tEnd :
#Fluid solver
fluid.implicit_euler(dt,reduced_gravity=1)
fluid.implicit_euler(dt,reduced_gravity=1, stab_param=0.025)
forces = fluid.compute_node_force(dt)
#Computation of the new velocities
vn = p.velocity() + forces * dt / p.mass()
......@@ -119,4 +137,4 @@ while t < tEnd :
p.write_vtk(outputdir, ioutput, t)
fluid.export_vtk(outputdir, t, ioutput)
ii += 1
print("%i : %.2g/%.2g (cpu %.6g)" % (ii, t, tEnd, time.clock() - tic))
print("%i : %.2g/%.2g (cpu %.6g)" % (ii, t, tEnd, time.time() - tic))
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