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

fluid solver with cython interface

parent f18ae7c6
import numpy
from cython.view cimport array as cvarray
cdef extern from "fluid_problem.h":
ctypedef struct FluidProblem :
double mu
double epsilon
double rho
pass
int fluid_problem_export(const FluidProblem *problem, const char *output_dir, double t, double it);
void fluid_problem_implicit_euler(FluidProblem *problem, double dt);
FluidProblem *fluid_problem_new(const char *mesh_file_name, double mu, double rho, double epsilon);
void fluid_problem_free(FluidProblem *problem);
void test_numpy_c(double *v);
cdef extern from "linear_system.h":
void initialize_linear_solver(int *argc, char ***argv);
include "openmpihack.pxi"
initialize_linear_solver(NULL, NULL)
cdef class fluid_problem :
cdef FluidProblem *_fp
def __cinit__(self, mesh_file_name, mu, rho, epsilon):
self._fp = fluid_problem_new(mesh_file_name.encode(), mu, rho, epsilon)
if self._fp == NULL :
raise NameError("cannot create fluid problem")
def __dealloc__(self):
if(self._fp != NULL) :
fluid_problem_free(self._fp)
def __getattr__(self, name) :
if name == "mu" :
return self._fp.mu
elif name == "rho" :
return self._fp.rho
elif name == "epsilon" :
return self._fp.epsilon
def __setattr__(self, name, val) :
if name == "mu" :
self._fp.mu = val
if name == "rho" :
self._fp.rho = val
elif name == "epsilon" :
self._fp.epsilon = val
def export(self, output_dir, double t, int it) :
fluid_problem_export(self._fp, output_dir.encode(), t, it)
def implicit_euler(self, double dt) :
fluid_problem_implicit_euler(self._fp, dt)
def test_numpy(array):
cdef double[::1] a = array
test_numpy_c(&a[0])
#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <math.h>
#include "fluid_problem.h"
#define n_fields 3
#define newton_max_it 20
#define newton_atol 1e-12
#define newton_rtol 1e-12
#define N_QP 3
static double QP[][2] = {{1./6,1./6},{1./6,2./3},{2./3,1./6}};
static double QW[] = {1./6,1./6,1./6};
#define VEXT 0.01
int fluid_problem_export(const FluidProblem *problem, const char *output_dir, double t, double iter)
{
if(mesh_write_msh_scalar(problem->mesh, output_dir, "porosity", t, iter, problem->porosity, 1, 0))
return -1;
if(mesh_write_msh_vector(problem->mesh, output_dir, "uv", t, iter, problem->solution, n_fields, 0, 1))
return -1;
if(mesh_write_msh_scalar(problem->mesh, output_dir, "p", t, iter, problem->solution, n_fields, 2))
return -1;
return 0;
}
static void fluid_problem_assemble_system(FluidProblem *problem, const double *solution_old, double dt)
{
LinearSystem *lsys = problem->linear_system;
const Mesh *mesh = problem->mesh;
const double *solution = problem->solution;
const double *porosity = problem->porosity;
double *local_vector = malloc(sizeof(double)*n_fields*3);
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;
for (int iel=0; iel < mesh->n_triangles; ++iel) {
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]}};
const double dxdxi[2][2] = {
{x[1][0]-x[0][0],x[2][0]-x[0][0]},
{x[1][1]-x[0][1],x[2][1]-x[0][1]}};
const double detj = dxdxi[0][0]*dxdxi[1][1]-dxdxi[0][1]*dxdxi[1][0];
const double dxidx[2][2] = {
{dxdxi[1][1]/detj,-dxdxi[0][1]/detj},
{-dxdxi[1][0]/detj,dxdxi[0][0]/detj}
};
for (int i=0; i < 3*n_fields; ++i)
local_vector[i] = 0;
for (int i=0; i < 9*n_fields*n_fields; ++i)
local_matrix[i] = 0;
const double dphi[][2] = {
{-dxidx[0][0]-dxidx[1][0],-dxidx[0][1]-dxidx[1][1]},
{dxidx[0][0],dxidx[0][1]},
{dxidx[1][0],dxidx[1][1]}
};
double C[] = {porosity[tri[0]],porosity[tri[1]],porosity[tri[2]]};
double U[] = {solution[tri[0]*3+0],solution[tri[1]*3+0],solution[tri[2]*3+0]};
double V[] = {solution[tri[0]*3+1],solution[tri[1]*3+1],solution[tri[2]*3+1]};
double U_old[] = {solution_old[tri[0]*3+0],solution_old[tri[1]*3+0],solution_old[tri[2]*3+0]};
double V_old[] = {solution_old[tri[0]*3+1],solution_old[tri[1]*3+1],solution_old[tri[2]*3+1]};
double P[] = {solution[tri[0]*3+2],solution[tri[1]*3+2],solution[tri[2]*3+2]};
double dc[2] = {
dphi[0][0]*C[0]+dphi[1][0]*C[1]+dphi[2][0]*C[2],
dphi[0][1]*C[0]+dphi[1][1]*C[1]+dphi[2][1]*C[2]
};
double du[2] = {
dphi[0][0]*U[0]+dphi[1][0]*U[1]+dphi[2][0]*U[2],
dphi[0][1]*U[0]+dphi[1][1]*U[1]+dphi[2][1]*U[2]
};
double dv[2] = {
dphi[0][0]*V[0]+dphi[1][0]*V[1]+dphi[2][0]*V[2],
dphi[0][1]*V[0]+dphi[1][1]*V[1]+dphi[2][1]*V[2]
};
double dp[2] = {
dphi[0][0]*P[0]+dphi[1][0]*P[1]+dphi[2][0]*P[2],
dphi[0][1]*P[0]+dphi[1][1]*P[1]+dphi[2][1]*P[2]
};
for (int iqp = 0; iqp< N_QP; ++iqp) {
const double xi = QP[iqp][0], eta = QP[iqp][1];
const double phi[] = {1-xi-eta, xi, eta};
double u = phi[0]*U[0]+phi[1]*U[1]+phi[2]*U[2];
double v = phi[0]*V[0]+phi[1]*V[1]+phi[2]*V[2];
double u_old = phi[0]*U_old[0]+phi[1]*U_old[1]+phi[2]*U_old[2];
double v_old = phi[0]*V_old[0]+phi[1]*V_old[1]+phi[2]*V_old[2];
double p = phi[0]*P[0]+phi[1]*P[1]+phi[2]*P[2];
double c = phi[0]*C[0]+phi[1]*C[1]+phi[2]*C[2];
const double jw = QW[iqp]*detj;
for (int iphi = 0; iphi < 3; ++iphi) {
double phii = phi[iphi];
double dphii[2] = {dphi[iphi][0],dphi[iphi][1]};
double tau[2][2] = {{du[0]-u*dc[0]/c,du[1]-u*dc[1]/c},
{dv[0]-v*dc[0]/c,dv[1]-v*dc[1]/c}};
local_vector[iphi+0] += jw*(
mu*(dphii[0]*tau[0][0]+dphii[1]*(tau[0][1]+tau[1][0])*0.5)
+rho*phii*((u-u_old)/dt+(u*tau[0][0]+v*tau[0][1])/c)
-p*dphii[0]
);
local_vector[iphi+3] += jw*(
mu*(dphii[0]*(tau[0][1]+tau[1][0])*0.5+dphii[1]*tau[1][1])
+rho*phii*((v-v_old)/dt+(u*tau[1][0]+v*tau[1][1])/c)
-p*dphii[1]
);
local_vector[iphi+6] += jw*(
epsilon*(dphii[0]*dp[0]+dphii[1]*dp[1])
+phii*(du[0]+dv[1])
-phii*VEXT*dc[1]
);
for (int jphi = 0; jphi < 3; ++jphi) {
double phij = phi[jphi];
double dphij[2] = {dphi[jphi][0],dphi[jphi][1]};
double dtau[2] = {dphij[0]-phij*dc[0]/c,dphij[1]-phij*dc[1]/c};
local_matrix[(iphi+0)*9+jphi+0] += jw*mu*(dphii[0]*dtau[0]+dphii[1]*dtau[1]*0.5) // U U
+jw*rho*phii*(phij/dt+(phij*tau[0][0]+u*dtau[0]+v*dtau[1])/c);
local_matrix[(iphi+0)*9+jphi+3] += jw*mu*dphii[1]*dtau[0]*0.5 // U V
+jw*rho*phii*phij*tau[0][1]/c;
local_matrix[(iphi+0)*9+jphi+6] +=-jw*dphii[0]*phij; // U P
local_matrix[(iphi+3)*9+jphi+0] += jw*mu*dphii[0]*dtau[1]*0.5 // V U
+jw*rho*phii*phij*tau[1][0]/c;
local_matrix[(iphi+3)*9+jphi+3] += jw*mu*(dphii[0]*dtau[0]*0.5+dphii[1]*dtau[1]) // V V
+jw*rho*phii*(phij/dt+(phij*tau[1][1]+u*dtau[0]+v*dtau[1])/c);
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
}
}
}
linear_system_add_to_matrix(lsys,iel,iel,local_matrix);
linear_system_add_to_rhs(lsys,iel,local_vector);
}
free(local_vector);
free(local_matrix);
}
static int apply_boundary_conditions(const Mesh *mesh, double *solution, const StrongBoundary *bnds)
{
for (const StrongBoundary *bnd = bnds; bnd->tag != NULL; bnd++) {
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)
solution[mesh->phys_nodes[iphys][i]*n_fields+bnd->field] = bnd->v;
}
}
void fluid_problem_implicit_euler(FluidProblem *problem, double dt)
{
const Mesh *mesh = problem->mesh;
apply_boundary_conditions(mesh, problem->solution, problem->strong_boundaries);
double *solution_old = malloc(mesh->n_nodes*n_fields*sizeof(double));
for (int i=0; i<mesh->n_nodes*n_fields; ++i)
solution_old[i] = problem->solution[i];
double *dx = malloc(sizeof(double)*mesh->n_nodes*n_fields);
double norm0 = 0;
for (int i=0; i<newton_max_it; ++i) {
linear_system_zero_matrix_and_rhs(problem->linear_system);
fluid_problem_assemble_system(problem, solution_old, dt);
double norm = linear_system_get_rhs_norm(problem->linear_system);
printf("iter %i norm = %g\n", i, norm);
if (norm < newton_atol)
break;
if (i == 0)
norm0 = norm;
else if(norm/norm0 < newton_rtol)
break;
linear_system_solve(problem->linear_system, dx);
for (int j=0; j<mesh->n_nodes*n_fields; ++j) {
problem->solution[j] -= dx[j];
}
}
free(dx);
free(solution_old);
}
FluidProblem *fluid_problem_new(const char *mesh_file_name, double mu, double rho, double epsilon)
{
FluidProblem *problem = malloc(sizeof(FluidProblem));
problem->epsilon = epsilon;
problem->rho = rho;
problem->mu = mu;
Mesh *mesh = mesh_load_msh(mesh_file_name);
if (!mesh)
return NULL;
problem->mesh = mesh;
//problem->strong_boundaries = boundaries;
problem->porosity = malloc(mesh->n_nodes*sizeof(double));
problem->solution = malloc(mesh->n_nodes*n_fields*sizeof(double));
// begin to remove
static StrongBoundary strong_boundaries[] = {{"Bottom",1,VEXT},{"Bottom",0,0.},{"Top",0,0.},{"PtFix",2,0.},{"Box",0,0.},{NULL, 0, 0.}};
problem->strong_boundaries = strong_boundaries;
for (int i = 0; i < problem->mesh->n_nodes; ++i){
double x = problem->mesh->x[i*3+0];
double y = problem->mesh->x[i*3+1];
problem->porosity[i] = 1-0.99*exp(-10*(x*x+y*y));
problem->solution[i*3+0] = 0.;
problem->solution[i*3+1] = 0;
problem->solution[i*3+2] = 0.;
}
// end to remove
problem->linear_system = linear_system_new(mesh, n_fields, problem->strong_boundaries);
return problem;
}
void fluid_problem_free(FluidProblem *problem)
{
free(problem->solution);
free(problem->porosity);
linear_system_free(problem->linear_system);
mesh_free(problem->mesh);
}
void fluid_problem_set_particles(FluidProblem *problem,
void test_numpy_c(double *v)
{
printf("%g %g\n", v[0], v[1]);
v[1] += 10;
printf("%g %g\n", v[0], v[1]);
}
#ifndef FLUID_PROBLEM_H
#define FLUID_PROBLEM_H
#include "mesh.h"
#include "linear_system.h"
typedef struct {
double epsilon;
double rho;
double mu;
Mesh *mesh;
double *porosity;
double *solution;
const StrongBoundary *strong_boundaries;
LinearSystem *linear_system;
} FluidProblem;
int fluid_problem_export(const FluidProblem *problem, const char *output_dir, double t, double iter);
void fluid_problem_implicit_euler(FluidProblem *problem, double dt);
FluidProblem *fluid_problem_new(const char *mesh_file_name, double mu, double rho, double epsilon);
void fluid_problem_free(FluidProblem *problem);
void test_numpy_c(double *v);
#endif
#include "linear_system.h"
#include "petscksp.h"
#include "petscvec.h"
struct LinearSystemStruct {
Mesh *mesh;
int n_fields;
int *map; //[inode*n_fields+ifield] >=0 -> dof_id, <0 -> fixed_id-1
int n_dofs;
double *fixed;
KSP ksp;
Mat a;
Vec b,x;
int assembly_needed;
};
#define CSIZE 10
typedef struct SparsityChunkStruct {
int n[CSIZE];
struct SparsityChunkStruct *next;
}SparsityChunk;
SparsityChunk *sparsity_chunk_new(void)
{
SparsityChunk *s = malloc(sizeof(SparsityChunk));
for (int i = 0; i < CSIZE; ++i)
s->n[i] = -1;
s->next = NULL;
return s;
}
void sparsity_chunk_insert(SparsityChunk *s, int id)
{
for (int i=0; i<CSIZE; ++i) {
if (s->n[i] == -1) s->n[i] = id;
if (s->n[i] == id) return;
}
if (s->next == NULL)
s->next = sparsity_chunk_new();
sparsity_chunk_insert(s->next, id);
}
void sparsity_chunk_free(SparsityChunk *s)
{
if (s->next)
sparsity_chunk_free(s->next);
free(s);
}
static void linear_system_pre_allocate_matrix(LinearSystem *lsys)
{
const Mesh *mesh = lsys->mesh;
SparsityChunk **sparsity = malloc(sizeof(SparsityChunk*)*lsys->n_dofs);
for (int i = 0; i < lsys->n_dofs; ++i)
sparsity[i] = sparsity_chunk_new();
for (int i = 0; i < mesh->n_triangles; ++i) {
for (int k = 0; k < 3; ++k) {
for (int kf = 0; kf < lsys->n_fields; ++kf) {
int row = lsys->map[mesh->triangles[i*3+k]*lsys->n_fields+kf];
if (row < 0)
continue;
for (int l = 0; l < 3; ++l){
for (int lf = 0; lf < lsys->n_fields; ++lf){
int col = lsys->map[mesh->triangles[i*3+l]*lsys->n_fields+lf];
if (col >= 0)
sparsity_chunk_insert(sparsity[row], col);
}
}
}
}
}
int *n_by_row = malloc(sizeof(int) * lsys->n_dofs);
for (int i = 0; i < lsys->n_dofs; ++i){
n_by_row[i] = 0;
SparsityChunk *c = sparsity[i];
while (c->next) {
n_by_row[i] += CSIZE;
c = c->next;
}
for (int j = 0; j < CSIZE; ++j)
if(c->n[j] != -1)
n_by_row[i]++;
sparsity_chunk_free(sparsity[i]);
}
free(sparsity);
MatSeqAIJSetPreallocation(lsys->a, 0, n_by_row);
free(n_by_row);
}
static void linear_system_create_map(LinearSystem *lsys, const StrongBoundary *bnds)
{
const Mesh *mesh = lsys->mesh;
int n = mesh->n_nodes*lsys->n_fields;
lsys->map = malloc(sizeof(int)*n);
for (int i = 0; i < n; ++i)
lsys->map[i] = 0;
int n_fixed = 0;
for (const StrongBoundary *bnd = bnds; bnd->tag != NULL; bnd++) {
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)
lsys->map[mesh->phys_nodes[iphys][i]*lsys->n_fields+bnd->field] = -(++n_fixed);
}
lsys->n_dofs = 0;
for (int i = 0; i < n; ++i)
if(lsys->map[i]==0)
lsys->map[i] = lsys->n_dofs++;
lsys->fixed = malloc(sizeof(double)*n_fixed);
}
LinearSystem *linear_system_new(Mesh *mesh, int n_fields, const StrongBoundary *boundaries)
{
LinearSystem *lsys = malloc(sizeof(LinearSystem));
lsys->mesh = mesh;
lsys->n_fields = n_fields;
linear_system_create_map(lsys, boundaries);
MatCreate(MPI_COMM_WORLD, &lsys->a);
MatSetSizes(lsys->a, lsys->n_dofs, lsys->n_dofs, PETSC_DETERMINE, PETSC_DETERMINE);
MatSetFromOptions(lsys->a);
linear_system_pre_allocate_matrix(lsys);
VecCreate(MPI_COMM_WORLD, &lsys->b);
VecSetSizes(lsys->b, lsys->n_dofs, PETSC_DETERMINE);
VecSetFromOptions(lsys->b);
VecSetOption(lsys->b, VEC_IGNORE_NEGATIVE_INDICES,PETSC_TRUE);
VecDuplicate(lsys->b, &lsys->x);
KSPCreate(MPI_COMM_WORLD, &lsys->ksp);
KSPSetFromOptions(lsys->ksp);
lsys->assembly_needed = 0;
return lsys;
}
void linear_system_free(LinearSystem *lsys)
{
KSPDestroy(&lsys->ksp);
VecDestroy(&lsys->b);
VecDestroy(&lsys->x);
MatDestroy(&lsys->a);
free(lsys->fixed);
free(lsys->map);
free(lsys);
}
void linear_system_add_to_matrix(LinearSystem *lsys, int el0, int el1, const double *local_matrix)
{
int *tri0 = &lsys->mesh->triangles[el0*3];
int *tri1 = &lsys->mesh->triangles[el1*3];
int nf = lsys->n_fields;
int *map0 = malloc(sizeof(int)*3*nf);
int *map1 = malloc(sizeof(int)*3*nf);
for(int i=0; i<3; ++i){
for (int j=0; j<nf; ++j){
map0[j*3+i] = lsys->map[tri0[i]*nf+j];
map1[j*3+i] = lsys->map[tri1[i]*nf+j];
}
}
MatSetValues(lsys->a,3*nf,map0,3*nf,map1,local_matrix,ADD_VALUES);
free(map0);
free(map1);
lsys->assembly_needed = 1;
}
void linear_system_add_to_rhs(LinearSystem *lsys, int el0, const double *local_vector)
{
int *tri0 = &lsys->mesh->triangles[el0*3];
int nf = lsys->n_fields;
int *map0 = malloc(sizeof(int)*3*nf);
for(int i=0; i<3; ++i){
for (int j=0; j<nf; ++j){
map0[j*3+i]= lsys->map[tri0[i]*nf+j];
}
}
VecSetValues(lsys->b,3*nf,map0,local_vector,ADD_VALUES);
free(map0);
}
void linear_system_zero_matrix_and_rhs(LinearSystem *lsys)
{
if(lsys->assembly_needed) {
MatAssemblyBegin(lsys->a, MAT_FINAL_ASSEMBLY);
MatAssemblyEnd(lsys->a, MAT_FINAL_ASSEMBLY);
lsys->assembly_needed = 0;
}
MatZeroEntries(lsys->a);
VecZeroEntries(lsys->b);
}
void linear_system_solve(LinearSystem *lsys, double *solution)
{
VecAssemblyBegin(lsys->b);
VecAssemblyEnd(lsys->b);
if(lsys->assembly_needed) {
MatAssemblyBegin(lsys->a, MAT_FINAL_ASSEMBLY);
MatAssemblyEnd(lsys->a, MAT_FINAL_ASSEMBLY);
lsys->assembly_needed = 0;
}
KSPSetOperators(lsys->ksp, lsys->a, lsys->a);
KSPSolve(lsys->ksp, lsys->b, lsys->x);
double *x;
VecGetArray(lsys->x, &x);
int n_row = lsys->n_fields*lsys->mesh->n_nodes;
for (int i=0; i<n_row; ++i) {
if (lsys->map[i] >= 0)
solution[i] = x[lsys->map[i]];
else
solution[i] = 0;
}
VecRestoreArray(lsys->x,&x);
}
double linear_system_get_rhs_norm(LinearSystem *lsys)
{
double norm;
VecNorm(lsys->b, NORM_2, &norm);
return norm;
}
void initialize_linear_solver(int *argc, char ***argv) {
PetscInitialize(argc, argv, NULL, NULL);
}
#ifndef LINEAR_SYSTEM_H
#define LINEAR_SYSTEM_H
#include "mesh.h"
typedef struct LinearSystemStruct LinearSystem;
typedef struct {
const char *tag;
int field;
double v;
}StrongBoundary;
LinearSystem *linear_system_new(Mesh *mesh, int n_fields, const StrongBoundary *boundaries);
void linear_system_add_to_matrix(LinearSystem *lsys, int el0, int el1, const double *local_matrix);
void linear_system_add_to_rhs(LinearSystem *lsys, int el0, const double *local_vector);
void linear_system_zero_matrix_and_rhs(LinearSystem *lsys);
void linear_system_solve(LinearSystem *lsys, double *solution);
void linear_system_free(LinearSystem *lsys);
double linear_system_get_rhs_norm(LinearSystem *lsys);
void initialize_linear_solver(int *argc, char ***argv);
#endif
#include "mesh.h"
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
static int check_word_in_file(FILE *f, const char *w) {
char word[256];
if(fscanf(f,"%255s", word) != 1){
printf("Invalid mesh file (\"%s\" expected).\n");
return -1;
}
if(strcmp(word, w) != 0){
printf("Invalid mesh file (\"%s\" expected).\n");
return -1;
}