Commit 8556a2ca authored by Jonathan Lambrechts's avatar Jonathan Lambrechts
Browse files

porosity in solution

parent a3c29a8c
......@@ -80,14 +80,14 @@ class fluid_problem :
def implicit_euler(self, dt) :
return lib.fluid_problem_implicit_euler(self._fp, c_double(dt))
def set_particles(self, mass, volume, position, velocity):
def set_particles(self, mass, volume, position, velocity, init=False):
def np2c(a) :
tmp = numpy.require(a,"float","C")
r = c_void_p(tmp.ctypes.data)
r.tmp = tmp
return r
self.n_particles = mass.shape[0]
lib.fluid_problem_set_particles(self._fp,c_int(mass.shape[0]),np2c(mass),np2c(volume),np2c(position),np2c(velocity),None)
lib.fluid_problem_set_particles(self._fp,c_int(mass.shape[0]),np2c(mass),np2c(volume),np2c(position),np2c(velocity),None,c_int(init))
def particle_element_id(self) :
......
......@@ -30,7 +30,7 @@
#include "mesh_find.h"
#include "mbxml.h"
#define n_fields (DIMENSION+1)
#define n_fields (DIMENSION+2)
#define newton_max_it 10
#define newton_atol 5e-7
......@@ -265,18 +265,22 @@ static void fluid_problem_assemble_system(FluidProblem *problem, double *rhs, co
du[k][j] += dphi[i][j]*solution[el[i]*n_fields+k];
}
}
int P = DIMENSION;
int Q = DIMENSION+1;
for (int iqp = 0; iqp< N_QP; ++iqp) {
double phi[N_SF];
shape_functions(QP[iqp],phi);
double u[DIMENSION]={0}, dudt[DIMENSION]={0}, p=0, c=0, dcdt=0;
double u[DIMENSION]={0}, dudt[DIMENSION]={0}, p=0, c=0, dcdt=0,q=0;
for (int i = 0; i < N_SF; ++i) {
for (int j = 0; j < DIMENSION; ++j){
u[j] += solution[el[i]*n_fields+j]*phi[i];
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];
p += solution[el[i]*n_fields+P]*phi[i];
q += solution[el[i]*n_fields+Q]*phi[i];
c += porosity[el[i]]*phi[i];
dcdt += (porosity[el[i]]-old_porosity[el[i]])/dt*phi[i];
dcdt += (solution[el[i]*n_fields+Q]-solution_old[el[i]*n_fields+Q])/dt*phi[i];
//dcdt += (problem->porosity[el[i]]-problem->old_porosity[el[i]])/dt*phi[i];
}
const double jw = QW[iqp]*detj;
for (int iphi = 0; iphi < N_SF; ++iphi) {
......@@ -286,6 +290,11 @@ static void fluid_problem_assemble_system(FluidProblem *problem, double *rhs, co
for (int i = 0; i < DIMENSION; ++i)
for (int j = 0; j < DIMENSION; ++j)
tau[i][j] = du[i][j]-u[i]*dc[j]/c;
double dphidp = 0, divu = 0;
for (int k = 0; k < DIMENSION; ++k){
dphidp += dphii[k]*dp[k];
divu += du[k][k];
}
for (int j = 0; j < DIMENSION; ++j) {
double utau = 0;
double dphisigma=0;
......@@ -297,16 +306,12 @@ static void fluid_problem_assemble_system(FluidProblem *problem, double *rhs, co
local_vector[iphi+N_SF*j] += jw*(
dphisigma
//+rho*phii*dudt[j]-rho*utau
+rho*phii*(dudt[j]-u[j]/c*dcdt+utau/c)
+rho*phii*(dudt[j]+u[j]/c*divu+utau/c)
-p*dphii[j]
);
}
double dphidp = 0, divu = 0;
for (int k = 0; k < DIMENSION; ++k){
dphidp += dphii[k]*dp[k];
divu += du[k][k];
}
local_vector[iphi+N_SF*DIMENSION] += jw*(epsilon*dphidp+phii*(divu+dcdt));
local_vector[iphi+N_SF*P] += jw*(epsilon*dphidp+phii*(divu+dcdt));
local_vector[iphi+N_SF*Q] += jw*phii*(c-q);
for (int jphi = 0; jphi < N_SF; ++jphi) {
double phij = phi[jphi];
const double *dphij = dphi[jphi];
......@@ -333,12 +338,17 @@ static void fluid_problem_assemble_system(FluidProblem *problem, double *rhs, co
for (int k = 0; k < DIMENSION; ++k) {
int V = k;
//LOCAL_MATRIX(U,V) += -jw*rho*u[j]*dphii[k]*phij/c+jw*2*mu*dphii[k]*dtau[j]*0.5;
LOCAL_MATRIX(U,V) += jw*rho*phii*phij*tau[j][k]/c+jw*2*mu*dphii[k]*dtau[j]*0.5;
LOCAL_MATRIX(U,V) +=
jw*rho*phii*phij*tau[j][k]/c
+jw*2*mu*dphii[k]*dtau[j]*0.5
+jw*rho*phii*u[j]*dphij[k]/c;
}
//LOCAL_MATRIX(U,U) += jw*mu*2*0.5*dphiidtau + jw*rho*(phii*phij/dt-phij*dphii[j]*u[j]/c);
LOCAL_MATRIX(U,U) += jw*mu*2*0.5*dphiidtau + jw*rho*phii*(phij/dt-phij/c*dcdt+udtau/c);
LOCAL_MATRIX(U,U) += jw*mu*2*0.5*dphiidtau + jw*rho*phii*(phij/dt+phij/c*divu+udtau/c);
}
LOCAL_MATRIX(P,P) += jw*epsilon*dphiidphij;
LOCAL_MATRIX(P,P) += jw*epsilon*dphiidphij;
LOCAL_MATRIX(Q,Q) += -jw*phii*phij;
LOCAL_MATRIX(P,Q) += jw*phii*phij/dt;
}
}
}
......@@ -479,7 +489,7 @@ static void fluid_problem_compute_porosity(FluidProblem *problem)
problem->porosity[i] = 1.;
}
else{
problem->porosity[i] = (1-problem->porosity[i]/problem->node_volume[i] /4.);
problem->porosity[i] = (1-problem->porosity[i]/problem->node_volume[i]);
}
}
}
......@@ -521,7 +531,7 @@ FluidProblem *fluid_problem_new(const char *mesh_file_name, double g, double mu,
problem->old_porosity[i] = 1.;
}
// begin to remove
for (int i = 0; i < problem->mesh->n_nodes*N_N; ++i){
for (int i = 0; i < problem->mesh->n_nodes*n_fields; ++i){
problem->solution[i] = 0.;
}
problem->node_volume = malloc(0);
......@@ -745,7 +755,7 @@ void fluid_problem_adapt_mesh(FluidProblem *problem, double lcmax, double lcmin,
Mesh *new_mesh = mesh_load_msh("adapt/mesh.msh");
clock_get(&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_solution = malloc(sizeof(double)*new_mesh->n_nodes*n_fields);
double *new_xi = malloc(sizeof(double)*new_mesh->n_nodes*DIMENSION);
clock_get(&tic);
int *new_eid = malloc(sizeof(int)*new_mesh->n_nodes);
......@@ -829,7 +839,7 @@ void fluid_problem_after_import(FluidProblem *problem)
}
void fluid_problem_set_particles(FluidProblem *problem, int n, double *mass, double *volume, double *position, double *velocity, long int *elid)
void fluid_problem_set_particles(FluidProblem *problem, int n, double *mass, double *volume, double *position, double *velocity, long int *elid, int init)
{
struct timespec tic, toc;
......@@ -877,6 +887,11 @@ void fluid_problem_set_particles(FluidProblem *problem, int n, double *mass, dou
for (int i = 0; i < problem->mesh->n_nodes; ++i){
problem->old_porosity[i] = problem->porosity[i];
}
if (init) {
int Q = DIMENSION+1;
for (int i = 0; i < problem->mesh->n_nodes; ++i)
problem->solution[i*n_fields+Q] = problem->porosity[i];
}
clock_get(&tic);
fluid_problem_compute_porosity(problem);
clock_get(&toc);
......
......@@ -73,7 +73,7 @@ int fluid_problem_import_vtk(FluidProblem *problem, const char *filename);
// complete force on the particle (including gravity)
void fluid_problem_compute_node_particle_force(FluidProblem *problem, double dt, double *particle_force);
int fluid_problem_implicit_euler(FluidProblem *problem, double dt);
void fluid_problem_set_particles(FluidProblem *problem, int n, double *mass, double *volume, double *position, double *velocity, long int *elid);
void fluid_problem_set_particles(FluidProblem *problem, int n, double *mass, double *volume, double *position, double *velocity, long int *elid, int init);
void fluid_problem_free(FluidProblem *problem);
FluidProblem *fluid_problem_new(const char *mesh_file_name, double g, double mu, double rho, double epsilon, int n_strong_boundaries, StrongBoundary *strong_boundaries);
void fluid_problem_adapt_mesh(FluidProblem *problem, double lcmax, double lcmin, double n_el);
......
......@@ -228,13 +228,13 @@ int fluid_problem_export_vtk(const FluidProblem *problem, const char *output_dir
fprintf(f,"\n</DataArray>\n");
fprintf(f,"<DataArray Name=\"Pressure\" NumberOfComponents = \"1\" type=\"Float64\" format=\"ascii\">\n");
for (int i = 0; i < mesh->n_nodes; ++i) {
fprintf(f, "%g ", problem->solution[i*(DIMENSION+1)+DIMENSION]);
fprintf(f, "%g ", problem->solution[i*(DIMENSION+2)+DIMENSION]);
}
fprintf(f,"\n</DataArray>\n");
fprintf(f,"<DataArray Name=\"Velocity\" NumberOfComponents = \"3\" type=\"Float64\" format=\"ascii\">\n");
for (int i = 0; i < mesh->n_nodes; ++i) {
for (int j = 0; j < DIMENSION; ++j)
fprintf(f, "%g ", problem->solution[i*(DIMENSION+1)+j]);
fprintf(f, "%g ", problem->solution[i*(DIMENSION+2)+j]);
if(DIMENSION==2)
fprintf(f, "0 ");
}
......
......@@ -66,7 +66,7 @@ t = 0
ii = 0
#Enable computation of the initial velocity
initialize = 1
initialize = 0
#physical parameters
g = -9.81 # gravity
......@@ -100,7 +100,7 @@ outf = 1 # number of iterations between o
#Format: strong_boundaries = [(Boundary tag, Fluid field, Value)
strong_boundaries = [("Top",2,0.),("Top",1,0.),("Bottom",1,0.),("Lateral",0,0.)]
fluid = fluid.fluid_problem("mesh.msh",g,nu*rho,rho,epsilon,strong_boundaries)
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity(),True)
ii = 0
jj = 0
......@@ -141,7 +141,7 @@ fluid.export_vtk(outputdir,0,0)
ii = 0
tic = time.clock()
#Computation loop
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity(),True)
while t < tEnd :
#Adaptation of the mesh. Args are minimal mesh radius, maximal mesh radius and number of elements
if (ii%5==0 and ii != 0):
......
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