Commit 77c0aa93 authored by Michel Henry's avatar Michel Henry
Browse files

ALE

parent ec2c77c0
......@@ -173,7 +173,7 @@ void fluid_problem_node_force_volume(FluidProblem *problem, const double *soluti
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 c, const double *dc, const double rho, double mu, const double dt, int eid, const double *data, double *f00, double *f01)
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)
{
const int vid = wbnd->vid;
const int pid = wbnd->pid;
......@@ -185,11 +185,13 @@ static void f_boundary(WeakBoundary *wbnd, FluidProblem *problem,const double *n
double s_c = 0;
double unold = 0;
double unext = 0;
double unmesh = 0;
double dcn = 0;
for (int iD = 0; iD < D; ++iD) {
u[iD] = s[U+iD];
dp[iD] = ds[P*D+iD];
unold += sold[U+iD]*n[iD];
unmesh += mesh_velocity[iD]*n[iD];
unext += (vid<0?s[U+iD]:data[vid+iD])*n[iD];
dcn += dc[iD]*n[iD];
s_c += vid<0?0:(u[iD]-data[vid+iD])*n[iD];
......@@ -230,9 +232,9 @@ static void f_boundary(WeakBoundary *wbnd, FluidProblem *problem,const double *n
double unbnd = unold;
if (wbnd->type == BND_WALL) unbnd = unold/2;
else if(vid>0) unbnd = (unold+unext)/2;
if (unbnd<0) {
if (unbnd - unmesh < 0) {
for (int id = 0; id < D; ++id) {
f0[U+id] += (unbnd*(vid<0?0:data[vid+id])-unold*u[id])*rho/c;
f0[U+id] += ((unbnd - unmesh)*(vid<0?0:data[vid+id])-(unold - unmesh)*u[id])*rho/c;
f00[(U+id)*n_fields+U+id] -= unold*rho/c;
}
}
......@@ -368,12 +370,12 @@ static void particle_force_f(FluidProblem *problem, double *f, double *dfdu, dou
*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 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) {
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];
double tauc = problem->tauc[iel];
const int n_fields = fluid_problem_n_fields(problem);
double u[D], uold[D], du[D][D], duold[D][D],dp[D];
double u[D], uold[D], du[D][D], duold[D][D], umesh[D], dp[D];
if(problem->stab_param != 0)
taup = tauc = 0;
double divu = 0;
......@@ -382,6 +384,7 @@ static void fluid_problem_f(const FluidProblem *problem, const double *sol, doub
u[iD] = sol[U+iD];
uold[iD] = solold[U+iD];
nold += uold[iD]*uold[iD];
umesh[iD] = mesh_velocity[iD];
dp[iD] = dsol[P*D+iD];
for (int jD = 0; jD < D; ++jD) {
du[iD][jD] = dsol[(U+iD)*D+jD];
......@@ -414,8 +417,10 @@ static void fluid_problem_f(const FluidProblem *problem, const double *sol, doub
// = ui/c * (dj(uj)- uj/c*dj(c)) + uj/c*dj(ui)
for (int j = 0; j < D; ++j) {
f0[U+i] += rho/c*(uold[j]*du[i][j] + u[i]*(duold[j][j]-uold[j]/c*dc[j]));
f0[U+i] -= rho/c*(umesh[j]*du[i][j]);
f00[(U+i)*n_fields+U+i] += rho/c*(duold[j][j]-uold[j]/c*dc[j]);
f01[(U+i)*n_fields*D+(U+i)*D+j] += rho/c*uold[j];
f01[(U+i)*n_fields*D+(U+i)*D+j] -= rho/c*umesh[j];
}
}
for (int j = 0; j < D; ++j) {
......@@ -518,7 +523,8 @@ static void compute_weak_boundary_conditions(FluidProblem *problem, const double
}
for (int iqp = 0; iqp < N_LQP; ++iqp) {
double *dataqp = &data[n_value*(N_LQP*iinterface+iqp)];
double dp[D]={0};
double dp[D] = {0};
double umesh[D] = {0};
double *local_vector = &all_local_vector[local_size*iel];
double *local_matrix = &all_local_matrix[local_size*local_size*iel];
double phi[DIMENSION];
......@@ -547,6 +553,10 @@ static void compute_weak_boundary_conditions(FluidProblem *problem, const double
s[j] += phi[i]*dof;
sold[j] += phi[i]*dofold;
}
for(int j = 0; j < D; ++j){
double dofmesh = problem->mesh_velocity[nodes[i]*D+j];
umesh[j] += phi[i]*dofmesh;
}
}
double rho, mu;
fluid_problem_interpolate_rho_and_nu(problem,a, &rho, &mu);
......@@ -563,7 +573,7 @@ static void compute_weak_boundary_conditions(FluidProblem *problem, const double
}
}
}
f_boundary(wbnd,problem,n,f0,s,ds,sold,dsold,c,dc,rho,mu,dt,iel,dataqp,f00,f01);
f_boundary(wbnd,problem,n,f0,s,ds,sold,dsold,umesh,c,dc,rho,mu,dt,iel,dataqp,f00,f01);
for (int ifield = 0; ifield < n_fields; ++ifield) {
for (int iphi = 0; iphi < D; ++iphi){
if (ifield<D){
......@@ -751,6 +761,7 @@ static void fluid_problem_volume(FluidProblem *problem, const double *solution_o
double *sold = malloc(sizeof(double)*n_fields);
double *ds = malloc(sizeof(double)*n_fields*D);
double *dsold = malloc(sizeof(double)*n_fields*D);
double *mesh_velocity = malloc(sizeof(double)*D);
double *f0 = malloc(sizeof(double)*n_fields);
double *f1 = malloc(sizeof(double)*n_fields*D);
double *f00 = malloc(sizeof(double)*n_fields*n_fields);
......@@ -796,6 +807,7 @@ static void fluid_problem_volume(FluidProblem *problem, const double *solution_o
double bf[D];
for (int j = 0; j < D; ++j) {
bf[j] = 0;
mesh_velocity[j] = 0;
}
for (int i = 0; i < n_fields; ++i) {
s[i] = 0;
......@@ -826,6 +838,7 @@ static void fluid_problem_volume(FluidProblem *problem, const double *solution_o
a += problem->concentration[iel*N_N+i]*phi[i];
}
for (int j=0; j<D; ++j) {
mesh_velocity[j] += phi[i]*problem->mesh_velocity[el[i]*D+j];
bf[j] += phi[i]*problem->bulk_force[el[i]*D+j];
}
}
......@@ -835,7 +848,7 @@ static void fluid_problem_volume(FluidProblem *problem, const double *solution_o
const double jw = QW[iqp]*detj;
problem->kinetic_energy += rho*sqrt(s[U]*s[U]+s[U+1]*s[U+1])*sqrt(s[U]*s[U]+s[U+1]*s[U+1])/2*jw;
fluid_problem_f(problem,s,ds,sold,dsold,c,dc,bf,cold,rho,mu,dt,iel,f0,f1,f00,f10,f01,f11);
fluid_problem_f(problem,s,ds,sold,dsold,mesh_velocity,c,dc,bf,cold,rho,mu,dt,iel,f0,f1,f00,f10,f01,f11);
for (int ifield = 0; ifield < n_fields; ++ifield) {
for (int iphi = 0; iphi < N_SF; ++iphi){
local_vector[iphi+N_SF*ifield] += phi[iphi]*f0[ifield]*jw;
......@@ -869,6 +882,7 @@ static void fluid_problem_volume(FluidProblem *problem, const double *solution_o
free(ds);
free(dsold);
free(sold);
free(mesh_velocity);
free(f0);
free(f1);
free(f00);
......@@ -1309,6 +1323,7 @@ FluidProblem *fluid_problem_new(double *g, int n_fluids, double *mu, double *rho
problem->porosity = NULL;
problem->oldporosity = NULL;
problem->solution = NULL;
problem->mesh_velocity = NULL;
problem->concentration = NULL;
problem->n_particles = 0;
problem->particle_uvw = NULL;
......@@ -1338,6 +1353,7 @@ void fluid_problem_free(FluidProblem *problem)
free(problem->mu);
free(problem->rho);
free(problem->solution);
free(problem->mesh_velocity);
free(problem->concentration);
free(problem->porosity);
free(problem->bulk_force);
......@@ -1636,6 +1652,9 @@ void fluid_problem_set_mesh(FluidProblem *problem, Mesh *mesh) {
free(problem->concentration);
problem->concentration = (double*)malloc(mesh->n_elements*N_N*sizeof(double));
}
free(problem->mesh_velocity);
problem->mesh_velocity = (double*)calloc(mesh->n_nodes*D,sizeof(double));
// for(int i = 0; i < mesh->n_nodes*D; ++i) problem->mesh_velocity[i] = 0;
problem->taup = malloc(sizeof(double)*mesh->n_elements);
for (int i = 0; i < problem->mesh->n_elements; ++i) problem->taup[i] = 0;
for (int i = 0; i < problem->mesh->n_nodes*n_fields; ++i){
......@@ -1867,6 +1886,10 @@ double *fluid_problem_solution(const FluidProblem *p)
return p->solution;
}
double *fluid_problem_mesh_velocity(const FluidProblem *p){
return p->mesh_velocity;
}
double *fluid_problem_concentration_dg_grad(const FluidProblem *p)
{
return p->grad_a_cg;
......
......@@ -69,6 +69,7 @@ struct FluidProblem {
double *oldporosity;
double *concentration;
double *solution;
double *mesh_velocity;
double *node_volume;
double *element_size;
double volume_drag;
......@@ -117,6 +118,7 @@ int fluid_problem_n_particles(FluidProblem *problem);
void fluid_problem_after_import(FluidProblem *problem);
int fluid_problem_n_fields(const FluidProblem *p);
double *fluid_problem_solution(const FluidProblem *p);
double *fluid_problem_mesh_velocity(const FluidProblem *p);
void fluid_problem_set_weak_boundary(FluidProblem *p, const char *tag, BoundaryType type, BoundaryCallback *cb, int vid, int pid, int cid, int aid);
void fluid_problem_set_strong_boundary(FluidProblem *p, const char *tag, int field, int row, BoundaryCallback *apply);
double *fluid_problem_old_porosity(const FluidProblem *p);
......
......@@ -484,6 +484,10 @@ class FluidProblem :
def velocity(self) :
"""Gives access to the fluid velocity value at the mesh nodes."""
return self.solution()[:,:-1]
def mesh_velocity(self) :
"""Give access to the mesh velocity value at the mesh nodes."""
return self._get_matrix("mesh_velocity", self.n_nodes(), self._dim)
def pressure(self) :
"""Gives access to the fluid field value at the mesh nodes."""
......
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