Skip to content
GitLab
Commit e8491b0c authored by Nathan Coppin's avatar Nathan Coppin
Browse files

adapting patankar

parent 12c190ff
Pipeline #9379 failed with stages
in 59 seconds
Hide whitespace changes
Inline Side-by-side
fluid/fluid_problem.c
View file @ e8491b0c
......@@ -133,7 +133,7 @@ void fluid_problem_node_force_volume(FluidProblem *problem, const double *soluti
dc[k] += problem->porosity[el[i]]*dphi[i][k];
}
}
double f[D],dfdu,dfddp;
double f[D],dfdu[D][D],dfddp[D][D];
particle_force_f(problem,f,&dfdu,&dfddp,s,ds,sold,c,dc,a,dt,iel,ip);
if (!local_vector)
continue;
......@@ -144,9 +144,9 @@ void fluid_problem_node_force_volume(FluidProblem *problem, const double *soluti
for (int iphi = 0; iphi < N_SF; ++iphi){
local_vector[iphi+N_SF*(U+iD)] += phi[iphi]*f[iD];
local_vector[iphi+N_SF*P] += pspg*dphi[iphi][iD]*f[iD];
for (int jD = 0; jD < D; ++jD) {
double supg = drag_in_stab?sold[U+jD]*problem->taup[iel]:0;
local_vector[iphi+N_SF*(U+iD)] += supg*dphi[iphi][jD]*f[iD];
for (int kD = 0; kD < D; ++kD) {
double supg = drag_in_stab?sold[U+kD]*problem->taup[iel]:0;
local_vector[iphi+N_SF*(U+iD)] += supg*dphi[iphi][kD]*f[iD];
}
}
}
......@@ -155,16 +155,18 @@ void fluid_problem_node_force_volume(FluidProblem *problem, const double *soluti
for (int iD = 0; iD < D; ++iD){
for (int iphi = 0; iphi < N_SF; ++iphi){
for (int jphi = 0; jphi < N_SF; ++jphi){
LOCAL_MATRIX(iphi,jphi,U+iD,U+iD) += phi[iphi]*phi[jphi]*dfdu;
LOCAL_MATRIX(iphi,jphi,U+iD,P) += phi[iphi]*dphi[jphi][iD]*dfddp;
LOCAL_MATRIX(iphi,jphi,P,U+iD) += pspg*dphi[iphi][iD]*phi[jphi]*dfdu;
LOCAL_MATRIX(iphi,jphi,P,P) += pspg*dphi[iphi][iD]*dphi[jphi][iD]*dfddp;
for (int jD = 0; jD < D; ++jD) {
double supg = drag_in_stab?sold[U+jD]*problem->taup[iel]:0;
LOCAL_MATRIX(iphi,jphi,U+iD,U+iD) += supg*dphi[iphi][jD]*phi[jphi]*dfdu;
LOCAL_MATRIX(iphi,jphi,U+iD,P) += supg*dphi[iphi][jD]*dphi[jphi][iD]*dfddp;
for(int jD = 0; jD < D; jD++){
LOCAL_MATRIX(iphi,jphi,U+iD,U+jD) += phi[iphi]*phi[jphi]*dfdu[iD][jD];
LOCAL_MATRIX(iphi,jphi,U+iD,P) += phi[iphi]*dphi[jphi][jD]*dfddp[iD][jD];
LOCAL_MATRIX(iphi,jphi,P,U+iD) += pspg*dphi[iphi][iD]*phi[jphi]*dfdu[iD][jD];
LOCAL_MATRIX(iphi,jphi,P,P) += pspg*dphi[iphi][iD]*dphi[jphi][jD]*dfddp[iD][jD];
for (int kD = 0; kD < D; ++kD) {
double supg = drag_in_stab?sold[U+kD]*problem->taup[iel]:0;
LOCAL_MATRIX(iphi,jphi,U+iD,U+jD) += supg*dphi[iphi][kD]*phi[jphi]*dfdu[iD][jD];
LOCAL_MATRIX(iphi,jphi,U+iD,P) += supg*dphi[iphi][kD]*dphi[jphi][jD]*dfddp[iD][jD];
}
}
}
}
}
}
......@@ -345,13 +347,13 @@ static void particle_force_f(FluidProblem *problem, double *f, double *dfdu, dou
for (int j = 0; j < D; ++j) {
Due[j] = problem->particle_velocity[ip*D+j]-uold[j]/c;
}
double mass = problem->particle_mass[ip];
//double mass = problem->particle_mass[ip];
double vol = problem->particle_volume[ip];
double rhop = mass/vol;
/*double rhop = mass/vol;
double massreduced = mass;
if(problem->reduced_gravity) {
massreduced = (rhop-problem->rho[0])*problem->particle_volume[ip];
}
}*/
double gamma;
if(problem->n_fluids == 2){
......@@ -363,15 +365,41 @@ static void particle_force_f(FluidProblem *problem, double *f, double *dfdu, dou
gamma = particle_drag_coeff(Due,mu,rho,vol,c);
}
gamma = gamma*problem->drag_coeff_factor;
double gammastar = gamma/(1+dt/mass*gamma);
double mat[D][D] = {[0][0] = 1, [1][1] = 1};
#if DIMENSION==3
mat[2][2] = 1;
#endif
double wdp[D] = {0.};
double wcontact[D] = {0.};
double upstarstar[D];
for(int i = 0;i<D;i++){
for(int j = 0;j<D;j++){
mat[i][j] += gamma*dt*problem->particle_delassus[ip*D*D + i*D + j];
wdp[i] += problem->particle_delassus[ip*D*D + i*D + j]*dp[j];
wcontact[i] += problem->particle_delassus[ip*D*D + i*D + j]*contact[j];
}
upstarstar[i] = up[i]+dt*(problem->g[i] + wcontact[i] -vol*wdp[i]);
}
double invmat[D][D];
double lol = invDD(mat,invmat);
for (int i = 0; i < D; ++i) {
double upstar = up[i]+dt/mass*(contact[i]+massreduced*problem->g[i]-vol*dp[i]);
double drag = gammastar*(upstar-u[i]/c);
problem->particle_force[ip*D+i] = -drag-vol*dp[i];
f[U+i] = -drag;//-vol*dp[i];
f[U+i] = 0;
for(int j = 0;j<D;j++){
f[U+i] += -gamma*invmat[i][j]*(upstarstar[j]-u[j]/c);
}
problem->particle_force[ip*D+i] = f[U+i]-vol*dp[i];
}
for(int i = 0; i<D;i++){
for(int j = 0;j<D;j++){
dfdu[i*D + j] = gamma*invmat[i][j]/c;
double invmatw = 0;
for(int k = 0;k<D;k++){
invmatw += invmat[i][k]*problem->particle_delassus[ip*D*D + k*D + j];
}
dfddp[i*D +j] = -vol*gamma*dt*invmatw;
}
}
*dfdu = gammastar/c;
*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* mesh_velocity, 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) {
......@@ -1333,7 +1361,7 @@ FluidProblem *fluid_problem_new(double *g, int n_fluids, double *mu, double *rho
problem->particle_uvw = NULL;
problem->particle_element_id = NULL;
problem->particle_position = NULL;
problem->particle_mass = NULL;
problem->particle_delassus = NULL;
problem->particle_force = NULL;
problem->particle_volume = NULL;
problem->particle_velocity = NULL;
......@@ -1366,7 +1394,7 @@ void fluid_problem_free(FluidProblem *problem)
free(problem->particle_uvw);
free(problem->particle_element_id);
free(problem->particle_position);
free(problem->particle_mass);
free(problem->particle_delassus);
free(problem->boundary_force);
free(problem->particle_force);
free(problem->particle_volume);
......@@ -1755,13 +1783,13 @@ void fluid_problem_set_coordinates(FluidProblem *problem, double *x) {
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, NULL);
}
void fluid_problem_set_particles(FluidProblem *problem, int n, double *mass, double *volume, double *position, double *velocity, double *contact) {
void fluid_problem_set_particles(FluidProblem *problem, int n, double *delassus, double *volume, double *position, double *velocity, double *contact) {
if (problem->n_particles != n) {
problem->n_particles = n;
free(problem->particle_uvw);
free(problem->particle_element_id);
free(problem->particle_position);
free(problem->particle_mass);
free(problem->particle_delassus);
free(problem->particle_force);
free(problem->particle_volume);
free(problem->particle_velocity);
......@@ -1769,19 +1797,20 @@ void fluid_problem_set_particles(FluidProblem *problem, int n, double *mass, dou
problem->particle_uvw = (double*)malloc(sizeof(double)*D*n);
problem->particle_element_id = (int*)malloc(sizeof(int)*n);
problem->particle_position = (double*)malloc(sizeof(double)*D*n);
problem->particle_mass = (double*)malloc(sizeof(double)*n);
problem->particle_delassus = (double*)malloc(sizeof(double)*n*D*D);
problem->particle_force = (double*)malloc(sizeof(double)*n*D);
problem->particle_volume = (double*)malloc(sizeof(double)*n);
problem->particle_velocity = (double*)malloc(sizeof(double)*n*D);
problem->particle_contact = (double*)malloc(sizeof(double)*n*D);
}
for (int i = 0; i < n; ++i) {
problem->particle_mass[i] = mass[i];
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_contact[i*D+k] = contact[i*D+k];
for(int l = 0; l < D; l++)
problem->particle_delassus[i*D*D+k*D+l] = delassus[i*D*D+k*D+l];
}
problem->particle_element_id[i]=-1;
}
......
fluid/fluid_problem.h
View file @ e8491b0c
......@@ -85,7 +85,7 @@ struct FluidProblem {
double *particle_volume;
double *particle_velocity;
double *particle_contact;
double *particle_mass;
double *particle_delassus;
double *particle_force;
double *grad_a_cg;
double *bulk_force;
......
python/fluid.py
View file @ e8491b0c
......@@ -509,18 +509,18 @@ 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, position, velocity, contact):
def set_particles(self, delassus, volume, position, velocity, contact):
"""Set location of the grains in the mesh and compute the porosity in each cell.
Keyword arguments:
mass -- List of particles mass
delassus -- List of particles delassus operators
volume -- List of particles volume
position -- List of particles centre positions
velocity -- List of particles velocity
contact -- List of particles contact resultant forces
"""
self.n_particles = mass.shape[0]
self._lib.fluid_problem_set_particles(self._fp,c_int(mass.shape[0]),_np2c(mass),_np2c(volume),_np2c(position),_np2c(velocity),_np2c(contact))
self.n_particles = delassus.shape[0]
self._lib.fluid_problem_set_particles(self._fp,c_int(delassus.shape[0]),_np2c(delassus),_np2c(volume),_np2c(position),_np2c(velocity),_np2c(contact))
def move_particles(self, position, velocity, contact):
"""Set location of the grains in the mesh and compute the porosity in each cell.
......
python/scontact.py
View file @ e8491b0c
......@@ -153,6 +153,12 @@ class ParticleProblem :
pvel.flags.writeable = False
return pvel
def delassus(self):
delassus = np.zeros((self.n_particles(),self._dim,self._dim))
self._lib.particle_problem_get_delassus(self._p,_np2c(delassus))
return delassus
def body_velocity(self):
return self._get_array("body", self._bodytype)["velocity"]
......@@ -163,6 +169,9 @@ class ParticleProblem :
else:
return self._get_array("body", self._bodytype)["omega"]
def body_theta(self) :
return self._get_array("body", self._bodytype)["theta"][:,None]
def save_state(self) :
self._saved_velocity = np.copy(self.body_velocity())
self._saved_position = np.copy(self.body_position())
......@@ -170,13 +179,13 @@ class ParticleProblem :
if self.dim() == 3 :
self._saved_triangles = np.copy(self.triangles())
self._saved_omega = np.copy(self.body_omega())
self._saved_theta = np.copy(self._bodies()["theta"])
self._saved_theta = np.copy(self.body_theta())
def restore_state(self) :
self.body_velocity()[:] = self._saved_velocity
self.body_position()[:] = self._saved_position
self.body_omega()[:] = self._saved_omega
self._bodies()["theta"] = _self.saved_theta
self.body_theta()[:] = self._saved_theta
self.segments()[:] = self._saved_segments
if self.dim() == 3 :
self.triangles()[:] = self._saved_triangles
......
python/time_integration.py
View file @ e8491b0c
......@@ -39,7 +39,7 @@ def _advance_particles(particles, f, dt, min_nsub,contact_tol,max_nsub=None,afte
f = np.zeros_like(v)
# Estimation of the solid sub-time steps
if particles.r() is not None and particles.r().size != 0:
vn = v + f*dt / particles.mass()
vn = v + f*dt * np.amin(particles.body_invert_mass())
vmax = np.max(np.linalg.norm(vn,axis=1))
nsub = max(min_nsub, int(np.ceil((vmax * dt * 8)/min(particles.r()))))
else:
......@@ -89,13 +89,6 @@ def predictor_corrector_iterate(fluid, particles, dt, min_nsub=1, contact_tol=1e
external_particles_forces = np.zeros_like(particles.velocity())
# Copy the fluid solution of the previous time step before computing the prediction
sf0 = np.copy(fluid.solution())
# Copy the solid solution of the previous time step before computing the prediction
x0 = np.copy(particles.position())
cp0 = np.copy(particles.contact_forces())
seg0 = np.copy(particles.segments())
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
......@@ -124,7 +117,7 @@ def predictor_corrector_iterate(fluid, particles, dt, min_nsub=1, contact_tol=1e
particles.restore_state()
_advance_particles(particles,(alpha*f0+(1-alpha)*f1),dt,min_nsub,contact_tol,after_sub_iter=after_sub_iter,max_nsub=max_nsub)
# Give to the fluid the solid information
fluid.move_particles(particles.position(), particles.velocity(), particles.contact_forces())
fluid.move_particles(particles.position(), particles.velocity(), np.zeros_like(particles.velocity()))
def iterate(fluid, particles, dt, min_nsub=1, contact_tol=1e-8, external_particles_forces=None, fixed_grains=False, after_sub_iter=None,max_nsub=None) :
......
scontact/scontact.c
View file @ e8491b0c
......@@ -186,10 +186,33 @@ static void body_get_point(const Body *b, const double *xlocal, double *xglobal)
xglobal[0] = b->position[0] + cost*xlocal[0] - sint*xlocal[1];
xglobal[1] = b->position[1] + sint*xlocal[0] + cost*xlocal[1];
#else
// TODO
// TODO fake here
xglobal = {0,0,0};
#endif
}
void particle_problem_get_delassus(ParticleProblem *p, double *delassus){
for (int i = 0; i<vectorSize(p->particles); i++){
double r[DIMENSION] = {0};
double w[DIMENSION][DIMENSION] = {0};
#if DIMENSION==2
Body *body = &p->bodies[p->particles[i].body];
body_get_point(body,&p->particles[i].x,r);
r[0] -= body->position[0];
r[1] -= body->position[1];
w[0][1] = w[1][0] = -r[0]*r[1]*body->invertI;
w[0][0] = body->invertM + r[1]*r[1]*body->invertI;
w[1][1] = body->invertM + r[0]*r[0]*body->invertI;
#else
//TODO
#endif
for(int j = 0;j<DIMENSION;j++)
for(int k = 0;k<DIMENSION;k++)
delassus[i*DIMENSION*DIMENSION + j*DIMENSION + k] = w[j][k];
}
}
void particle_problem_get_particles_position(ParticleProblem *p, double *x){
for(int i = 0;i<vectorSize(p->particles);i++){
const Particle *p0 = &p->particles[i];;
......
testcases/depot-2d/depot.py
View file @ e8491b0c
......@@ -59,7 +59,7 @@ def genInitialPosition(filename, r, H, ly, lx, rhop) :
y = y.flat
# Add a grain at each centre position
for i in range(len(x)) :
p.add_particle((x[i], y[i]), r, r**2 * np.pi * rhop)
p.add_particle(np.array([x[i], y[i]]), r, r**2 * np.pi * rhop)
p.write_vtk(filename,0,0)
# Define output directory
......@@ -91,9 +91,8 @@ H = 0.6 # domain height
genInitialPosition(outputdir, r, H, ly, lx, rhop)
p = scontact.ParticleProblem(2)
p.read_vtk(outputdir,0)
print(p.body_position())
print("r = %g, m = %g\n" % (p.r()[0], p.mass()[0]))
print("RHOP = %g" % rhop)
# Initial time and iteration
t = 0
......@@ -109,7 +108,7 @@ fluid.set_wall_boundary("Bottom")
fluid.set_wall_boundary("Lateral")
fluid.set_wall_boundary("Top",pressure=0)
# Set location of the particles in the mesh and compute the porosity in each computation cell
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity(), p.contact_forces())
fluid.set_particles(p.delassus(), p.volume(), p.position(), p.velocity(), np.zeros_like(p.velocity()))
fluid.write_vtk(outputdir,0,0)
p.write_vtk(outputdir, 0, 0)
......@@ -118,7 +117,7 @@ tic = time.time()
# COMPUTATION LOOP
#
while t < tEnd :
time_integration.iterate(fluid, p, dt, min_nsub=5, external_particles_forces=g*p.mass())
time_integration.iterate(fluid, p, dt, min_nsub=5, external_particles_forces=g*p.r()**2*rhop*np.pi)
t += dt
# Output files writting
......
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