Commit 02115f3b authored by Matthieu Constant's avatar Matthieu Constant
Browse files

debut deux fluide (rho+mu+nfields+debut assemble system)

parent 244a8681
This diff is collapsed.
......@@ -30,7 +30,6 @@
#include "mesh_find.h"
#include "mbxml.h"
#define n_fields (DIMENSION+2)
#define newton_max_it 10
#define newton_atol 5e-7
......@@ -157,6 +156,7 @@ void fluid_problem_compute_node_particle_force(FluidProblem *problem, double dt,
static void compute_node_force_implicit(FluidProblem *problem, double dt, double *all_local_vector, double *all_local_matrix) {
double *porosity = problem->porosity;
Mesh *mesh = problem->mesh;
int n_fields = problem->nFields;
for (int i = 0; i < problem->n_particles; ++i) {
int iel = problem->particle_element_id[i];
double mass = problem->particle_mass[i];
......@@ -232,12 +232,14 @@ 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;
int n_fields = problem->nFields;
const double *compacity = problem->compacity;
const double *old_compacity = problem->old_compacity;
const double mu = problem->mu;
const double rho = problem->rho;
const double epsilon = problem->epsilon;
size_t local_size = N_SF*n_fields;
int n_fluids = problem->nFluids;
size_t local_size = N_SF*n_fields*problem->nFluids;
size_t N_E = mesh->n_elements;
double *all_local_vector = malloc(N_E*local_size*sizeof(double));
double *all_local_matrix = malloc(N_E*local_size*local_size*sizeof(double));
......@@ -249,119 +251,123 @@ static void fluid_problem_assemble_system(FluidProblem *problem, double *rhs, co
for (int iel=0; iel < mesh->n_elements; ++iel) {
const unsigned int *el = &mesh->elements[iel*N_N];
double dxdxi[DIMENSION][DIMENSION], dxidx[DIMENSION][DIMENSION], dphi[N_N][DIMENSION];
for (int i = 0; i < DIMENSION; ++i)
for (int j = 0; j < DIMENSION; ++j)
dxdxi[i][j] = mesh->x[el[j+1]*3+i] - mesh->x[el[0]*3+i];
const double detj = invDD(dxdxi, dxidx);
grad_shape_functions(dxidx, dphi);
double *local_vector = &all_local_vector[N_SF*n_fields*iel];
double *local_matrix = &all_local_matrix[N_SF*N_SF*n_fields*n_fields*iel];
double dc[DIMENSION]={0}, du[DIMENSION][DIMENSION]={{0}}, dp[DIMENSION]={0}, dcdq[DIMENSION]={0};
for (int i = 0; i< N_SF; ++i) {
for (int j = 0; j < DIMENSION; ++j) {
dc[j] += dphi[i][j]*solution[el[i]*n_fields+DIMENSION+1];
dp[j] += dphi[i][j]*solution[el[i]*n_fields+DIMENSION];
dcdq[j] += dphi[i][j];
for (int k = 0; k < DIMENSION; ++k)
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,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];
for (int ff = 0; ff < n_fluids; ++ff)
{
int P = n_fields;
int Q = DIMENSION+ff*(DIMENSION+1);
for (int i = 0; i < DIMENSION; ++i)
for (int j = 0; j < DIMENSION; ++j)
dxdxi[i][j] = mesh->x[el[j+1]*3+i] - mesh->x[el[0]*3+i];
const double detj = invDD(dxdxi, dxidx);
grad_shape_functions(dxidx, dphi);
double *local_vector = &all_local_vector[N_SF*n_fields*n_fluids*iel+ff*N_SF*(DIMENSION+1)];
double *local_matrix = &all_local_matrix[N_SF*N_SF*n_fields*n_fields*n_fluids*n_fluids*iel+ff*N_SF*N_E*(DIMENSION+1)*(DIMENSION+1)];
double dc[DIMENSION]={0}, du[DIMENSION][DIMENSION]={{0}}, dp[DIMENSION]={0}, dcdq[DIMENSION]={0};
for (int i = 0; i< N_SF; ++i) {
for (int j = 0; j < DIMENSION; ++j) {
dc[j] += dphi[i][j]*solution[el[i]*n_fields+P];
dp[j] += dphi[i][j]*solution[el[i]*n_fields+Q];
dcdq[j] += dphi[i][j];
for (int k = 0; k < DIMENSION; ++k)
du[k][j] += dphi[i][j]*solution[el[i]*n_fields+k+ff*(DIMENSION+1)];
}
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 += (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) {
double phii = phi[iphi];
const double *dphii = dphi[iphi];
double tau[DIMENSION][DIMENSION];
for (int i = 0; i < DIMENSION; ++i)
for (int j = 0; j < DIMENSION; ++j)
tau[i][j] = du[i][j]-u[i]*dc[j]/q;
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;
for (int k = 0; k < DIMENSION; ++k) {
//utau += u[k]*dphii[k]*u[j]/c;
utau += u[k]*tau[j][k];
dphisigma += 2*mu*dphii[k]*0.5*(tau[k][j]+tau[j][k]);
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,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+ff*(DIMENSION+1)]*phi[i];
dudt[j] += (solution[el[i]*n_fields+j+ff*(DIMENSION+1]-solution_old[el[i]*n_fields+j+ff*(DIMENSION+1])/dt*phi[i];
}
local_vector[iphi+N_SF*j] += jw*(
dphisigma
//+rho*phii*dudt[j]-rho*utau
+rho*phii*(dudt[j]+u[j]/q*divu+utau/q)
-p*dphii[j]
);
p += solution[el[i]*n_fields+P]*phi[i];
q += solution[el[i]*n_fields+Q]*phi[i];
c += compacity[el[i]]*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];
}
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];
int P = DIMENSION;
//int N2 = N_SF*N_SF;
//int IDX = iphi*N2+jphi;
int N2 = n_fields*N_SF;
int IDX = (iphi*N2+jphi)*n_fields;
#define LOCAL_MATRIX(a,b) local_matrix[IDX+N2*a+b]
double dphiidphij = 0;
double dtau[DIMENSION];
double udtau = 0;
double dphiidtau = 0;
for (int j = 0; j < DIMENSION; ++j) {
dtau[j] = dphij[j]-phij*dc[j]/q;
dphiidphij += dphii[j]*dphij[j];
udtau += u[j]*dtau[j];
dphiidtau += dphii[j]*dtau[j];
/////////////
const double jw = QW[iqp]*detj;
for (int iphi = 0; iphi < N_SF; ++iphi) {
double phii = phi[iphi];
const double *dphii = dphi[iphi];
double tau[DIMENSION][DIMENSION];
for (int i = 0; i < DIMENSION; ++i)
for (int j = 0; j < DIMENSION; ++j)
tau[i][j] = du[i][j]-u[i]*dc[j]/q;
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) {
int U = j;
LOCAL_MATRIX(U,P) += -jw*dphii[j]*phij;
LOCAL_MATRIX(P,U) += jw*phii*dphij[j];
double utau =0;
double dphisigmadq = 0;
double dutaudq = 0;
double utau = 0;
double dphisigma=0;
for (int k = 0; k < DIMENSION; ++k) {
int V = k;
//utau += u[k]*dphii[k]*u[j]/c;
utau += u[k]*tau[j][k];
dutaudq += u[k]*u[j]*(dc[k]*phij/(q*q)-dcdq[k]/q);
dphisigmadq += 2*mu*dphii[k]*0.5*(u[k]*dc[j]*phij/(q*q)+u[j]*dc[k]*phij/(q*q)-u[k]*dcdq[j]/q-u[j]*dcdq[k]/q);
//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]/q
+jw*2*mu*dphii[k]*dtau[j]*0.5
+jw*rho*phii*u[j]*dphij[k]/q;
dphisigma += 2*mu*dphii[k]*0.5*(tau[k][j]+tau[j][k]);
}
local_vector[iphi+N_SF*j] += jw*(
dphisigma
//+rho*phii*dudt[j]-rho*utau
+rho*phii*(dudt[j]+u[j]/q*divu+utau/q)
-p*dphii[j]
);
}
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];
int P = DIMENSION;
//int N2 = N_SF*N_SF;
//int IDX = iphi*N2+jphi;
int N2 = n_fields*N_SF;
int IDX = (iphi*N2+jphi)*n_fields;
#define LOCAL_MATRIX(a,b) local_matrix[IDX+N2*a+b]
double dphiidphij = 0;
double dtau[DIMENSION];
double udtau = 0;
double dphiidtau = 0;
for (int j = 0; j < DIMENSION; ++j) {
dtau[j] = dphij[j]-phij*dc[j]/q;
dphiidphij += dphii[j]*dphij[j];
udtau += u[j]*dtau[j];
dphiidtau += dphii[j]*dtau[j];
}
for (int j = 0; j < DIMENSION; ++j) {
int U = j;
LOCAL_MATRIX(U,P) += -jw*dphii[j]*phij;
LOCAL_MATRIX(P,U) += jw*phii*dphij[j];
double utau =0;
double dphisigmadq = 0;
double dutaudq = 0;
for (int k = 0; k < DIMENSION; ++k) {
int V = k;
utau += u[k]*tau[j][k];
dutaudq += u[k]*u[j]*(dc[k]*phij/(q*q)-dcdq[k]/q);
dphisigmadq += 2*mu*dphii[k]*0.5*(u[k]*dc[j]*phij/(q*q)+u[j]*dc[k]*phij/(q*q)-u[k]*dcdq[j]/q-u[j]*dcdq[k]/q);
//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]/q
+jw*2*mu*dphii[k]*dtau[j]*0.5
+jw*rho*phii*u[j]*dphij[k]/q;
}
//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/q*divu+udtau/q);
LOCAL_MATRIX(U,Q) += jw*rho*phii*(-u[j]*divu*phij/(q*q)-utau*phij/(q*q)+dutaudq/q)+jw*dphisigmadq;
}
//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/q*divu+udtau/q);
LOCAL_MATRIX(U,Q) += jw*rho*phii*(-u[j]*divu*phij/(q*q)-utau*phij/(q*q)+dutaudq/q)+jw*dphisigmadq;
LOCAL_MATRIX(P,P) += jw*epsilon*dphiidphij;
LOCAL_MATRIX(Q,Q) += -jw*phii*phij;
LOCAL_MATRIX(P,Q) += jw*phii*phij/dt;
}
LOCAL_MATRIX(P,P) += jw*epsilon*dphiidphij;
LOCAL_MATRIX(Q,Q) += -jw*phii*phij;
LOCAL_MATRIX(P,Q) += jw*phii*phij/dt;
}
}
hxtLinearSystemAddToMatrix(lsys,iel,iel,local_matrix);
hxtLinearSystemAddToRhs(lsys,rhs,iel,local_vector);
}
hxtLinearSystemAddToMatrix(lsys,iel,iel,local_matrix);
hxtLinearSystemAddToRhs(lsys,rhs,iel,local_vector);
}
free(all_local_matrix);
free(all_local_vector);
......@@ -384,7 +390,7 @@ static void fluid_problem_assemble_system(FluidProblem *problem, double *rhs, co
static void apply_boundary_conditions(const Mesh *mesh, double *solution, int nBnd, StrongBoundary *bnds)
static void apply_boundary_conditions(const Mesh *mesh, int n_fields, double *solution, int nBnd, StrongBoundary *bnds)
{
for (int ibnd = 0; ibnd < nBnd; ++ibnd) {
StrongBoundary *bnd = bnds + ibnd;
......@@ -416,7 +422,8 @@ int fluid_problem_implicit_euler(FluidProblem *problem, double dt)
static double totaltime = 0;
struct timespec tic, toc;
const Mesh *mesh = problem->mesh;
apply_boundary_conditions(mesh, problem->solution, problem->n_strong_boundaries, problem->strong_boundaries);
apply_boundary_conditions(mesh, problem->nFields, problem->solution, problem->n_strong_boundaries, problem->strong_boundaries);
int n_fields = problem->nFields;
double *solution_old = malloc(mesh->n_nodes*n_fields*sizeof(double));
double *rhs = malloc(mesh->n_nodes*n_fields*sizeof(double));
for (int i=0; i<mesh->n_nodes*n_fields; ++i)
......@@ -482,7 +489,7 @@ static void fluid_problem_compute_porosity(FluidProblem *problem)
Mesh *mesh = problem->mesh;
double *volume = problem->particle_volume;
for (int i = 0; i < mesh->n_nodes; ++i){
problem->porosity[i] = 0;
problem->compacity[i] = 0;
}
for (int i = 0; i < problem->n_particles; ++i) {
if(problem->particle_element_id[i] == -1) continue;
......@@ -490,19 +497,11 @@ static void fluid_problem_compute_porosity(FluidProblem *problem)
shape_functions(&problem->particle_uvw[i*DIMENSION],sf);
const uint32_t *el = &mesh->elements[problem->particle_element_id[i]*N_N];
for (int j=0; j<N_N;++j)
problem->porosity[el[j]] += sf[j]*volume[i];
}
for (int i = 0; i < mesh->n_nodes; ++i){
if(problem->node_volume[i]==0.){
problem->porosity[i] = 1.;
}
else{
problem->porosity[i] = (1-problem->porosity[i]/problem->node_volume[i]);
}
problem->compacity[el[j]] += sf[j]*volume[i];
}
}
FluidProblem *fluid_problem_new(const char *mesh_file_name, double g, double mu, double rho, double epsilon, int n_strong_boundaries, StrongBoundary *strong_boundaries)
FluidProblem *fluid_problem_new(const char *mesh_file_name, double g, double *mu, double *rho, double epsilon, int n_fluids, int n_strong_boundaries, StrongBoundary *strong_boundaries)
{
static int initialized = 0;
if (!initialized) {
......@@ -517,9 +516,16 @@ FluidProblem *fluid_problem_new(const char *mesh_file_name, double g, double mu,
if (!mesh)
return NULL;
problem->mesh_tree = mesh_tree_create(mesh);
problem->rho = rho;
problem->mu = mu;
problem->rho = malloc(n_fluids*sizeof(double));
problem->mu = malloc(n_fluids*sizeof(double));
for (int i = 0; i < n_fluids; ++i)
{
problem->rho[i] = rho[i];
problem->mu[i] = mu[i];
}
problem->nFluids = n_fluids;
problem->nFields = (DIMENSION+1)*n_fluids+1;
int n_fields = problem->nFields;
problem->g = g;
problem->mesh = mesh;
problem->epsilon = epsilon;
......@@ -530,8 +536,8 @@ FluidProblem *fluid_problem_new(const char *mesh_file_name, double g, double mu,
problem->strong_boundaries[i].apply = strong_boundaries[i].apply;
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->compacity = malloc(mesh->n_nodes*sizeof(double));
problem->old_compacity = 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)
{
......@@ -565,8 +571,8 @@ FluidProblem *fluid_problem_new(const char *mesh_file_name, double g, double mu,
void fluid_problem_free(FluidProblem *problem)
{
free(problem->solution);
free(problem->porosity);
free(problem->old_porosity);
free(problem->compacity);
free(problem->old_compacity);
hxtLinearSystemDelete(&problem->linear_system);
free(problem->particle_uvw);
free(problem->particle_element_id);
......@@ -600,12 +606,13 @@ void fluid_problem_adapt_gen_mesh(FluidProblem *problem, double lcmax, double lc
double gradM = 0.;
double gradPM = 0.;
double volT = 0.;
double *porosity = problem->porosity;
int n_fields = nFields;
double *compacity = problem->compacity;
for (int iel = 0; iel < mesh->n_elements; ++iel) {
const unsigned int *el = &mesh->elements[iel*N_N];
double C[N_N],U[DIMENSION][N_N], P[N_N];
for (int i = 0; i < N_N; ++i){
C[i] = porosity[el[i]];
C[i] = compacity[el[i]];
P[i] = solution[el[i]*n_fields+DIMENSION];
for (int j=0; j<DIMENSION; ++j) {
U[j][i] = solution[el[i]*n_fields+j];
......@@ -761,6 +768,7 @@ void fluid_problem_adapt_mesh(FluidProblem *problem, double lcmax, double lcmin,
fluid_problem_adapt_gen_mesh(problem, lcmax, lcmin, n_el);
clock_get(&tic);
Mesh *new_mesh = mesh_load_msh("adapt/mesh.msh");
int n_fields = problem->nFields;
clock_get(&toc);
printf("Time mesh_load : %.6e \n",clock_diff(&tic,&toc));
double *new_solution = malloc(sizeof(double)*new_mesh->n_nodes*n_fields);
......@@ -801,14 +809,14 @@ void fluid_problem_adapt_mesh(FluidProblem *problem, double lcmax, double lcmin,
free(newx);
free(problem->solution);
problem->solution = new_solution;
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);
free(problem->compacity);
problem->compacity = malloc(sizeof(double)*new_mesh->n_nodes);
free(problem->old_compacity);
problem->old_compacity = 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.;
problem->compacity[i] = 1.;
problem->old_compacity[i] = 1.;
}
mesh_free(problem->mesh);
problem->mesh = new_mesh;
......@@ -850,7 +858,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, int init)
{
struct timespec tic, toc;
int n_fields = problem->nFields;
if(problem->n_particles != n) {
problem->n_particles = n;
......@@ -879,9 +887,6 @@ 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));
clock_get(&tic);
for (int i = 0; i < n; ++i) {
problem->particle_mass[i] = mass[i];
problem->particle_volume[i] = volume[i];
......@@ -890,20 +895,16 @@ void fluid_problem_set_particles(FluidProblem *problem, int n, double *mass, dou
problem->particle_velocity[i*DIMENSION+k] = velocity[i*DIMENSION+k];
}
}
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];
problem->old_compacity[i] = problem->compacity[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];
problem->solution[i*n_fields+Q] = 1-problem->compacity[i];
}
clock_get(&tic);
fluid_problem_compute_porosity(problem);
clock_get(&toc);
printf("Time compute_porosity : %.6e \n",clock_diff(&tic,&toc));
}
int *fluid_problem_particle_element_id(FluidProblem *problem)
......
......@@ -42,14 +42,16 @@ typedef struct {
typedef struct {
double epsilon;
double rho;
double mu;
double *rho;
double *mu;
int nFluids;
int nFields;
double alpha;
double g;
Mesh *mesh;
MeshTree *mesh_tree;
double *porosity;
double *old_porosity;
double *compacity;
double *old_compacity;
double *solution;
double *solution_explicit;
double *node_volume;
......
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