Commit b717baec authored by Jonathan Lambrechts's avatar Jonathan Lambrechts
Browse files

weak bc

parent 73a154c7
......@@ -32,7 +32,7 @@ signal.signal(signal.SIGINT, signal.SIG_DFL)
BNDCB = CFUNCTYPE(None,c_int,POINTER(c_double),POINTER(c_double))
class StrongBoundary(Structure):
_fields_ = [("tag",c_char_p),("field",c_int),("apply",BNDCB)]
_fields_ = [("tag",c_char_p),("row",c_int),("field",c_int),("apply",BNDCB)]
class fluid_problem :
......@@ -53,9 +53,7 @@ class fluid_problem :
r = c_void_p(tmp.ctypes.data)
r.tmp = tmp
return r
asb = (StrongBoundary*nsb)(*[StrongBoundary(i[0].encode(),i[1],BNDCB(Bnd(i[2]).apply)) for i in strong_boundaries])
asb = (StrongBoundary*nsb)(*[StrongBoundary(i[0].encode(),i[1],i[2],BNDCB(Bnd(i[3]).apply)) for i in strong_boundaries])
self.asb = asb
lib.fluid_problem_new.restype = c_void_p
self._fp = c_void_p(lib.fluid_problem_new(mesh_file_name.encode(), c_double(g), np2c(mu), np2c(rho), c_double(epsilon), nsb, asb, c_int(noGrains), c_int(n_fluids)))
......
......@@ -281,7 +281,56 @@ static void f_q_1(double f[D], const double dp[D], const double epsilon)
}
}
static void compute_weak_boundary_conditions(FluidProblem *problem, double dt, double *all_local_vector, double *all_local_matrix)
{
const Mesh *mesh = problem->mesh;
const double *solution = problem->solution;
const int n_fields = problem_n_fields(problem);
const size_t local_size = N_SF*n_fields;
for (int iphys = 0; iphys < mesh->n_phys; ++iphys) {
int forced = 0;
if (strcmp(mesh->phys_name[iphys],"Right") == 0)
forced = 1;
for (int ibnd = 0; ibnd < mesh->phys_n_boundaries[iphys]; ++ibnd) {
const int *bnd = &mesh->phys_boundaries[iphys][ibnd*2];
const int iel = bnd[0];
const int i0 = bnd[1];
const int i1 = (bnd[1]+1) % 3;
const int *el = &mesh->elements[iel*N_N];
const int nodes[2] = {el[i0],el[i1]};
const double *x[2] = {&mesh->x[nodes[0]*3], &mesh->x[nodes[1]*3]};
const double dx[2] = {x[1][0]-x[0][0], x[1][1]-x[0][1]};
const double l = sqrt(dx[0]*dx[0] + dx[1]*dx[1]);
const double n[2] = {-dx[1]/l,dx[0]/l};
double *local_matrix = &all_local_matrix[local_size*local_size*iel];
double *local_vector = &all_local_vector[local_size*iel];
const int U = 0;
const int Q = 2;
const int P = 3;
double p0 = solution[el[i0]*n_fields+P];
double p1 = solution[el[i1]*n_fields+P];
/*for (int k = 0; k < D; ++k) {
if (!forced) {
local_matrix[(i0*n_fields+U+k)*local_size + (i0*n_fields+P)] += l/3*n[k];
local_matrix[(i0*n_fields+U+k)*local_size + (i1*n_fields+P)] += l/6*n[k];
local_matrix[(i1*n_fields+U+k)*local_size + (i0*n_fields+P)] += l/6*n[k];
local_matrix[(i1*n_fields+U+k)*local_size + (i1*n_fields+P)] += l/3*n[k];
local_vector[i0*n_fields+U+k] += l*(p0/3+p1/6)*n[k];
local_vector[i1*n_fields+U+k] += l*(p0/6+p1/3)*n[k];
}
}*/
if (forced) {
double a = -1;
local_vector[i0+Q*N_SF] += l*(p0/3+p1/6)*a;
local_vector[i1+Q*N_SF] += l*(p0/6+p1/3)*a;
local_matrix[(i0*n_fields+Q)*local_size + (i0*n_fields+P)] += l/3*a;
local_matrix[(i0*n_fields+Q)*local_size + (i1*n_fields+P)] += l/6*a;
local_matrix[(i1*n_fields+Q)*local_size + (i0*n_fields+P)] += l/6*a;
local_matrix[(i1*n_fields+Q)*local_size + (i1*n_fields+P)] += l/3*a;
}
}
}
}
static void fluid_problem_assemble_system(FluidProblem *problem, double *rhs, const double *solution_old, double dt)
{
......@@ -304,7 +353,26 @@ static void fluid_problem_assemble_system(FluidProblem *problem, double *rhs, co
all_local_vector[i] = 0;
for (size_t i = 0; i < N_E*local_size*local_size; ++i)
all_local_matrix[i] = 0;
char *forced_row = malloc(sizeof(char)*n_fields*mesh->n_nodes);
for (int i = 0; i < n_fields*mesh->n_nodes; ++i)
forced_row[i] = -1;
for (int ibnd = 0; ibnd < problem->n_strong_boundaries; ++ibnd) {
const StrongBoundary *bnd = problem->strong_boundaries + ibnd;
int iphys;
for (iphys = 0; iphys < mesh->n_phys; ++iphys) {
if (strcmp(bnd->tag, mesh->phys_name[iphys]) == 0)
break;
}
if (iphys == mesh->n_phys)
printf("Boundary tag \"%s\" not found.", bnd->tag);
for (int i = 0; i < mesh->phys_n_nodes[iphys]; ++i){
forced_row[mesh->phys_nodes[iphys][i]*n_fields + bnd->row] = bnd->field;
}
}
//compute_node_force_implicit(problem, dt, all_local_vector, all_local_matrix);
compute_weak_boundary_conditions(problem, dt, all_local_vector, all_local_matrix);
for (int iel=0; iel < mesh->n_elements; ++iel) {
const unsigned int *el = &mesh->elements[iel*N_N];
double dxdxi[D][D], dxidx[D][D], dphi[N_N][D];
......@@ -510,26 +578,27 @@ static void fluid_problem_assemble_system(FluidProblem *problem, double *rhs, co
}
}
}
for (int inode = 0; inode < N_SF; ++inode) {
for(int ifield = 0; ifield < n_fields; ++ifield) {
char forced_field = forced_row[el[inode]*n_fields+ifield];
if (forced_field == -1) continue;
int i = inode*n_fields + ifield;
for (int jnode = 0; jnode< N_SF; ++jnode){
for (int jfield= 0; jfield< n_fields; ++jfield){
local_matrix[(inode*n_fields+ifield)*local_size+(jnode*n_fields+jfield)] =
(inode==jnode && jfield == forced_field) ? 1. : 0.;
}
}
local_vector[inode+ifield*N_SF] = 0;
}
}
hxtLinearSystemAddToMatrix(lsys,iel,iel,local_matrix);
hxtLinearSystemAddToRhs(lsys,rhs,iel,local_vector);
}
free(all_local_matrix);
free(all_local_vector);
for (int ibnd = 0; ibnd < problem->n_strong_boundaries; ++ibnd) {
const StrongBoundary *bnd = problem->strong_boundaries + ibnd;
int iphys;
for (iphys = 0; iphys < mesh->n_phys; ++iphys) {
if (strcmp(bnd->tag, mesh->phys_name[iphys]) == 0)
break;
}
if (iphys == mesh->n_phys)
printf("Boundary tag \"%s\" not found.", bnd->tag);
for (int i = 0; i < mesh->phys_n_nodes[iphys]; ++i){
hxtLinearSystemSetMatrixRowIdentity(lsys, mesh->phys_nodes[iphys][i], bnd->field);
hxtLinearSystemSetRhsEntry(lsys, rhs,mesh->phys_nodes[iphys][i], bnd->field, 0.);
}
}
free(forced_row);
}
static void apply_boundary_conditions(const Mesh *mesh, double *solution, int n_fields, int nBnd, StrongBoundary *bnds)
......@@ -663,7 +732,7 @@ FluidProblem *fluid_problem_new(const char *mesh_file_name, double g, double *mu
hxtInitializeLinearSystems(0, NULL);
initialized = 1;
#ifdef HXT_HAVE_PETSC
hxtPETScInsertOptions("-pc_type lu -ksp_max_it 30", "fluid");
hxtPETScInsertOptions("-pc_type ilu -ksp_max_it 30", NULL);
#endif
}
FluidProblem *problem = (FluidProblem*)malloc(sizeof(FluidProblem));
......@@ -689,6 +758,7 @@ FluidProblem *fluid_problem_new(const char *mesh_file_name, double g, double *mu
problem->strong_boundaries[i].tag = strdup(strong_boundaries[i].tag);
problem->strong_boundaries[i].apply = strong_boundaries[i].apply;
problem->strong_boundaries[i].field = strong_boundaries[i].field;
problem->strong_boundaries[i].row = strong_boundaries[i].row;
}
problem->compacity = (double*)malloc(mesh->n_nodes*sizeof(double));
problem->old_compacity = (double*)malloc(mesh->n_nodes*sizeof(double));
......
......@@ -36,6 +36,7 @@
typedef struct {
char *tag;
int row;
int field;
void (*apply)(int n, const double *x, double *bnd);
} StrongBoundary;
......
......@@ -85,6 +85,18 @@ static int int_cmp(const void *p0, const void *p1)
return 1;
}
typedef struct {
int tag;
int node[2];
} BoundarySegment;
int segment_cmp(const void *ps0, const void *ps1) {
BoundarySegment *s0 = (BoundarySegment*) ps0;
BoundarySegment *s1 = (BoundarySegment*) ps1;
if (s0->node[0] != s1->node[0]) return s0->node[0]-s1->node[0];
return s0->node[1]-s1->node[1];
}
Mesh *mesh_load_msh(const char *filename)
{
Mesh *m = malloc(sizeof(Mesh));
......@@ -133,6 +145,7 @@ Mesh *mesh_load_msh(const char *filename)
return NULL;
}
m->x = malloc(sizeof(double)*m->n_nodes*3);
for (int i = 0; i < m->n_nodes; ++i) {
int nid;
if (fscanf(f, "%d %le %le %le", &nid, &m->x[i*3], &m->x[i*3+1], &m->x[i*3+2]) != 4){
......@@ -151,6 +164,8 @@ Mesh *mesh_load_msh(const char *filename)
printf("Cannot read elements.\n");
return NULL;
}
BoundarySegment *segments = malloc(sizeof(BoundarySegment)*n_el);
size_t n_segments = 0;
m->elements = malloc(sizeof(int)*(DIMENSION+1)*n_el);
m->n_elements = 0;
for (int i = 0; i < n_el; ++i) {
......@@ -175,6 +190,10 @@ Mesh *mesh_load_msh(const char *filename)
printf("Cannot read elements.\n");
return NULL;
}
segments[n_segments].tag = phys;
segments[n_segments].node[0] = (n[0] < n[1] ? n[0] : n[1]) -1;
segments[n_segments].node[1] = (n[0] < n[1] ? n[1] : n[0]) -1;
n_segments++;
mesh_add_physical_node(m, 1, phys, n[0]-1);
mesh_add_physical_node(m, 1, phys, n[1]-1);
break;
......@@ -230,6 +249,53 @@ Mesh *mesh_load_msh(const char *filename)
}
m->phys_n_nodes[i] = l;
}
// boundary segments
// assign (eid,side) to segment and sort by tag;
// count by tag
m->phys_n_boundaries = malloc(m->n_phys*sizeof(int));
m->phys_boundaries = malloc(m->n_phys*sizeof(int*));
for (int j = 0; j < m->n_phys; ++j)
m->phys_n_boundaries[j] = 0;
for (int j = 0; j < n_segments; ++j) {
int i = 0;
int id = segments[j].tag;
segments[j].tag = -1;
for (i = 0; i < m->n_phys; ++i) {
if(m->phys_id[i] == id && m->phys_dimension[i] == 1) {
segments[j].tag = i;
break;
}
}
m->phys_n_boundaries[segments[j].tag]++;
}
for (int j = 0; j < m->n_phys; ++j){
m->phys_boundaries[j] = malloc(sizeof(int)*2*m->phys_n_boundaries[j]);
m->phys_n_boundaries[j] = 0;
}
qsort(segments,n_segments, sizeof(BoundarySegment), segment_cmp);
for (int i = 0; i < m->n_elements; ++i) {
for (int j = 0; j < 3; j++) {
int v0 = m->elements[i*3+j];
int v1 = m->elements[i*3+(j+1)%3];
BoundarySegment key;
if (v0 > v1) {
key.node[0] = v1;
key.node[1] = v0;
}
else {
key.node[1] = v1;
key.node[0] = v0;
}
BoundarySegment *found = bsearch(&key, segments, n_segments, sizeof(BoundarySegment), segment_cmp);
if (found == NULL) continue;
int tag = found->tag;
if (tag == -1) continue;
m->phys_boundaries[tag][m->phys_n_boundaries[tag]*2+0] = i;
m->phys_boundaries[tag][m->phys_n_boundaries[tag]*2+1] = j;
m->phys_n_boundaries[tag]++;
}
}
free(segments);
return m;
}
......@@ -261,12 +327,15 @@ void mesh_free(Mesh *m)
{
free(m->elements);
free(m->x);
for (int i = 0; i < m->n_phys; ++i){
for (int i = 0; i < m->n_phys; ++i) {
free(m->phys_name[i]);
free(m->phys_nodes[i]);
free(m->phys_boundaries[i]);
}
free(m->phys_name);
free(m->phys_nodes);
free(m->phys_n_boundaries);
free(m->phys_boundaries);
free(m->phys_dimension);
free(m->phys_id);
free(m->phys_n_nodes);
......
......@@ -36,6 +36,8 @@ typedef struct {
char **phys_name;
int *phys_n_nodes;
int **phys_nodes;
int *phys_n_boundaries;
int **phys_boundaries; // eid, side
int *phys_dimension;
int *phys_id;
} Mesh;
......
......@@ -48,7 +48,7 @@ tEnd = 100000 # final time
#numerical parameters
lcmin = .1 # mesh size
dt = 5 # time step
dt = 0.5 # time step
alpha = 1e-4 # stabilization coefficient
epsilon = alpha*lcmin**2 /nu # stabilization parametre
print('epsilon',epsilon)
......@@ -59,7 +59,7 @@ outf = 10 # number of iterations between o
#Object fluid creation + Boundary condition of the fluid (field 0 is horizontal velocity; field 1 is vertical velocity; field 2 is pressure)
#Format: strong_boundaries = [(Boundary tag, Fluid field, Value)
strong_boundaries = [("Top",1,0.),("Bottom",1,0.),("Left",1,0),("Right",1,0),("Top",0,0),("Bottom",0,0),("Right",3,0),("Left",0,lambda x : 1/(20*mu)*x[:,1]*(1-x[:, 1]))]
strong_boundaries = [("Top",1,1,0.),("Bottom",1,1,0.),("Left",1,1,0.),("Right",1,1,0.),("Top",0,0,0.),("Bottom",0,0,0.),("Left",0,0,lambda x : 1/(20*mu)*x[:,1]*(1-x[:, 1]))]
fluid = fluid.fluid_problem("mesh.msh",g,[nu*rho],[rho],epsilon,strong_boundaries,1,1)
ii = 0
t = 0
......
......@@ -79,4 +79,4 @@ while t < tEnd :
fluid.export_vtk(outputdir, t, ioutput)
ii += 1
print("%i : %.2g/%.2g (cpu %.6g)" % (ii, t, tEnd, time.clock() - tic))
\ No newline at end of file
Markdown is supported
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