Commit 9a3dc018 authored by Michel Henry's avatar Michel Henry
Browse files

wip

parent 695429b5
Pipeline #9672 failed with stages
in 4 minutes and 42 seconds
......@@ -116,101 +116,101 @@ void fluid_problem_solution_at_particles(FluidProblem *problem, double *v_p, dou
}
}
// void fluid_problem_node_force_volume(FluidProblem *problem, const double *solution_old, double dt, double *all_local_vector, double *all_local_matrix)
// {
// const Mesh *mesh = problem->mesh;
// const double *porosity = problem->porosity;
// const double *solution = problem->solution;
// int drag_in_stab = problem->drag_in_stab;
// int n_fields = fluid_problem_n_fields(problem);
// size_t local_size = N_SF*n_fields;
// double *s = malloc(sizeof(double)*n_fields);
// double *sold = malloc(sizeof(double)*n_fields);
// double *ds = malloc(sizeof(double)*n_fields*D);
// for (int ip=0; ip < problem->n_particles; ++ip) {
// int iel = problem->particle_element_id[ip];
// if(iel < 0){
// for (int d = 0; d < D; ++d)
// problem->particle_force[ip*D+d] = 0;
// continue;
// }
// double *xi = problem->particle_uvw + D*ip;
// double phi[N_SF];
// shape_functions(xi,phi);
// const int *el = &mesh->elements[iel*N_N];
// double dxidx[D][D], dphi[N_SF][D];
// double detj = mesh_dxidx(mesh, iel, dxidx);
// grad_shape_functions(dxidx, dphi);
// double *local_vector = all_local_vector ? &all_local_vector[local_size*iel] : NULL;
// double *local_matrix = all_local_matrix ? &all_local_matrix[local_size*local_size*iel] : NULL;
// for (int i = 0; i < n_fields; ++i) {
// s[i] = 0;
// sold[i] = 0;
// for (int j = 0; j < D; ++j) {
// ds[i*D+j] = 0;
// }
// }
// 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];
// }
// 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];
// s[j] += phi[i]*dof;
// sold[j] += phi[i]*dofold;
// for (int k = 0; k < D; ++k)
// ds[j*D+k] += dphi[i][k]*dof;
// }
// for (int k=0; k<D; ++k){
// dc[k] += problem->porosity[el[i]]*dphi[i][k];
// }
// }
// double f[D],dfdu,dfddp;
// particle_force_f(problem,f,&dfdu,&dfddp,s,ds,sold,c,dc,a,dt,iel,ip);
// if (!local_vector)
// continue;
// double rho,nu;
// fluid_problem_interpolate_rho_and_nu(problem,a, &rho,&nu);
// double pspg = drag_in_stab?problem->taup[iel]/rho:0;
// for (int iD = 0; iD < D; ++iD) {
// 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];
// }
// }
// }
// if (!local_matrix)
// continue;
// 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;
// }
// }
// }
// }
// }
// free(s);
// free(ds);
// free(sold);
// }
void fluid_problem_node_force_volume(FluidProblem *problem, const double *solution_old, double dt, double *all_local_vector, double *all_local_matrix)
{
const Mesh *mesh = problem->mesh;
const double *porosity = problem->porosity;
const double *solution = problem->solution;
int drag_in_stab = problem->drag_in_stab;
int n_fields = fluid_problem_n_fields(problem);
size_t local_size = N_SF*n_fields;
double *s = malloc(sizeof(double)*n_fields);
double *sold = malloc(sizeof(double)*n_fields);
double *ds = malloc(sizeof(double)*n_fields*D);
for (int ip=0; ip < problem->n_particles; ++ip) {
int iel = problem->particle_element_id[ip];
if(iel < 0){
for (int d = 0; d < D; ++d)
problem->particle_force[ip*D+d] = 0;
continue;
}
double *xi = problem->particle_uvw + D*ip;
double phi[N_SF];
shape_functions(xi,phi);
const int *el = &mesh->elements[iel*N_N];
double dxidx[D][D], dphi[N_SF][D];
double detj = mesh_dxidx(mesh, iel, dxidx);
grad_shape_functions(dxidx, dphi);
double *local_vector = all_local_vector ? &all_local_vector[local_size*iel] : NULL;
double *local_matrix = all_local_matrix ? &all_local_matrix[local_size*local_size*iel] : NULL;
for (int i = 0; i < n_fields; ++i) {
s[i] = 0;
sold[i] = 0;
for (int j = 0; j < D; ++j) {
ds[i*D+j] = 0;
}
}
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];
}
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];
s[j] += phi[i]*dof;
sold[j] += phi[i]*dofold;
for (int k = 0; k < D; ++k)
ds[j*D+k] += dphi[i][k]*dof;
}
for (int k=0; k<D; ++k){
dc[k] += problem->porosity[el[i]]*dphi[i][k];
}
}
double f[D],dfdu,dfddp;
particle_force_f(problem,f,&dfdu,&dfddp,s,ds,sold,c,dc,a,dt,iel,ip);
if (!local_vector)
continue;
double rho,nu;
fluid_problem_interpolate_rho_and_nu(problem,a, &rho,&nu);
double pspg = drag_in_stab?problem->taup[iel]/rho:0;
for (int iD = 0; iD < D; ++iD) {
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];
}
}
}
if (!local_matrix)
continue;
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;
}
}
}
}
}
free(s);
free(ds);
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 *mesh_velocity, const double c, const double *dc, const double rho, double mu, const double dt, int eid, const double *data, double *f00, double *f01)
{
......@@ -381,7 +381,7 @@ static void particle_force_f(FluidProblem *problem, double *f, double *dfdu, dou
double Due[D];
const double *up = &problem->particle_velocity[ip*D];
for (int j = 0; j < D; ++j) {
Due[j] = (problem->particle_velocity[ip*D+j]-uold[j]);///c;
Due[j] = (problem->particle_velocity[ip*D+j]-uold[j])/c;
}
double mass = problem->particle_mass[ip];
double vol = problem->particle_volume[ip];
......@@ -402,25 +402,21 @@ static void particle_force_f(FluidProblem *problem, double *f, double *dfdu, dou
}
gamma = gamma*problem->drag_coeff_factor;
double gammastar = gamma/(1+dt/mass*gamma);
double r = D==2? sqrt(problem->particle_volume[ip]/M_PI):pow(3*problem->particle_volume[ip]/(4*M_PI),1/3);
double h = problem->element_size[iel];
double penalty_coeff = 1/(5*dt)*(r*r/(h*h))*(1-c)*0;
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];
// double upstar = up[i]+dt/mass*(contact[i]+massreduced*problem->g[i]-vol*dp[i]);
// double drag = gammastar*(up[i]-u[i]);//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] = -penalty_coeff*(upstar-u[i]);
f[U+i] = -contact[i];
}
*dfdu = penalty_coeff;
// *dfdu += gammastar/c;
*dfddp = penalty_coeff*dt/mass*vol;//-vol;
*dfddp = 0;
*dfdu = 0;
}
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) {
double p = sol[P];
double taup = problem->taup[iel]/16.0;
double tauc = problem->tauc[iel]/16.0;
double taup = problem->taup[iel];///16.0;
double tauc = problem->tauc[iel];///16.0;
const int n_fields = fluid_problem_n_fields(problem);
double u[D], uold[D], du[D][D], duold[D][D], umesh[D], dp[D];
if(problem->stab_param != 0)
......@@ -449,7 +445,7 @@ static void fluid_problem_f(const FluidProblem *problem, const double *sol, doub
double drag = problem->volume_drag + problem->quadratic_drag*nold;
for (int i = 0; i < D; ++i) {
f0[U+i] =
+ dp[i] // remove pressure here
+ dp[i]
- g[i]*rhoreduced
- bf[i]
+ drag*u[i];
......@@ -474,11 +470,6 @@ static void fluid_problem_f(const FluidProblem *problem, const double *sol, doub
f1[(U+i)*D+j] = mu*(du[i][j] + du[j][i]);
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;
// divide mu by two and add the symmetric contribution
// f1[(U+j)*D+i] = mu/2.0*(du[i][j] + du[j][i]);
// f11[((U+j)*D+i)*n_fields*D+(U+i)*D+j] += mu/2.0;
// f11[((U+j)*D+i)*n_fields*D+(U+j)*D+i] += mu/2.0;
// SUPG
double supg = uold[j]*taup;
f1[(U+i)*D+j] += supg*f0[U+i];
......@@ -1211,7 +1202,7 @@ void fluid_problem_assemble_system(FluidProblem *problem, double *rhs, const dou
forced_row[problem->boundaries[iphys].nodes[i]*n_fields+bnd->row] = bnd->field;
}
}
// fluid_problem_node_force_volume(problem, solution_old, dt, all_local_vector, all_local_matrix);
fluid_problem_node_force_volume(problem, solution_old, dt, all_local_vector, all_local_matrix);
compute_weak_boundary_conditions(problem, solution_old, dt, all_local_vector, all_local_matrix);
fluid_problem_volume(problem, solution_old, dt, all_local_vector, all_local_matrix);
if (problem->n_fluids == 2){
......
......@@ -130,7 +130,7 @@ def RequestData():
if vals_s is not None :
output.PointData.append(np.repeat(vals_s,2*nf),"Reaction_s")
n = points.shape[0]
pattern = np.ndarray([nf,4], np.int)
pattern = np.ndarray([nf,4], int)
for i in range(nf) :
j = (i+1)%nf
pattern[i,:] = (i*2,i*2+1,j*2+1,j*2)
......
......@@ -107,18 +107,18 @@ def _load_msh(mesh_file_name, lib, dim) :
cbname = (c_char_p*len(bname))(*(i.encode() for i in bname))
# remove nodes and boundaries not connected to elements
keepnodes = np.full(x.shape[0], False, np.bool)
keepnodes = np.full(x.shape[0], False, bool)
keepnodes[el] = True
newidx = np.cumsum(keepnodes)-1
el = newidx[el]
x = x[keepnodes,:]
periodic_nodes = newidx[periodic_nodes]
keepbnd = np.full(bnd.shape[0], True, np.bool)
keepbnd = np.full(bnd.shape[0], True, bool)
for i in range(dim) :
keepbnd = np.logical_and(keepbnd, keepnodes[bnd[:,i]])
bnd = newidx[bnd][keepbnd,:]
btag = btag[keepbnd]
is_parent = np.full([x.shape[0]],True,np.bool)
is_parent = np.full([x.shape[0]],True,bool)
is_parent[periodic_nodes[:,0]] = False
periodic_parent = np.cumsum(is_parent)-1
periodic_parent[periodic_nodes[:,0]] = periodic_parent[periodic_nodes[:,1]]
......
......@@ -75,14 +75,14 @@ rho = 1000 # fluid density
nu = 1e-6 # kinematic viscosity
# Numerical parameters
outf = 10 # number of iterations between output files
dt = 5e-4 # time step
outf = 50 # number of iterations between output files
dt = 1e-4 # time step
tEnd = 100 # final time
# Geometrical parameters
ly = 1e-1 # particles area height
lx = 4e-1 # particles area widht
H = 0.6 # domain height
H = 0.2 # domain height
#
# PARTICLE PROBLEM
......@@ -102,7 +102,7 @@ ii = 0
#
# FLUID PROBLEM
#
fluid = fluid.FluidProblem(2,g,[nu*rho],[rho])
fluid = fluid.FluidProblem(2,g,[nu*rho],[rho], drag_in_stab=0)
# Set the mesh geometry for the fluid computation
fluid.load_msh("mesh.msh")
fluid.set_wall_boundary("Bottom")
......@@ -117,6 +117,7 @@ tic = time.time()
def compute_drag():
pv, pdp, pc = fluid.solution_at_particles()
contacts = p.contact_forces()
dv = (p.velocity()-pv)/pc
mu = nu/rho
d = 2*(p.volume()/np.pi)**0.5
......@@ -140,46 +141,22 @@ def compute_drag():
# -[v+g*dt-Vdp*dt/m]*gamma/(1+gamma*dt/m) - Vdp
# -[v+g*dt-dp*dt*V/m]*1/(1/gamma+dt/m) - Vdp
gammastar = 1/(1/gamma+dt/p.mass())
dvstar = dv+dt*(g-pdp/rhop)
dvstar = dv+dt*(contacts+g-pdp/rhop)
print(np.mean(np.linalg.norm(gammastar*dvstar, axis=1)), np.mean(np.linalg.norm(gamma*dv,axis=1)))
return -gammastar*dvstar-p.volume()*pdp
#
# COMPUTATION LOOP
#
while t < tEnd :
# fluid.move_particles(p.position(), p.velocity(), p.contact_forces())
# fluid.implicit_euler(dt)
# fext = g*p.mass() + compute_drag()#fluid.compute_node_force(dt)
# p.iterate(dt, fext)
# time_integration.iterate(fluid, p, dt, min_nsub=5, external_particles_forces=g*p.mass())
fluid.move_particles(p.position(), p.velocity(), p.contact_forces())
fluid.implicit_euler(dt)
fext = g*p.mass() + compute_drag()#fluid.compute_node_force(dt)
fext = g*p.mass() + compute_drag() #+ fluid.compute_node_force(dt) #g*p.mass() + compute_drag()
time_integration._advance_particles(p, fext, dt, 5, 1e-8)
# predicteur correcteur
# alpha = 0.5
# sf0 = np.copy(fluid.solution())
# p.save_state()
# predicteur
# fluid.implicit_euler(dt)
# sf1 = np.copy(fluid.solution())
# f0 = g*p.mass() + compute_drag()
# time_integration._advance_particles(p,f0,dt,min_nsub=5,contact_tol=1e-8)
# fluid.particle_velocity()[:] = p.velocity()
# correcteur
# fluid.solution()[:,:] = sf0
# fluid.implicit_euler(dt)
# fluid.solution()[:,:] = (alpha*fluid.solution()+(1-alpha)*sf1)
# f1 = g*p.mass() + compute_drag()
# p.restore_state()
# time_integration._advance_particles(p,f1,dt,min_nsub=5,contact_tol=1e-8)
t += dt
# Output files writting
if ii %outf == 0 :
......
L = 0.2;
H = 0.3;
H = 0.1;
y = 0;
lc = 0.01;
......
......@@ -187,7 +187,7 @@ while t < tEnd :
fluid.move_particles(p.position(), p.velocity(), p.contact_forces())
fluid.implicit_euler(dt)
fext = g*p.mass() + compute_drag()#fluid.compute_node_force(dt)
fext = g*p.mass() + fluid.compute_node_force(dt) #g*p.mass() + compute_drag()
time_integration._advance_particles(p, fext, dt, 5, 1e-8)
t += dt
......
......@@ -66,7 +66,7 @@ def genInitialPosition(filename, r, rout, rhop, compacity) :
p.write_vtk(filename,0,0)
# Define output directory
outputdir = "output_adapt_tau_by_16"
outputdir = "output_2"
if not os.path.isdir(outputdir) :
os.makedirs(outputdir)
......@@ -82,8 +82,8 @@ mu = nu*rho # dynamic viscosity
print("RHOP = %g, r = %g, RHO = %g, MU = %g" % (rhop,r,rho, mu))
# Numerical parameters
outf = 10 # number of iterations between output files
dt = 2.5e-3 # time step
outf = 2 # number of iterations between output files
dt = 10e-3 # time step
tEnd = 100 # final time
#
......@@ -171,17 +171,17 @@ def compute_drag():
while t < tEnd :
if ((ii+1)%200==0):
if ((ii+1)%50==0):
number_p = fluid.n_particles
position_p = fluid.particle_position()
volume_p = fluid.particle_volume()
if (ii%200==0 and ii != 0):
if (ii%50==0 and ii != 0):
write_points(p)
fluid.adapt_mesh(2e-2,1e-3,10000,old_n_particles=number_p,old_particle_position=position_p,old_particle_volume=volume_p)
fluid.move_particles(p.position(), p.velocity(), p.contact_forces())
fluid.implicit_euler(dt)
fext = g*p.mass() + compute_drag()#fluid.compute_node_force(dt)
fext = g*p.mass() + fluid.compute_node_force(dt) #g*p.mass() + compute_drag()
time_integration._advance_particles(p, fext, dt, 5, 1e-8)
t += dt
......
......@@ -110,7 +110,7 @@ class Darcy(unittest.TestCase) :
U = 0.00001 # averaged fluid velocity
V = U/(1-compacity) # fluid velocity
fluid.solution()[:,0] = U
fluid.set_open_boundary("Left",velocity=[V,0])
fluid.set_open_boundary("Left",velocity=[U,0])
fluid.set_open_boundary("Right",pressure=0)
fluid.set_wall_boundary("Bottom",velocity=[0,0])
fluid.set_wall_boundary("Top",velocity=[0,0])
......@@ -118,12 +118,36 @@ class Darcy(unittest.TestCase) :
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity(), p.contact_forces())
fluid.write_vtk(outputdir,0,0)
def compute_drag():
pv, pdp, pc = fluid.solution_at_particles()
contacts = p.contact_forces()
dv = (p.velocity()-pv)/pc
mu = nu/rho
d = 2*(p.volume()/np.pi)**0.5
normu = np.hypot(dv[0],dv[1])
#Reynolds/|u_p-u_f|
Re_O_u = d*pc/mu*rho
Cd_u = 0.63*np.sqrt(normu)+4.8/np.sqrt(Re_O_u)
Cd_u = Cd_u*Cd_u
f = pc**(-(1.8-0.65*np.exp(-(1.5-np.log(Re_O_u*normu))*(1.5-np.log(Re_O_u*normu))/2.)))
gamma = f*Cd_u*rho/2*d
gammastar = 1/(1/gamma+dt/p.mass())
dvstar = dv+dt*(contacts+g-pdp/rhop)
print(np.mean(np.linalg.norm(gammastar*dvstar, axis=1)), np.mean(np.linalg.norm(gamma*dv,axis=1)))
return -gammastar*dvstar-p.volume()*pdp
t = 0
ii = 0
tic = time.time()
#Computation loop
fext = 0
while t < tEnd :
#Fluid solver
fluid.move_particles(p.position(), p.velocity(), -fext)
fluid.implicit_euler(dt)
fext = g*p.mass() + compute_drag() #+ fluid.compute_node_force(dt) #g*p.mass() + compute_drag()
# time_integration._advance_particles(p, fext, dt, 5, 1e-8)
time_integration.iterate(fluid,p,dt,fixed_grains=True)
t += dt
#Output files writting
......@@ -134,10 +158,10 @@ class Darcy(unittest.TestCase) :
ii += 1
print("%i : %.2g/%.2g (cpu %.6g)" % (ii, t, tEnd, time.time() - tic))
gamma = (0.63*U/(1-compacity)+4.8*(mu/((1-compacity)*rho*2*r))**0.5)**2*(1-compacity)**-1.8*rho*r
gamma = (0.63*U + 4.8*(mu/((1-compacity)*rho*2*r))**0.5)**2*(1-compacity)**-1.8*rho*r
vol = (np.pi*r**2)
N = compacity/vol
Dp_drag[j] = gamma*N*U/(1-compacity)**2
Dp_drag[j] = gamma*N*U**2
pressure = fluid.solution()[:,2]
pM[j] = np.mean(fluid.solution()[fluid.coordinates()[:,0]<-.99,2])
......
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