Skip to content
GitLab
Commit 050ca3b0 authored by Jonathan Lambrechts's avatar Jonathan Lambrechts
Browse files

dcdt in n+1 - n

parent a35928ca
Hide whitespace changes
Inline Side-by-side
scontact/scontact.c
View file @ 050ca3b0
......@@ -686,7 +686,6 @@ static void _particleProblemGenContacts(ParticleProblem *p, const double alert)
}
}
}
printf("time old %i\n", clock()-tic);
//exit(0);
// Disks
for (size_t i = 0; i < vectorSize(p->disks); ++i) {
......
src/fluid_problem.c
View file @ 050ca3b0
......@@ -161,9 +161,8 @@ static void compute_node_force_implicit(FluidProblem *problem, double dt, double
gradp[j] += dphi[iphi][j]*si[el[iphi]*n_fields+DIMENSION];
}
}
double du[DIMENSION], due[DIMENSION];
double due[DIMENSION];
for (int j = 0; j < DIMENSION; ++j) {
du[j] = problem->particle_velocity[i*DIMENSION+j]-vf[j]/c;
due[j] = problem->particle_velocity[i*DIMENSION+j]-vfe[j]/c;
}
double vol = problem->particle_volume[i];
......@@ -173,18 +172,27 @@ static void compute_node_force_implicit(FluidProblem *problem, double dt, double
double gamma = particle_drag_coeff(due,problem->mu,problem->rho,vol,c);
double g = problem->g*drho/rhop;
double fcoeff = mass/(mass+dt*gamma);
double drag[DIMENSION];
for (int j = 0; j < DIMENSION; ++j) {
problem->particle_new_velocity[i*DIMENSION+j] =
fcoeff * (
problem->particle_velocity[i*DIMENSION+j]
+gamma*vf[j]*dt/c/mass-vol*gradp[j]*dt/mass
+(j == 1 ? g*dt : 0.)
);
drag[j] = gamma*(problem->particle_new_velocity[i*DIMENSION+j]-vf[j]/c);
}
for (int d = 0; d < DIMENSION; ++d)
problem->particle_force[i*DIMENSION+d] = fcoeff*(-gamma*du[d]-vol*gradp[d]);
problem->particle_force[i*DIMENSION+1] += fcoeff*g*mass;
problem->particle_force[i*DIMENSION+d] = -drag[d] -vol*gradp[d];
problem->particle_force[i*DIMENSION+1] += g*mass;
double *local_vector = all_local_vector+iel*n_fields*N_SF;
double *local_matrix = all_local_matrix+iel*n_fields*n_fields*N_SF*N_SF;
for (int iphi = 0; iphi < N_SF; ++iphi) {
for (int d = 0; d < DIMENSION; ++d)
local_vector[iphi+N_SF*d] -= fcoeff*gamma*(du[d]-dt/mass*vol*gradp[d])*phi[iphi];
local_vector[iphi+N_SF*1] -= fcoeff*gamma*dt*g*phi[iphi];
local_vector[iphi+N_SF*d] -= drag[d]*phi[iphi];
}
for (int iphi = 0; iphi < N_SF; ++iphi) {
for (int jphi = 0; jphi < N_SF; ++jphi) {
......@@ -199,6 +207,37 @@ static void compute_node_force_implicit(FluidProblem *problem, double dt, double
}
}
}
//compute new porosity
double *new_position = malloc(sizeof(double)*problem->n_particles*DIMENSION);
double *new_uvw = malloc(sizeof(double)*problem->n_particles*DIMENSION);
int *new_eid = malloc(sizeof(int)*problem->n_particles);
for (int i = 0; i < problem->n_particles; ++i) {
new_eid[i] = problem->particle_element_id[i];
for (int d = 0; d < DIMENSION; ++d) {
new_position[i*DIMENSION+d] = problem->particle_position[i*DIMENSION+d]+dt*problem->particle_new_velocity[i*DIMENSION+d];
}
}
mesh_tree_particle_to_mesh(problem->mesh_tree,problem->n_particles, new_position, new_eid, new_uvw);
{
double *volume = problem->particle_volume;
for (int i = 0; i < mesh->n_nodes; ++i){
problem->new_porosity[i] = 0;
}
for (int i = 0; i < problem->n_particles; ++i) {
if(new_eid[i] == -1) continue;
double sf[N_SF];
shape_functions(&new_uvw[i*DIMENSION],sf);
const uint32_t *el = &mesh->elements[new_eid[i]*N_N];
for (int j=0; j<N_N;++j)
problem->new_porosity[el[j]] += sf[j]*volume[i];
}
for (int i = 0; i < mesh->n_nodes; ++i){
problem->new_porosity[i] = (1-problem->new_porosity[i]/problem->node_volume[i]);
}
}
free(new_eid);
free(new_uvw);
free(new_position);
}
static void fluid_problem_assemble_system(FluidProblem *problem, double *rhs, const double *solution_old, double dt)
......@@ -206,8 +245,6 @@ static void fluid_problem_assemble_system(FluidProblem *problem, double *rhs, co
HXTLinearSystem *lsys = problem->linear_system;
const Mesh *mesh = problem->mesh;
const double *solution = problem->solution;
const double *porosity = problem->porosity;
const double *old_porosity = problem->old_porosity;
const double mu = problem->mu;
const double rho = problem->rho;
const double epsilon = problem->epsilon;
......@@ -233,7 +270,7 @@ static void fluid_problem_assemble_system(FluidProblem *problem, double *rhs, co
double dc[DIMENSION]={0}, du[DIMENSION][DIMENSION]={{0}}, dp[DIMENSION]={0};
for (int i = 0; i< N_SF; ++i) {
for (int j = 0; j < DIMENSION; ++j) {
dc[j] += dphi[i][j]*porosity[el[i]];
dc[j] += dphi[i][j]*problem->porosity[el[i]];
dp[j] += dphi[i][j]*solution[el[i]*n_fields+DIMENSION];
for (int k = 0; k < DIMENSION; ++k)
du[k][j] += dphi[i][j]*solution[el[i]*n_fields+k];
......@@ -249,8 +286,8 @@ static void fluid_problem_assemble_system(FluidProblem *problem, double *rhs, co
dudt[j] += (solution[el[i]*n_fields+j]-solution_old[el[i]*n_fields+j])/dt*phi[i];
}
p += solution[el[i]*n_fields+DIMENSION]*phi[i];
c += porosity[el[i]]*phi[i];
dcdt += (porosity[el[i]]-old_porosity[el[i]])/dt*phi[i];
c += problem->porosity[el[i]]*phi[i];
dcdt += (problem->new_porosity[el[i]]-problem->porosity[el[i]])/dt*phi[i];
}
const double jw = QW[iqp]*detj;
for (int iphi = 0; iphi < N_SF; ++iphi) {
......@@ -395,10 +432,10 @@ int fluid_problem_implicit_euler(FluidProblem *problem, double dt)
norm0 = norm;
else if(norm/norm0 < newton_rtol)
break;
clock_get(&tic);
clock_get(&tic);
hxtLinearSystemSolve(problem->linear_system, rhs, dx);
clock_get(&toc);
totaltime += clock_diff(&tic, &toc);
clock_get(&toc);
totaltime += clock_diff(&tic, &toc);
for (int j=0; j<mesh->n_nodes*n_fields; ++j) {
problem->solution[j] -= dx[j];
}
......@@ -482,18 +519,17 @@ FluidProblem *fluid_problem_new(const char *mesh_file_name, double g, double mu,
problem->strong_boundaries[i].field = strong_boundaries[i].field;
}
problem->porosity = malloc(mesh->n_nodes*sizeof(double));
problem->old_porosity = malloc(mesh->n_nodes*sizeof(double));
problem->new_porosity = malloc(mesh->n_nodes*sizeof(double));
problem->solution = malloc(mesh->n_nodes*n_fields*sizeof(double));
for (int i = 0; i < problem->mesh->n_nodes; ++i)
{
problem->porosity[i] = 1.;
problem->old_porosity[i] = 1.;
}
// begin to remove
for (int i = 0; i < problem->mesh->n_nodes*N_N; ++i){
problem->solution[i] = 0.;
}
problem->node_volume = malloc(0);
problem->node_volume = NULL;
fluid_problem_compute_node_volume(problem);
// end to remove
problem->linear_system = NULL;
......@@ -503,13 +539,14 @@ FluidProblem *fluid_problem_new(const char *mesh_file_name, double g, double mu,
hxtLinearSystemCreateLU(&problem->linear_system,mesh->n_elements,N_N,n_fields,mesh->elements);
#endif
problem->n_particles = 0;
problem->particle_uvw = malloc(0);
problem->particle_element_id = malloc(0);
problem->particle_position = malloc(0);
problem->particle_mass = malloc(0);
problem->particle_force = malloc(0);
problem->particle_volume = malloc(0);
problem->particle_velocity = malloc(0);
problem->particle_uvw = NULL;
problem->particle_element_id = NULL;
problem->particle_position = NULL;
problem->particle_mass = NULL;
problem->particle_force = NULL;
problem->particle_volume = NULL;
problem->particle_velocity = NULL;
problem->particle_new_velocity = NULL;
return problem;
}
......@@ -517,7 +554,7 @@ void fluid_problem_free(FluidProblem *problem)
{
free(problem->solution);
free(problem->porosity);
free(problem->old_porosity);
free(problem->new_porosity);
hxtLinearSystemDelete(&problem->linear_system);
free(problem->particle_uvw);
free(problem->particle_element_id);
......@@ -713,7 +750,7 @@ void fluid_problem_adapt_mesh(FluidProblem *problem, double lcmax, double lcmin,
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));
//printf("Time mesh_load : %.6e \n",clock_diff(&tic,&toc));
double *new_solution = malloc(sizeof(double)*new_mesh->n_nodes*N_N);
double *new_xi = malloc(sizeof(double)*new_mesh->n_nodes*DIMENSION);
clock_get(&tic);
......@@ -721,18 +758,18 @@ void fluid_problem_adapt_mesh(FluidProblem *problem, double lcmax, double lcmin,
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));
//printf("Time new_eid : %.6e \n",clock_diff(&tic,&toc));
clock_get(&tic);
double *newx = malloc(sizeof(double)*DIMENSION*new_mesh->n_nodes);
for (int i = 0;i < new_mesh->n_nodes;++i)
for (int k = 0; k < DIMENSION; ++k)
newx[i*DIMENSION+k] = new_mesh->x[i*3+k];
clock_get(&toc);
printf("Time newx : %.6e \n",clock_diff(&tic,&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));
//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) {
......@@ -752,14 +789,13 @@ void fluid_problem_adapt_mesh(FluidProblem *problem, double lcmax, double lcmin,
free(newx);
free(problem->solution);
problem->solution = new_solution;
free(problem->new_porosity);
problem->new_porosity = malloc(sizeof(double)*new_mesh->n_nodes);
free(problem->porosity);
problem->porosity = malloc(sizeof(double)*new_mesh->n_nodes);
free(problem->old_porosity);
problem->old_porosity = malloc(sizeof(double)*new_mesh->n_nodes);
for (int i = 0; i < new_mesh->n_nodes; ++i)
{
problem->porosity[i] = 1.;
problem->old_porosity[i] = 1.;
}
mesh_free(problem->mesh);
problem->mesh = new_mesh;
......@@ -782,8 +818,8 @@ void fluid_problem_after_import(FluidProblem *problem)
{
mesh_tree_free(problem->mesh_tree);
problem->mesh_tree = mesh_tree_create(problem->mesh);
free(problem->old_porosity);
problem->old_porosity = malloc(sizeof(double)*problem->mesh->n_nodes);
free(problem->new_porosity);
problem->new_porosity = malloc(sizeof(double)*problem->mesh->n_nodes);
for (int i = 0; i < problem->n_particles; ++i)
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);
......@@ -812,6 +848,7 @@ void fluid_problem_set_particles(FluidProblem *problem, int n, double *mass, dou
free(problem->particle_force);
free(problem->particle_volume);
free(problem->particle_velocity);
free(problem->particle_new_velocity);
problem->particle_uvw = malloc(sizeof(double)*DIMENSION*n);
problem->particle_element_id = malloc(sizeof(int)*n);
for (int i = 0; i < n; ++i) {
......@@ -822,6 +859,7 @@ void fluid_problem_set_particles(FluidProblem *problem, int n, double *mass, dou
problem->particle_force = malloc(sizeof(double)*n*DIMENSION);
problem->particle_volume = malloc(sizeof(double)*n);
problem->particle_velocity = malloc(sizeof(double)*n*DIMENSION);
problem->particle_new_velocity = malloc(sizeof(double)*n*DIMENSION);
}
if (elid) {
for (int i = 0; i < n; ++i) {
......@@ -831,7 +869,7 @@ void fluid_problem_set_particles(FluidProblem *problem, int n, double *mass, dou
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));
//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];
......@@ -842,14 +880,11 @@ void fluid_problem_set_particles(FluidProblem *problem, int n, double *mass, dou
}
}
clock_get(&toc);
printf("Time initialize particle_fluid : %.6e \n",clock_diff(&tic,&toc));
for (int i = 0; i < problem->mesh->n_nodes; ++i){
problem->old_porosity[i] = problem->porosity[i];
}
//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));
//printf("Time compute_porosity : %.6e \n",clock_diff(&tic,&toc));
}
int *fluid_problem_particle_element_id(FluidProblem *problem)
......
src/fluid_problem.h
View file @ 050ca3b0
......@@ -26,7 +26,7 @@ typedef struct {
Mesh *mesh;
MeshTree *mesh_tree;
double *porosity;
double *old_porosity;
double *new_porosity;
double *solution;
double *solution_explicit;
double *node_volume;
......@@ -39,6 +39,8 @@ typedef struct {
double *particle_position;
double *particle_volume;
double *particle_velocity;
double *particle_new_velocity;
double *particle_drag;
double *particle_mass;
double *particle_force;
int *particle_element_id;
......
testcases/drop-2d/drop.py
View file @ 050ca3b0
......@@ -76,14 +76,12 @@ forces = g*p.mass()
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())
while t < tEnd :
if (ii%5==0 and ii != 0):
#fluid.adapt_mesh(0.01,200,1e-4)
fluid.adapt_mesh(1e-3,1e-4,10000, epsilon, 1)
fluid.adapt_mesh(1e-3,1e-4,10000)
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())
else :
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())
fluid.implicit_euler(dt)
forces = fluid.compute_node_force(dt)
print(forces)
vn = p.velocity() + forces * dt / p.mass()
vmax = np.max(np.hypot(vn[:, 0], vn[:, 1]))
nsub = max(1, int(np.ceil((vmax * dt * 4)/min(10*p.r()))))
......
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