Skip to content

GitLab

Commit 7c6b89b4 authored by Nathan Coppin's avatar Nathan Coppin
Browse files

close to ok

parent 2bab93d2
Pipeline #8637 passed with stages
in 3 minutes and 5 seconds
Hide whitespace changes
Inline Side-by-side
fluid/fluid_problem.c
View file @ 7c6b89b4
......@@ -113,9 +113,11 @@ 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];
}
......@@ -366,7 +368,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 dcdt, 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 cold, 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];
......@@ -388,7 +390,7 @@ static void fluid_problem_f(const FluidProblem *problem, const double *sol, doub
divu += du[iD][iD];
}
nold = sqrt(nold);
f0[P] = dcdt;// + (sol[P]-solold[P])/dt*0.1; // dcdt
f0[P] = (c-cold)/dt;// + (sol[P]-solold[P])/dt*0.1;
//f00[P*n_fields+P] = 1/dt*0.1;
double *g = problem->g;
double rhoreduced = (problem->reduced_gravity ? (rho-problem->rho[0]) : rho);
......@@ -431,7 +433,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] += (dcdt+divu)*lsic;//dcdt
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;
}
......@@ -810,10 +812,10 @@ static void fluid_problem_volume(FluidProblem *problem, const double *solution_o
f11[i] = 0;
}
double c = 0, a = 0;
double dcdt = 0, aold = 0;
double cold = 0, aold = 0;
for (int i = 0; i < N_SF; ++i) {
c += problem->porosity[el[i]]*phi[i];
dcdt += problem->dcdt[el[i]]*phi[i];
cold += problem->oldporosity[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];
......@@ -833,7 +835,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,dcdt,rho,mu,dt,iel,f0,f1,f00,f10,f01,f11);
fluid_problem_f(problem,s,ds,sold,dsold,c,dc,bf,cold,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;
......@@ -1223,10 +1225,9 @@ static void fluid_problem_compute_node_volume(FluidProblem *problem) {
}
}
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)
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)
{
for (int i = 0; i < mesh->n_nodes; ++i){
dcdt[i] = porosity[i];
porosity[i] = 0;
}
for (int i = 0; i < n_particles; ++i) {
......@@ -1244,7 +1245,6 @@ 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->dcdt = NULL;
problem->oldporosity = 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->dcdt);
free(problem->oldporosity);
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->dcdt);
problem->dcdt = (double*)malloc(mesh->n_nodes*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->dcdt[i] = 0.;
problem->oldporosity[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, double dt)
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)
{
struct timespec tic, toc;
fluid_problem_adapt_gen_mesh(problem, lcmax, lcmin, n_el, cmax, cmin);
......@@ -1697,8 +1697,10 @@ 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->dcdt, dt, 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->n_particles, problem->particle_position, problem->particle_volume, problem->particle_element_id, problem->particle_uvw);
free(old_particle_element_id);
free(old_particle_uvw);
......@@ -1713,7 +1715,7 @@ void fluid_problem_load_msh(FluidProblem *problem, const char *mesh_file_name) {
}
void fluid_problem_after_import(FluidProblem *problem, double dt)
void fluid_problem_after_import(FluidProblem *problem)
{
if(problem->mesh_tree)
mesh_tree_free(problem->mesh_tree);
......@@ -1723,10 +1725,10 @@ void fluid_problem_after_import(FluidProblem *problem, double dt)
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->dcdt, dt, 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->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, double dt)
void fluid_problem_set_particles(FluidProblem *problem, int n, double *mass, double *volume, double *state, double *contact, long int *elid, int init, int reload)
{
struct timespec tic, toc;
......@@ -1772,7 +1774,7 @@ 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->dcdt[i] = 0;
problem->oldporosity[i] = problem->porosity[i];
}
}
if (elid) {
......@@ -1780,12 +1782,13 @@ void fluid_problem_set_particles(FluidProblem *problem, int n, double *mass, dou
if(elid[i] == -1){
problem->particle_element_id[i]=elid[i];
}
}
}
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);
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);
if (init){
for (int i=0; i<problem->mesh->n_nodes; ++i){
problem->dcdt[i] = 0;
problem->oldporosity[i] = problem->porosity[i];
}
}
}
......@@ -1895,9 +1898,9 @@ double *fluid_problem_stab_param(FluidProblem *problem)
return problem->taup;
}
double *fluid_problem_dcdt(const FluidProblem *p)
double *fluid_problem_old_porosity(const FluidProblem *p)
{
return p->dcdt;
return p->oldporosity;
}
......
fluid/fluid_problem.h
View file @ 7c6b89b4
......@@ -66,7 +66,7 @@ struct FluidProblem {
MeshTree *mesh_tree;
MeshBoundary *boundaries;
double *porosity;
double *dcdt;
double *oldporosity;
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, double dt);
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_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, double dt);
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);
int *fluid_problem_particle_element_id(FluidProblem *problem);
int fluid_problem_n_particles(FluidProblem *problem);
void fluid_problem_after_import(FluidProblem *problem, double dt);
void fluid_problem_after_import(FluidProblem *problem);
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_dcdt(const FluidProblem *p);
double *fluid_problem_old_porosity(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);
......
fluid/fluid_problem_concentration.c
View file @ 7c6b89b4
......@@ -55,10 +55,12 @@ 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];
}
......
python/fluid.py
View file @ 7c6b89b4
......@@ -265,7 +265,7 @@ class FluidProblem :
("pressure",self.solution()[:,[self._dim]]),
("velocity",v),
("porosity",self.porosity()),
("dcdt",self.dcdt()),
("old_porosity",self.old_porosity()),
("grad",da)
]
cell_data.append(("curvature",self.curvature()))
......@@ -303,7 +303,7 @@ class FluidProblem :
if self._n_fluids == 2 :
self.concentration_dg()[:] = cdata["concentration"]
self.porosity()[:] = data["porosity"]
self.dcdt()[:] = data["dcdt"]
self.old_porosity()[:] = data["old_porosity"]
self._lib.fluid_problem_after_import(self._fp,c_double(dt))
self.sys = None
......@@ -407,7 +407,7 @@ class FluidProblem :
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))
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))
def _get_matrix(self, f_name, nrow, ncol,typec=c_double) :
......@@ -479,9 +479,9 @@ class FluidProblem :
"""Returns the porosity at previous time step"""
return self._get_matrix("curvature", self.n_elements(), 1)
def dcdt(self):
def old_porosity(self):
"""Returns the porosity time derivative at current time step"""
return self._get_matrix("dcdt", self.n_nodes(), 1)
return self._get_matrix("old_porosity", self.n_nodes(), 1)
def volume_flux(self, bnd_tag):
"""Computes the integral of the (outgoing) normal velocity through boundary with tag bnd_tag"""
......
python/time_integration.py
View file @ 7c6b89b4
......@@ -128,10 +128,10 @@ def predictor_corrector_iterate(fluid, particles, dt, min_nsub=1, contact_tol=1e
_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(), dt,elid=flags_inside)
dcdt1 = fluid.dcdt().copy()
dcdt1 = fluid.old_porosity().copy()
porosity1 = fluid.porosity().copy()
fluid.set_particles(particles.mass(), particles.volume(), particles.state(), particles.contact_forces(), dt)
fluid.dcdt()[:] = dcdt1
fluid.old_porosity()[:] = dcdt1 + fluid.porosity() - porosity1
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) :
......@@ -176,9 +176,9 @@ def particle_changed(fluid, p, dt) :
""" Modifies the porosity and fluid velocity fields when particles have been removed or added in the fluid domain
"""
c = fluid.porosity().copy()
dcdt0 = fluid.dcdt().copy()
dcdt0 = fluid.old_porosity().copy()
fluid.set_particles(p.mass(),p.volume(), p.state(), p.contact_forces(), dt)
c2 = fluid.porosity()
fluid.dcdt()[:] = dcdt0
fluid.old_porosity()[:] = dcdt0 + c2 -c
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