Commit 5d41ddd1 authored by Matthieu Constant's avatar Matthieu Constant Committed by Jonathan Lambrechts
Browse files

Discontinuous stabilization

parent 9bf7d9d2
......@@ -15,11 +15,11 @@ class StrongBoundary(Structure):
class fluid_problem :
def __init__(self, mesh_file_name, mu, rho, epsilon, strong_boundaries):
def __init__(self, mesh_file_name, mu, rho, alpha, autoEpsilon, strong_boundaries):
nsb = len(strong_boundaries)
asb = (StrongBoundary*nsb)(*[StrongBoundary(i[0].encode(),i[1],i[2]) for i in strong_boundaries])
lib.fluid_problem_new.restype = c_void_p
self._fp = c_void_p(lib.fluid_problem_new(mesh_file_name.encode(), c_double(mu), c_double(rho), c_double(epsilon), nsb, asb))
self._fp = c_void_p(lib.fluid_problem_new(mesh_file_name.encode(), c_double(mu), c_double(rho), c_double(alpha), c_bool(autoEpsilon), nsb, asb))
if self._fp == None :
raise NameError("cannot create fluid problem")
......@@ -27,8 +27,8 @@ class fluid_problem :
if(self._fp != None) :
lib.fluid_problem_free(self._fp)
def adapt_mesh(self, gradmin, gradmax, lcmin) :
lib.fluid_problem_adapt_mesh(self._fp, c_double(gradmin), c_double(gradmax), c_double(lcmin))
def adapt_mesh(self, gradmin, gradmax, lcmin, autoEpsilon) :
lib.fluid_problem_adapt_mesh(self._fp, c_double(gradmin), c_double(gradmax), c_double(lcmin), c_bool(autoEpsilon))
def export(self, output_dir, t, it) :
lib.fluid_problem_export(self._fp, output_dir.encode(), c_double(t), c_int(it))
......
#include "tools.h"
#include <stdlib.h>
#include <stdio.h>
#include <stdbool.h>
#include <float.h>
#include <math.h>
#include "fluid_problem.h"
......@@ -136,7 +137,7 @@ static void fluid_problem_assemble_system(FluidProblem *problem, const double *s
double *local_matrix = malloc(sizeof(double)*n_fields*n_fields*9);
const double mu = problem->mu;
const double rho = problem->rho;
const double epsilon = problem->epsilon;
const double *epsilon = problem->epsilon;
for (int iel=0; iel < mesh->n_triangles; ++iel) {
const int *tri = &mesh->triangles[iel*3];
const double x[3][2] = {
......@@ -216,7 +217,7 @@ static void fluid_problem_assemble_system(FluidProblem *problem, const double *s
);
local_vector[iphi+6] += jw*(
epsilon*(dphii[0]*dp[0]+dphii[1]*dp[1])
epsilon[iel]*(dphii[0]*dp[0]+dphii[1]*dp[1])
+phii*((du[0]+dv[1])+dcdt)
);
for (int jphi = 0; jphi < 3; ++jphi) {
......@@ -235,7 +236,7 @@ static void fluid_problem_assemble_system(FluidProblem *problem, const double *s
local_matrix[(iphi+3)*9+jphi+6] +=-jw*dphii[1]*phij; // V P
local_matrix[(iphi+6)*9+jphi+0] += jw*phii*dphij[0]; // P U
local_matrix[(iphi+6)*9+jphi+3] += jw*phii*dphij[1]; // P V
local_matrix[(iphi+6)*9+jphi+6] += jw*epsilon*(dphii[0]*dphij[0]+dphii[1]*dphij[1]);// P P
local_matrix[(iphi+6)*9+jphi+6] += jw*epsilon[iel]*(dphii[0]*dphij[0]+dphii[1]*dphij[1]);// P P
}
}
}
......@@ -321,13 +322,54 @@ static void fluid_problem_compute_node_volume(FluidProblem *problem) {
}
}
FluidProblem *fluid_problem_new(const char *mesh_file_name, double mu, double rho, double epsilon, int n_strong_boundaries, StrongBoundary *strong_boundaries)
static void computeEpsilon(FluidProblem *problem, bool autoEpsilon)
{
double alpha = problem->alpha;
double nu = problem->mu/problem->rho;
const Mesh *mesh = problem->mesh;
double a,b,c,h;
double epsMax = 0;
double epsMin = 10000;
for(int iel = 0; iel < mesh->n_triangles; ++iel)
{
if (autoEpsilon)
{
const int *tri = &mesh->triangles[iel*3];
const double x[3][2] = {
{mesh->x[tri[0]*3+0],mesh->x[tri[0]*3+1]},
{mesh->x[tri[1]*3+0],mesh->x[tri[1]*3+1]},
{mesh->x[tri[2]*3+0],mesh->x[tri[2]*3+1]}};
a = sqrt(pow((x[0][0]-x[1][0]),2)+pow((x[0][1]-x[1][1]),2));
b = sqrt(pow((x[1][0]-x[2][0]),2)+pow((x[1][1]-x[2][1]),2));
c = sqrt(pow((x[2][0]-x[0][0]),2)+pow((x[2][1]-x[0][1]),2));
h = a*b*c/sqrt((a+b+c)*(b+c-a)*(c+a-b)*(a+b-c));
problem->epsilon[iel] = alpha*h*h/nu;
if (problem->epsilon[iel] > epsMax)
{
epsMax = problem->epsilon[iel];
}
if (problem->epsilon[iel] < epsMin)
{
epsMin = problem->epsilon[iel];
}
}
else
{
problem->epsilon[iel] = alpha;
}
}
printf("Epsilon Max = %g\n",epsMax);
printf("Epsilon Min = %g\n",epsMin);
}
FluidProblem *fluid_problem_new(const char *mesh_file_name, double mu, double rho, double alpha, bool autoEpsilon, int n_strong_boundaries, StrongBoundary *strong_boundaries)
{
FluidProblem *problem = malloc(sizeof(FluidProblem));
problem->epsilon = epsilon;
problem->rho = rho;
problem->alpha = alpha;
problem->mu = mu;
Mesh *mesh = mesh_load_msh(mesh_file_name);
problem->epsilon = malloc(sizeof(double)*mesh->n_triangles);
if (!mesh)
return NULL;
problem->mesh = mesh;
......@@ -364,6 +406,7 @@ FluidProblem *fluid_problem_new(const char *mesh_file_name, double mu, double rh
problem->particle_mass = malloc(0);
problem->particle_volume = malloc(0);
problem->particle_velocity = malloc(0);
computeEpsilon(problem, autoEpsilon);
return problem;
}
......@@ -373,6 +416,7 @@ void fluid_problem_free(FluidProblem *problem)
free(problem->node_force);
free(problem->porosity);
free(problem->old_porosity);
free(problem->epsilon);
linear_system_free(problem->linear_system);
free(problem->particle_uvw);
free(problem->particle_element_id);
......@@ -406,7 +450,7 @@ static void fluid_problem_compute_porosity(FluidProblem *problem)
}
}
void fluid_problem_adapt_mesh(FluidProblem *problem, double gradmin, double gradmax, double lcmin)
void fluid_problem_adapt_mesh(FluidProblem *problem, double gradmin, double gradmax, double lcmin, bool autoEpsilon)
{
Mesh *mesh = problem->mesh;
double *solution = problem->solution;
......@@ -509,9 +553,12 @@ void fluid_problem_adapt_mesh(FluidProblem *problem, double gradmin, double grad
problem->node_force = malloc(sizeof(double)*2*new_mesh->n_nodes);
mesh_free(problem->mesh);
problem->mesh = new_mesh;
free(problem->epsilon);
problem->epsilon = malloc(sizeof(double)*new_mesh->n_triangles);
mesh_particles_to_mesh(problem->mesh, problem->n_particles, problem->particle_position, problem->particle_element_id, problem->particle_uvw);
fluid_problem_compute_node_volume(problem);
fluid_problem_compute_porosity(problem);
computeEpsilon(problem, autoEpsilon);
linear_system_free(problem->linear_system);
problem->linear_system = linear_system_new(problem->mesh, n_fields, problem->n_strong_boundaries, problem->strong_boundaries);
}
......
......@@ -5,9 +5,10 @@
#include "linear_system.h"
typedef struct {
double epsilon;
double *epsilon;
double rho;
double mu;
double alpha;
Mesh *mesh;
double *porosity;
double *old_porosity;
......@@ -31,8 +32,8 @@ int fluid_problem_export(const FluidProblem *problem, const char *output_dir, do
// complete force on the particle (including gravity)
void fluid_problem_compute_node_particle_force(FluidProblem *problem, double dt, double *particle_force);
void fluid_problem_implicit_euler(FluidProblem *problem, double dt);
FluidProblem *fluid_problem_new(const char *mesh_file_name, double mu, double rho, double epsilon, int n_strong_boundaries, StrongBoundary *strong_boundaries);
FluidProblem *fluid_problem_new(const char *mesh_file_name, double mu, double rho, double alpha, bool autoEpsilon, int n_strong_boundaries, StrongBoundary *strong_boundaries);
void fluid_problem_set_particles(FluidProblem *problem, int n, double *mass, double *volume, double *position, double *velocity);
void fluid_problem_free(FluidProblem *problem);
void fluid_problem_adapt_mesh(FluidProblem *problem, double gradmin, double gradmax, double lcmin);
void fluid_problem_adapt_mesh(FluidProblem *problem, double gradmin, double gradmax, double lcmin, bool autoEpsilon);
#endif
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