Commit 474031a2 authored by Nathan Coppin's avatar Nathan Coppin
Browse files

allowing entering particles to enter+dcdt instead of oldporos

parent cadabba3
Pipeline #8532 failed with stages
in 6 minutes and 50 seconds
......@@ -113,11 +113,9 @@ void fluid_problem_node_force_volume(FluidProblem *problem, const double *soluti
}
double c = 0;
double a = 0;
double cold = 0;
double dc[D] = {0};
for (int i = 0; i < N_SF; ++i) {
c += problem->porosity[el[i]]*phi[i];
cold += problem->oldporosity[el[i]]*phi[i];
if (problem->n_fluids == 2) {
a += problem->concentration[iel*N_N+i]*phi[i];
}
......@@ -368,7 +366,7 @@ 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 *bf, 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 *bf, const double dcdt, const double rho, const double mu, double dt, int iel, double *f0, double *f1, double *f00, double *f10, double *f01, double *f11) {//rempalcer cold par dcdt
double p = sol[P];
double taup = problem->taup[iel];
double tauc = problem->tauc[iel];
......@@ -390,7 +388,7 @@ static void fluid_problem_f(const FluidProblem *problem, const double *sol, doub
divu += du[iD][iD];
}
nold = sqrt(nold);
f0[P] = (c-cold)/dt;// + (sol[P]-solold[P])/dt*0.1;
f0[P] = dcdt;// + (sol[P]-solold[P])/dt*0.1; // dcdt
//f00[P*n_fields+P] = 1/dt*0.1;
double *g = problem->g;
double rhoreduced = (problem->reduced_gravity ? (rho-problem->rho[0]) : rho);
......@@ -433,7 +431,7 @@ static void fluid_problem_f(const FluidProblem *problem, const double *sol, doub
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] += (dcdt+divu)*lsic;//dcdt
for (int j = 0; j < D; ++j) {
f11[((U+i)*D+i)*(n_fields*D)+(U+j)*D+j] += lsic;
}
......@@ -812,10 +810,10 @@ static void fluid_problem_volume(FluidProblem *problem, const double *solution_o
f11[i] = 0;
}
double c = 0, a = 0;
double cold = 0, aold = 0;
double dcdt = 0, aold = 0; //cold par dcdt
for (int i = 0; i < N_SF; ++i) {
c += problem->porosity[el[i]]*phi[i];
cold += problem->oldporosity[el[i]]*phi[i];
dcdt += problem->dcdt[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];
......@@ -835,7 +833,7 @@ static void fluid_problem_volume(FluidProblem *problem, const double *solution_o
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,bf,cold,rho,mu,dt,iel,f0,f1,f00,f10,f01,f11);
fluid_problem_f(problem,s,ds,sold,dsold,c,dc,bf,dcdt,rho,mu,dt,iel,f0,f1,f00,f10,f01,f11);//dcdt a la palce de cold
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;
......@@ -1225,9 +1223,10 @@ static void fluid_problem_compute_node_volume(FluidProblem *problem) {
}
}
static void compute_porosity(Mesh *mesh, double *node_volume, double *porosity, int n_particles, double *particle_position, double *particle_volume, int *particle_element_id, double *particle_uvw)
static void compute_porosity(Mesh *mesh, double *node_volume, double *porosity, double *dcdt, double dt, int n_particles, double *particle_position, double *particle_volume, int *particle_element_id, double *particle_uvw)
{
for (int i = 0; i < mesh->n_nodes; ++i){
dcdt[i] = porosity[i];
porosity[i] = 0;
}
for (int i = 0; i < n_particles; ++i) {
......@@ -1245,6 +1244,7 @@ static void compute_porosity(Mesh *mesh, double *node_volume, double *porosity,
else{
porosity[i] = 1-porosity[i]/node_volume[i];
}
dcdt[i] = (porosity[i]-dcdt[i])/dt;
}
}
......@@ -1289,7 +1289,7 @@ FluidProblem *fluid_problem_new(double *g, int n_fluids, double *mu, double *rho
problem->tauc = NULL;
problem->element_size = NULL;
problem->porosity = NULL;
problem->oldporosity = NULL;
problem->dcdt = NULL;
problem->solution = NULL;
problem->concentration = NULL;
problem->n_particles = 0;
......@@ -1323,7 +1323,7 @@ void fluid_problem_free(FluidProblem *problem)
free(problem->concentration);
free(problem->porosity);
free(problem->bulk_force);
free(problem->oldporosity);
free(problem->dcdt);
free(problem->node_volume);
free(problem->particle_uvw);
free(problem->particle_element_id);
......@@ -1599,11 +1599,11 @@ static void fluid_problem_set_mesh(FluidProblem *problem, Mesh *mesh) {
problem->grad_a_cg = malloc(sizeof(double)*mesh->n_nodes*D);
problem->all_kappa = malloc(sizeof(double)*mesh->n_elements);
free(problem->oldporosity);
problem->oldporosity = (double*)malloc(mesh->n_nodes*sizeof(double));
free(problem->dcdt);
problem->dcdt = (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->dcdt[i] = 1.;
}
free(problem->boundary_force);
......@@ -1633,7 +1633,7 @@ double * fluid_problem_bulk_force(FluidProblem *problem) {
}
void fluid_problem_adapt_mesh(FluidProblem *problem, double lcmax, double lcmin, double n_el, double cmax, double cmin, int old_n_particles, double *old_particle_position, double *old_particle_volume)
void fluid_problem_adapt_mesh(FluidProblem *problem, double lcmax, double lcmin, double n_el, double cmax, double cmin, int old_n_particles, double *old_particle_position, double *old_particle_volume, double dt)
{
struct timespec tic, toc;
fluid_problem_adapt_gen_mesh(problem, lcmax, lcmin, n_el, cmax, cmin);
......@@ -1697,10 +1697,8 @@ void fluid_problem_adapt_mesh(FluidProblem *problem, double lcmax, double lcmin,
}
double *old_particle_uvw = malloc(sizeof(double)*D*old_n_particles);
mesh_tree_particle_to_mesh(problem->mesh_tree, old_n_particles, old_particle_position, old_particle_element_id, old_particle_uvw);
compute_porosity(problem->mesh, problem->node_volume, problem->oldporosity, old_n_particles, old_particle_position, old_particle_volume, old_particle_element_id, old_particle_uvw);
mesh_tree_particle_to_mesh(problem->mesh_tree, problem->n_particles, problem->particle_position, problem->particle_element_id, problem->particle_uvw);
compute_porosity(problem->mesh, problem->node_volume, problem->porosity, problem->n_particles, problem->particle_position, problem->particle_volume, problem->particle_element_id, problem->particle_uvw);
compute_porosity(problem->mesh, problem->node_volume, problem->porosity, problem->dcdt, dt, problem->n_particles, problem->particle_position, problem->particle_volume, problem->particle_element_id, problem->particle_uvw);
free(old_particle_element_id);
free(old_particle_uvw);
......@@ -1715,7 +1713,7 @@ void fluid_problem_load_msh(FluidProblem *problem, const char *mesh_file_name) {
}
void fluid_problem_after_import(FluidProblem *problem)
void fluid_problem_after_import(FluidProblem *problem, double dt)
{
if(problem->mesh_tree)
mesh_tree_free(problem->mesh_tree);
......@@ -1725,10 +1723,10 @@ void fluid_problem_after_import(FluidProblem *problem)
problem->particle_element_id[i] = -1;
mesh_tree_particle_to_mesh(problem->mesh_tree, problem->n_particles, problem->particle_position, problem->particle_element_id, problem->particle_uvw);
fluid_problem_compute_node_volume(problem);
compute_porosity(problem->mesh, problem->node_volume, problem->porosity, problem->n_particles, problem->particle_position, problem->particle_volume, problem->particle_element_id, problem->particle_uvw);
compute_porosity(problem->mesh, problem->node_volume, problem->porosity, problem->dcdt, dt, problem->n_particles, problem->particle_position, problem->particle_volume, problem->particle_element_id, problem->particle_uvw);
}
void fluid_problem_set_particles(FluidProblem *problem, int n, double *mass, double *volume, double *state, double *contact, long int *elid, int init, int reload)
void fluid_problem_set_particles(FluidProblem *problem, int n, double *mass, double *volume, double *state, double *contact, long int *elid, int init, int reload, double dt)
{
struct timespec tic, toc;
......@@ -1774,13 +1772,20 @@ void fluid_problem_set_particles(FluidProblem *problem, int n, double *mass, dou
mesh_tree_particle_to_mesh(problem->mesh_tree, n, problem->particle_position, problem->particle_element_id, problem->particle_uvw);
if (!reload){
for (int i=0; i<problem->mesh->n_nodes; ++i){
problem->oldporosity[i] = problem->porosity[i];
problem->dcdt[i] = 0;
}
}
if (elid) {
for (int i = 0; i < n; ++i) {
if(elid[i] == -1){
problem->particle_element_id[i]=elid[i];
}
}
}
compute_porosity(problem->mesh, problem->node_volume, problem->porosity, problem->n_particles, problem->particle_position, problem->particle_volume, problem->particle_element_id, problem->particle_uvw);
compute_porosity(problem->mesh, problem->node_volume, problem->porosity, problem->dcdt, dt, problem->n_particles, problem->particle_position, problem->particle_volume, problem->particle_element_id, problem->particle_uvw);
if (init){
for (int i=0; i<problem->mesh->n_nodes; ++i){
problem->oldporosity[i] = problem->porosity[i];
problem->dcdt[i] = 0;
}
}
}
......@@ -1890,9 +1895,9 @@ double *fluid_problem_stab_param(FluidProblem *problem)
return problem->taup;
}
double *fluid_problem_old_porosity(const FluidProblem *p)
double *fluid_problem_dcdt(const FluidProblem *p)
{
return p->oldporosity;
return p->dcdt;
}
......
......@@ -66,7 +66,7 @@ struct FluidProblem {
MeshTree *mesh_tree;
MeshBoundary *boundaries;
double *porosity;
double *oldporosity;
double *dcdt;
double *concentration;
double *solution;
double *node_volume;
......@@ -104,21 +104,21 @@ struct FluidProblem {
// complete force on the particle (including gravity)
void fluid_problem_compute_node_particle_force(FluidProblem *problem, double dt, double *particle_force);
void fluid_problem_set_particles(FluidProblem *problem, int n, double *mass, double *volume, double *state, double *contact, long int *elid, int init, int reload);
void fluid_problem_set_particles(FluidProblem *problem, int n, double *mass, double *volume, double *state, double *contact, long int *elid, int init, int reload, double dt);
void fluid_problem_set_bulk_force(FluidProblem *problem, double *force);
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 sigma, double coeffStab, double volume_drag, double quadratic_drag, int drag_in_stab, double drag_coeff_factor, double ip_factor,int temporal, int advection);
void fluid_problem_load_msh(FluidProblem *problem, const char *mesh_file_name);
void fluid_problem_adapt_mesh(FluidProblem *problem, double lcmax, double lcmin, double n_el, double cmax, double cmin, int old_n_particles, double *old_particle_position, double *old_particle_volume);
void fluid_problem_adapt_mesh(FluidProblem *problem, double lcmax, double lcmin, double n_el, double cmax, double cmin, int old_n_particles, double *old_particle_position, double *old_particle_volume, double dt);
int *fluid_problem_particle_element_id(FluidProblem *problem);
int fluid_problem_n_particles(FluidProblem *problem);
void fluid_problem_after_import(FluidProblem *problem);
void fluid_problem_after_import(FluidProblem *problem, double dt);
int fluid_problem_n_fields(const FluidProblem *p);
double *fluid_problem_solution(const FluidProblem *p);
void fluid_problem_set_weak_boundary(FluidProblem *p, const char *tag, BoundaryType type, BoundaryCallback *cb, int vid, int pid, int cid, int aid);
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_dcdt(const FluidProblem *p);
double *fluid_problem_porosity(const FluidProblem *p);
double *fluid_problem_concentration_dg(const FluidProblem *p);
double *fluid_problem_element_size(const FluidProblem *p);
......
......@@ -55,12 +55,10 @@ void fluid_problem_advance_concentration(FluidProblem *problem, double dt) {
shape_functions(QP[iqp],phi);
double u[D] = {0};
double c = 0;
double cold = 0;
double a = 0;
for (int i = 0; i < N_SF; ++i) {
a += phi[i]*problem->concentration[iel*N_N+i];
c += phi[i]*problem->porosity[el[i]];
cold += phi[i]*problem->oldporosity[el[i]];
for (int id = 0; id < D; ++id) {
u[id] += phi[i]*velocity[el[i]*D+id];
}
......
......@@ -217,7 +217,7 @@ class FluidProblem :
self.strong_cb_ref.append(bndcb)
self._lib.fluid_problem_set_strong_boundary(self._fp, tag.encode(), c_int(field), c_int(row), bndcb)
def adapt_mesh(self, lcmax, lcmin, n_el, old_n_particles, old_particle_position, old_particle_volume, cmax=1, cmin=1) :
def adapt_mesh(self, lcmax, lcmin, n_el, old_n_particles, old_particle_position, old_particle_volume, dt, cmax=1, cmin=1) :
"""Adapts the mesh.
Keyword arguments:
......@@ -227,10 +227,11 @@ class FluidProblem :
old_n_particles -- Number of particles at the previous time step
old_particle_position -- Position of the particles at the previous time step
old_particle_volume -- Volume of the particles at the previous time step
dt -- time step
cmax -- Optional argument to weigh maximum gradient used in the adaptation formula
cmin -- Optional argument to weigh minimum gradient used in the adaptation formula
"""
self._lib.fluid_problem_adapt_mesh(self._fp, c_double(lcmax), c_double(lcmin), c_double(n_el), c_double(cmax), c_double(cmin), c_int(old_n_particles), _np2c(old_particle_position), _np2c(old_particle_volume))
self._lib.fluid_problem_adapt_mesh(self._fp, c_double(lcmax), c_double(lcmin), c_double(n_el), c_double(cmax), c_double(cmin), c_int(old_n_particles), _np2c(old_particle_position), _np2c(old_particle_volume), c_double(dt))
self.sys = None
def _mesh_boundaries(self) :
......@@ -264,7 +265,7 @@ class FluidProblem :
("pressure",self.solution()[:,[self._dim]]),
("velocity",v),
("porosity",self.porosity()),
("old_porosity",self.old_porosity()),
("dcdt",self.dcdt()),
("grad",da)
]
cell_data.append(("curvature",self.curvature()))
......@@ -277,7 +278,7 @@ class FluidProblem :
offsets = np.cumsum(np.repeat([self._dim+1],connectivity.shape[0])).astype(np.int32)
VTK.write(output_dir+"/fluid",it,t,(types,connectivity,offsets),self.coordinates(),data,field_data,cell_data)
def import_vtk(self, filename) :
def import_vtk(self, filename, dt = 1) :
"""Reads output file to reload computation.
Keyword arguments:
......@@ -302,8 +303,8 @@ class FluidProblem :
if self._n_fluids == 2 :
self.concentration_dg()[:] = cdata["concentration"]
self.porosity()[:] = data["porosity"]
self.old_porosity()[:] = data["old_porosity"]
self._lib.fluid_problem_after_import(self._fp)
self.dcdt()[:] = data["dcdt"]
self._lib.fluid_problem_after_import(self._fp,c_double(dt))
self.sys = None
def compute_node_force(self, dt) :
......@@ -390,7 +391,7 @@ class FluidProblem :
for i in range(nsub) :
self._lib.fluid_problem_advance_concentration(self._fp,c_double(dt/nsub))
def set_particles(self, mass, volume, state, contact, init=False, reload = False):
def set_particles(self, mass, volume, state, contact, dt, elid = None, init=False, reload = False):
"""Set location of the grains in the mesh and compute the porosity in each cell.
Keyword arguments:
......@@ -398,11 +399,15 @@ class FluidProblem :
volume -- List of particles volume
position -- List of particles centre positions
velocity -- List of particles velocity
init -- if True sets the old_porosity to zero (has to be used at the first call)
dt -- time step
init -- if True sets the porosity derivative to zero (has to be used at the first call)
reload -- Optional arguments specifying if the simulation is reloaded or if it is the first iteration
"""
self.n_particles = mass.shape[0]
self._lib.fluid_problem_set_particles(self._fp,c_int(mass.shape[0]),_np2c(mass),_np2c(volume),_np2c(state),_np2c(contact),None,c_int(init),c_int(reload))
elemid = None
if elid is not None:
elemid = _np2c(elid,dtype=int)
self._lib.fluid_problem_set_particles(self._fp,c_int(mass.shape[0]),_np2c(mass),_np2c(volume),_np2c(state),_np2c(contact),elemid,c_int(init),c_int(reload),c_double(dt))
def _get_matrix(self, f_name, nrow, ncol,typec=c_double) :
......@@ -474,9 +479,9 @@ class FluidProblem :
"""Returns the porosity at previous time step"""
return self._get_matrix("curvature", self.n_elements(), 1)
def old_porosity(self):
"""Returns the porosity at previous time step"""
return self._get_matrix("old_porosity", self.n_nodes(), 1)
def dcdt(self):
"""Returns the porosity time derivative at current time step"""
return self._get_matrix("dcdt", self.n_nodes(), 1)
def volume_flux(self, bnd_tag):
"""Computes the integral of the (outgoing) normal velocity through boundary with tag bnd_tag"""
......
......@@ -78,6 +78,7 @@ def predictor_corrector_iterate(fluid, particles, dt, min_nsub=1, contact_tol=1e
after_sub_iter -- callback to execute once a sub iteration has been made
max_split -- maximum number of times the time step can be further split if conergence is not reached
"""
flags_inside = fluid.particle_element_id().copy()
state0 = particles.state()
if external_particles_forces is None:
external_particles_forces = np.zeros_like(state0.v)
......@@ -89,8 +90,6 @@ def predictor_corrector_iterate(fluid, particles, dt, min_nsub=1, contact_tol=1e
disk0 = np.copy(particles.disks())
if particles.dim() == 3 :
tri0 = np.copy(particles.triangles())
# Predictor
#
# Compute the fluid solution and keep a copy of this solution and the forces applied by the fluid on the grains
......@@ -104,7 +103,7 @@ def predictor_corrector_iterate(fluid, particles, dt, min_nsub=1, contact_tol=1e
particles.contact_forces()[:] = cp0
# Keep same contact forces and positions while giving to the fluid the predicted solid velocities to compute the correction
fluid.set_particles(particles.mass(), particles.volume(), particles.state(), particles.contact_forces(), reload=True)
fluid.set_particles(particles.mass(), particles.volume(), particles.state(), particles.contact_forces(), dt, reload=True)
# Corrector
#
......@@ -127,7 +126,12 @@ def predictor_corrector_iterate(fluid, particles, dt, min_nsub=1, contact_tol=1e
particles.triangles()[:] = tri0
_advance_particles(particles,(alpha*f0+(1-alpha)*f1),dt,min_nsub,contact_tol,after_sub_iter=after_sub_iter,max_split=max_split)
# Give to the fluid the solid information
fluid.set_particles(particles.mass(), particles.volume(), particles.state(), particles.contact_forces())
fluid.set_particles(particles.mass(), particles.volume(), particles.state(), particles.contact_forces(), dt,elid=flags_inside)
dcdt1 = fluid.dcdt().copy()
porosity1 = fluid.porosity().copy()
fluid.set_particles(particles.mass(), particles.volume(), particles.state(), particles.contact_forces(), dt)
fluid.dcdt()[:] = dcdt1
fluid.velocity()[:] *= fluid.porosity()/porosity1
def iterate(fluid, particles, dt, min_nsub=1, contact_tol=1e-8, external_particles_forces=None, fixed_grains=False, after_sub_iter=None,max_split=1) :
"""Migflow solver: the solver type depends on the fluid and particles given as arguments.
......@@ -163,17 +167,17 @@ def iterate(fluid, particles, dt, min_nsub=1, contact_tol=1e-8, external_particl
if fluid is not None and particles is not None and particles.r() is not None and particles.r().size !=0:
if fixed_grains:
fluid.implicit_euler(dt)
fluid.set_particles(particles.mass(), particles.volume(), particles.state(),-fluid.compute_node_force(dt))
fluid.set_particles(particles.mass(), particles.volume(), particles.state(),-fluid.compute_node_force(dt), dt)
else:
predictor_corrector_iterate(fluid,particles,dt,min_nsub,contact_tol,external_particles_forces,after_sub_iter=after_sub_iter,max_split=max_split)
def particle_changed(fluid,p) :
def particle_changed(fluid, p, dt) :
""" Modifies the porosity and fluid velocity fields when particles have been removed or added in the fluid domain
"""
co = fluid.old_porosity().copy()
c = fluid.porosity().copy()
fluid.set_particles(p.mass(),p.volume(), p.state(), p.contact_forces())
dcdt0 = fluid.dcdt().copy()
fluid.set_particles(p.mass(),p.volume(), p.state(), p.contact_forces(), dt)
c2 = fluid.porosity()
fluid.old_porosity()[:] = co+c2-c
fluid.dcdt()[:] = dcdt0
fluid.velocity()[:] *= c2/c
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