Skip to content
GitLab
Commit 032e3c70 authored by Jonathan Lambrechts's avatar Jonathan Lambrechts
Browse files

c fluid solver with particles

parent edac09ce
Hide whitespace changes
Inline Side-by-side
fluid/fluid.pyx
View file @ 032e3c70
import numpy
import signal
from cython.view cimport array as cvarray
signal.signal(signal.SIGINT, signal.SIG_DFL)
cdef extern from "linear_system.h":
ctypedef struct StrongBoundary :
const char *tag;
int field;
double v;
cdef extern from "fluid_problem.h":
ctypedef struct FluidProblem :
double mu
double epsilon
double rho
int n_particles
double *particle_force
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);
FluidProblem *fluid_problem_new(const char *mesh_file_name, double mu, double rho, double epsilon, int n_boundary, StrongBoundary *boundaries);
void fluid_problem_free(FluidProblem *problem);
void test_numpy_c(double *v);
void fluid_problem_set_particles(FluidProblem *problem, int n, double *mass, double *volume, double *position, double *velocity);
cdef extern from "linear_system.h":
void initialize_linear_solver(int *argc, char ***argv);
void initialize_linear_solver(int argc, char **argv);
include "openmpihack.pxi"
initialize_linear_solver(NULL, NULL)
import sys
from libc.stdlib cimport malloc, free
cdef char **argv = <char**>malloc(sizeof(char*)*(len(sys.argv)));
argvbyte = [c.encode() for c in sys.argv]
for i,c in enumerate(argvbyte) :
argv[i] = argvbyte[i]
initialize_linear_solver(len(sys.argv), argv)
free(argv)
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)
def __cinit__(self, mesh_file_name, mu, rho, epsilon, strong_boundaries):
cdef StrongBoundary* bnds = <StrongBoundary*> malloc(sizeof(StrongBoundary)*len(strong_boundaries));
tags = [bnd[0].encode() for bnd in strong_boundaries]
for i, bnd in enumerate(strong_boundaries) :
bnds[i].tag = <char*>tags[i]
bnds[i].field = bnd[1]
bnds[i].v = bnd[2]
self._fp = fluid_problem_new(mesh_file_name.encode(), mu, rho, epsilon, len(strong_boundaries),bnds)
free(bnds)
if self._fp == NULL :
raise NameError("cannot create fluid problem")
def __dealloc__(self):
......@@ -46,7 +71,15 @@ cdef class fluid_problem :
fluid_problem_export(self._fp, output_dir.encode(), t, it)
def implicit_euler(self, double dt) :
fluid_problem_implicit_euler(self._fp, dt)
forces = numpy.ndarray([self._fp.n_particles,2],"d")
for i in range(self._fp.n_particles) :
forces[i,0] = self._fp.particle_force[i*2+0]
forces[i,1] = self._fp.particle_force[i*2+1]
return forces
def test_numpy(array):
cdef double[::1] a = array
test_numpy_c(&a[0])
def set_particles(self, mass, volume, position, velocity):
cdef double[:,::1] _mass = numpy.ascontiguousarray(mass)
cdef double[:,::1] _volume = numpy.ascontiguousarray(volume)
cdef double[:,::1] _position = numpy.ascontiguousarray(position)
cdef double[:,::1] _velocity = numpy.ascontiguousarray(velocity)
fluid_problem_set_particles(self._fp,mass.shape[0],&_mass[0][0],&_volume[0][0],&_position[0][0],&_velocity[0][0])
fluid/fluid_problem.c
View file @ 032e3c70
......@@ -3,19 +3,18 @@
#include <string.h>
#include <math.h>
#include "fluid_problem.h"
#include "mesh_find.h"
#define n_fields 3
#define newton_max_it 20
#define newton_atol 1e-12
#define newton_rtol 1e-12
#define newton_rtol 1e-5
#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))
......@@ -27,12 +26,109 @@ int fluid_problem_export(const FluidProblem *problem, const char *output_dir, do
return 0;
}
static void particle_drag(double u[2], double mu, double rho, double d, double c, double drag[2], double dt, double mass)
{
double normu = hypot(u[0],u[1]);
//Reynolds/|u_p-u_f|
double Re_O_u = d*c/mu*rho;
double Cd_u = 0.63*sqrt(normu)+4.8/sqrt(Re_O_u);
Cd_u = Cd_u*Cd_u;
double f = pow(c,-1.8);
double GoU = f*Cd_u*rho/2*d;
drag[0] = GoU*u[0]/(1+dt/mass*GoU);
drag[1] = GoU*u[1]/(1+dt/mass*GoU);
}
static void fluid_problem_add_particle_force(FluidProblem *problem, const double *solution, double dt) {
double *porosity = problem->porosity;
Mesh *mesh = problem->mesh;
for (int i = 0; i < problem->n_particles; ++i) {
double *xi = problem->particle_uvw + 2*i;
double phi[3] = {1-xi[0]-xi[1],xi[0],xi[1]};
int iel = problem->particle_element_id[i];
int *tri = problem->mesh->triangles+iel*3;
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 P[] = {solution[tri[0]*3+2],solution[tri[1]*3+2],solution[tri[2]*3+2]};
double vf[2] = {
phi[0]*U[0]+phi[1]*U[1]+phi[2]*U[2],
phi[0]*V[0]+phi[1]*V[1]+phi[2]*V[2]
};
double c = phi[0]*C[0]+phi[1]*C[1]+phi[2]*C[2];
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}
};
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 gradp[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]
};
double du[2] = {
problem->particle_velocity[i*2+0]-vf[0]/c,
problem->particle_velocity[i*2+1]-vf[1]/c
};
double vol = problem->particle_volume[i];
double mass = problem->particle_mass[i];
double d = 2*sqrt(vol/M_PI);
double drag[2];
particle_drag(du,problem->mu,problem->rho,d,c,drag, dt, mass);
double ddragdu[2],ddragdv[2];
double eps = 1e-8;
double duu[2] = {du[0]+eps, du[1]};
particle_drag(duu,problem->mu,problem->rho,d,c,ddragdu, dt, mass);
double duv[2] = {du[0], du[1]+eps};
particle_drag(duv,problem->mu,problem->rho,d,c,ddragdv, dt, mass);
ddragdu[0] = -(ddragdu[0]-drag[0])/(eps*c);
ddragdu[1] = -(ddragdu[1]-drag[1])/(eps*c);
ddragdv[0] = -(ddragdv[0]-drag[0])/(eps*c);
ddragdv[1] = -(ddragdv[1]-drag[1])/(eps*c);
problem->particle_force[i*2+0] = -drag[0]-vol*gradp[0];
problem->particle_force[i*2+1] = -drag[1]-9.81*(mass-problem->rho*vol)-vol*gradp[1];
double local_vector[3*n_fields];
double local_matrix[9*n_fields*n_fields];
for (int iphi = 0; iphi < 3; ++iphi) {
local_vector[iphi+0] = -drag[0]*phi[iphi];
local_vector[iphi+3] = -drag[1]*phi[iphi];
local_vector[iphi+6] = 0;
for (int jphi = 0; jphi < 3; ++jphi) {
local_matrix[(iphi+0)*9+jphi+0] = -phi[iphi]*ddragdu[0]*phi[jphi];
local_matrix[(iphi+0)*9+jphi+3] = -phi[iphi]*ddragdv[0]*phi[jphi];
local_matrix[(iphi+0)*9+jphi+6] = 0;//phi[iphi]*vol*dphi[jphi][0];
local_matrix[(iphi+3)*9+jphi+0] = -phi[iphi]*ddragdu[1]*phi[jphi];
local_matrix[(iphi+3)*9+jphi+3] = -phi[iphi]*ddragdv[1]*phi[jphi];
local_matrix[(iphi+3)*9+jphi+6] = 0;//phi[iphi]*vol*dphi[jphi][1];
local_matrix[(iphi+6)*9+jphi+0] = 0;
local_matrix[(iphi+6)*9+jphi+3] = 0;
local_matrix[(iphi+6)*9+jphi+6] = 0;
}
}
linear_system_add_to_matrix(problem->linear_system,iel,iel,local_matrix);
linear_system_add_to_rhs(problem->linear_system,iel,local_vector);
}
}
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;
const double *old_porosity = problem->old_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;
......@@ -52,16 +148,18 @@ static void fluid_problem_assemble_system(FluidProblem *problem, const double *s
{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]}
};
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;
double C[] = {porosity[tri[0]],porosity[tri[1]],porosity[tri[2]]};
double OLDC[] = {old_porosity[tri[0]],old_porosity[tri[1]],old_porosity[tri[2]]};
double DCDT[] = {(C[0]-OLDC[0])/dt,(C[1]-OLDC[1])/dt,(C[2]-OLDC[2])/dt};
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]};
......@@ -92,6 +190,7 @@ static void fluid_problem_assemble_system(FluidProblem *problem, const double *s
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];
double dcdt = phi[0]*DCDT[0]+phi[1]*DCDT[1]+phi[2]*DCDT[2];
const double jw = QW[iqp]*detj;
for (int iphi = 0; iphi < 3; ++iphi) {
double phii = phi[iphi];
......@@ -110,8 +209,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])
+phii*(du[0]+dv[1])
-phii*VEXT*dc[1]
+phii*((du[0]+dv[1])+dcdt)
);
for (int jphi = 0; jphi < 3; ++jphi) {
double phij = phi[jphi];
......@@ -138,11 +236,13 @@ static void fluid_problem_assemble_system(FluidProblem *problem, const double *s
}
free(local_vector);
free(local_matrix);
fluid_problem_add_particle_force(problem, solution, dt);
}
static int apply_boundary_conditions(const Mesh *mesh, double *solution, const StrongBoundary *bnds)
static void apply_boundary_conditions(const Mesh *mesh, double *solution, int nBnd, StrongBoundary *bnds)
{
for (const StrongBoundary *bnd = bnds; bnd->tag != NULL; bnd++) {
for (int ibnd = 0; ibnd < nBnd; ++ibnd) {
StrongBoundary *bnd = bnds + ibnd;
int iphys;
for (iphys = 0; iphys < mesh->n_phys; ++iphys) {
if (strcmp(bnd->tag, mesh->phys_name[iphys]) == 0)
......@@ -159,7 +259,7 @@ static int apply_boundary_conditions(const Mesh *mesh, double *solution, const S
void fluid_problem_implicit_euler(FluidProblem *problem, double dt)
{
const Mesh *mesh = problem->mesh;
apply_boundary_conditions(mesh, problem->solution, problem->strong_boundaries);
apply_boundary_conditions(mesh, problem->solution, problem->n_strong_boundaries, 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];
......@@ -185,7 +285,7 @@ void fluid_problem_implicit_euler(FluidProblem *problem, double dt)
free(solution_old);
}
FluidProblem *fluid_problem_new(const char *mesh_file_name, double mu, double rho, double epsilon)
FluidProblem *fluid_problem_new(const char *mesh_file_name, double mu, double rho, double epsilon, int n_strong_boundaries, StrongBoundary *strong_boundaries)
{
FluidProblem *problem = malloc(sizeof(FluidProblem));
problem->epsilon = epsilon;
......@@ -195,22 +295,34 @@ FluidProblem *fluid_problem_new(const char *mesh_file_name, double mu, double rh
if (!mesh)
return NULL;
problem->mesh = mesh;
//problem->strong_boundaries = boundaries;
problem->strong_boundaries = malloc(sizeof(StrongBoundary)*n_strong_boundaries);
problem->n_strong_boundaries = n_strong_boundaries;
for (int i = 0; i < n_strong_boundaries; ++i) {
problem->strong_boundaries[i].tag = strdup(strong_boundaries[i].tag);
problem->strong_boundaries[i].v = strong_boundaries[i].v;
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->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;
//static StrongBoundary strong_boundaries[] = {{"Bottom",1,VEXT},{"Bottom",0,0.},{"Top",0,0.},{"PtFix",2,0.},{"Box",0,0.},{NULL, 0, 0.}};
//static StrongBoundary strong_boundaries[] = {{"Top",2,0},{"Top",1,0.},{"Top",0,0.},{"Box",1,0.},{"Box",0,0.},{NULL, 0, 0.}};
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);
problem->linear_system = linear_system_new(mesh, n_fields, problem->n_strong_boundaries, problem->strong_boundaries);
problem->n_particles = 0;
problem->particle_uvw = malloc(0);
problem->particle_element_id = malloc(0);
problem->particle_position = malloc(0);
problem->particle_mass = malloc(0);
problem->particle_force = malloc(0);
problem->particle_volume = malloc(0);
problem->particle_velocity = malloc(0);
return problem;
}
......@@ -218,15 +330,77 @@ void fluid_problem_free(FluidProblem *problem)
{
free(problem->solution);
free(problem->porosity);
free(problem->old_porosity);
linear_system_free(problem->linear_system);
free(problem->particle_uvw);
free(problem->particle_element_id);
free(problem->particle_position);
free(problem->particle_mass);
free(problem->particle_force);
free(problem->particle_volume);
free(problem->particle_velocity);
for (int i = 0; i < problem->n_strong_boundaries; ++i)
free(problem->strong_boundaries[i].tag);
free(problem->strong_boundaries);
mesh_free(problem->mesh);
}
void fluid_problem_set_particles(FluidProblem *problem,
void test_numpy_c(double *v)
void fluid_problem_set_particles(FluidProblem *problem, int n, double *mass, double *volume, double *position, double *velocity)
{
printf("%g %g\n", v[0], v[1]);
v[1] += 10;
printf("%g %g\n", v[0], v[1]);
problem->n_particles = n;
free(problem->particle_uvw);
free(problem->particle_element_id);
free(problem->particle_position);
free(problem->particle_mass);
free(problem->particle_force);
free(problem->particle_volume);
free(problem->particle_velocity);
problem->particle_uvw = malloc(sizeof(double)*2*n);
problem->particle_element_id = malloc(sizeof(int)*n);
problem->particle_position = malloc(sizeof(double)*2*n);
problem->particle_mass = malloc(sizeof(double)*n);
problem->particle_force = malloc(sizeof(double)*2*n);
problem->particle_volume = malloc(sizeof(double)*n);
problem->particle_velocity = malloc(sizeof(double)*n*2);
mesh_particles_to_mesh(problem->mesh, n, position, problem->particle_element_id, problem->particle_uvw);
for (int i = 0; i < n; ++i) {
problem->particle_position[i*2+0] = position[i*2+0];
problem->particle_position[i*2+1] = position[i*2+1];
problem->particle_mass[i] = mass[i];
problem->particle_volume[i] = volume[i];
problem->particle_velocity[i*2+0] = velocity[i*2+0];
problem->particle_velocity[i*2+1] = velocity[i*2+1];
}
Mesh *mesh = problem->mesh;
double *vertex_volume = malloc(sizeof(double)*mesh->n_nodes);
for (int i = 0; i < mesh->n_nodes; ++i){
problem->old_porosity[i] = problem->porosity[i];
problem->porosity[i] = 0;
vertex_volume[i] = 0;
}
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];
vertex_volume[tri[0]] += detj/6;
vertex_volume[tri[1]] += detj/6;
vertex_volume[tri[2]] += detj/6;
}
for (int i = 0; i < n; ++i) {
const int *tri = &mesh->triangles[problem->particle_element_id[i]*3];
double u = problem->particle_uvw[i*2+0], v = problem->particle_uvw[i*2+1];
problem->porosity[tri[0]] += (1-u-v)*volume[i];
problem->porosity[tri[1]] += u*volume[i];
problem->porosity[tri[2]] += v*volume[i];
}
for (int i = 0; i < mesh->n_nodes; ++i){
problem->porosity[i] = (1-problem->porosity[i]/vertex_volume[i]);
}
free(vertex_volume);
}
fluid/fluid_problem.h
View file @ 032e3c70
......@@ -10,14 +10,25 @@ typedef struct {
double mu;
Mesh *mesh;
double *porosity;
double *old_porosity;
double *solution;
const StrongBoundary *strong_boundaries;
int n_strong_boundaries;
StrongBoundary *strong_boundaries;
LinearSystem *linear_system;
int n_particles;
double *particle_uvw;
double *particle_position;
double *particle_volume;
double *particle_velocity;
double *particle_mass;
double *particle_force; // complete force on the particle (including gravity)
int *particle_element_id;
} 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);
FluidProblem *fluid_problem_new(const char *mesh_file_name, double mu, double rho, double epsilon, 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 test_numpy_c(double *v);
#endif
fluid/linear_system.c
View file @ 032e3c70
......@@ -87,7 +87,7 @@ static void linear_system_pre_allocate_matrix(LinearSystem *lsys)
free(n_by_row);
}
static void linear_system_create_map(LinearSystem *lsys, const StrongBoundary *bnds)
static void linear_system_create_map(LinearSystem *lsys, int n_boundaries, const StrongBoundary *bnds)
{
const Mesh *mesh = lsys->mesh;
int n = mesh->n_nodes*lsys->n_fields;
......@@ -95,7 +95,8 @@ static void linear_system_create_map(LinearSystem *lsys, const StrongBoundary *b
for (int i = 0; i < n; ++i)
lsys->map[i] = 0;
int n_fixed = 0;
for (const StrongBoundary *bnd = bnds; bnd->tag != NULL; bnd++) {
for (int ibnd = 0; ibnd < n_boundaries; ++ibnd) {
const StrongBoundary *bnd = bnds + ibnd;
int iphys;
for (iphys = 0; iphys < mesh->n_phys; ++iphys) {
if (strcmp(bnd->tag, mesh->phys_name[iphys]) == 0)
......@@ -113,12 +114,12 @@ static void linear_system_create_map(LinearSystem *lsys, const StrongBoundary *b
lsys->fixed = malloc(sizeof(double)*n_fixed);
}
LinearSystem *linear_system_new(Mesh *mesh, int n_fields, const StrongBoundary *boundaries)
LinearSystem *linear_system_new(Mesh *mesh, int n_fields, int n_boundaries, const StrongBoundary *boundaries)
{
LinearSystem *lsys = malloc(sizeof(LinearSystem));
lsys->mesh = mesh;
lsys->n_fields = n_fields;
linear_system_create_map(lsys, boundaries);
linear_system_create_map(lsys, n_boundaries, boundaries);
MatCreate(MPI_COMM_WORLD, &lsys->a);
MatSetSizes(lsys->a, lsys->n_dofs, lsys->n_dofs, PETSC_DETERMINE, PETSC_DETERMINE);
MatSetFromOptions(lsys->a);
......@@ -219,6 +220,7 @@ double linear_system_get_rhs_norm(LinearSystem *lsys)
return norm;
}
void initialize_linear_solver(int *argc, char ***argv) {
PetscInitialize(argc, argv, NULL, NULL);
void initialize_linear_solver(int argc, char **argv) {
MPI_Init(NULL, NULL);
PetscInitialize(&argc, &argv, NULL, NULL);
}
fluid/linear_system.h
View file @ 032e3c70
......@@ -6,18 +6,18 @@
typedef struct LinearSystemStruct LinearSystem;
typedef struct {
const char *tag;
char *tag;
int field;
double v;
}StrongBoundary;
LinearSystem *linear_system_new(Mesh *mesh, int n_fields, const StrongBoundary *boundaries);
LinearSystem *linear_system_new(Mesh *mesh, int n_fields, int n_boundaries, 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);
void initialize_linear_solver(int argc, char **argv);
#endif
fluid/mesh.c
View file @ 032e3c70
......@@ -6,11 +6,11 @@
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");
printf("Invalid mesh file (\"%s\" expected).\n", w);
return -1;
}
if(strcmp(word, w) != 0){
printf("Invalid mesh file (\"%s\" expected).\n");
printf("Invalid mesh file (\"%s\" expected).\n", w);
return -1;
}
return 0;
......@@ -131,7 +131,7 @@ Mesh *mesh_load_msh(const char *filename)
m->n_triangles = 0;
for (int i = 0; i < n_el; ++i) {
int eid, etype, nflags, flag;
int phys;
int phys=0;
if (fscanf(f, "%d %d %d", &eid, &etype, &nflags) != 3){
printf("Cannot read elements.\n");
return NULL;
......
fluid/mesh_find.c 0 → 100644
View file @ 032e3c70
#include <math.h>
#include <stdlib.h>
#include "quadtree.h"
#include "mesh.h"
#define DIMENSION 2
#include "vector.h"
static void _bboxtol(double *bbmin, double *bbmax) {
double lmax = 0;
for (int i = 0; i < DIMENSION; ++i) {
lmax = fmax(lmax, bbmax[i] - bbmin[i]);
}
double tol = 1e-8 * lmax;
for (int i = 0; i < DIMENSION; ++i) {
bbmax[i] += tol;
bbmin[i] -= tol;
}
}
double volTet(double x0[3], double x1[3], double x2[3], double x3[3]) {
// x1 * (y2 * (z3 - z4) - z2 * (y3 - y4) + y3 * z4 - y4 * z3)
//- y1 * (x2 * (z3 - z4) - z2 * (x3 - x4) + x3 * z4 - x4 * z3)
//+ z1 * (x2 * (y3 - y4) - y2 * (x3 - x4) + x3 * y4 - x4 * y3)
//- (x2 * (y3 * z4 - y4 * z3) - y2 * (x3 * z4 - x4 * z3) + z2 * (x3 * y4 - x4 * y3))
return x0[0] * (x1[1] * (x2[2] - x3[2]) - x1[2] * (x2[1] - x3[1]) + x2[1] * x3[2] - x3[1] * x2[2])
- x0[1] * (x1[0] * (x2[2] - x3[2]) - x1[2] * (x2[0] - x3[0]) + x2[0] * x3[2] - x3[0] * x2[2])
+ x0[2] * (x1[0] * (x2[1] - x3[1]) - x1[1] * (x2[0] - x3[0]) + x2[0] * x3[1] - x3[0] * x2[1])
- (x1[0] * (x2[1] * x3[2] - x3[1] * x2[2]) - x1[1] * (x2[0] * x3[2] - x3[0] * x2[2]) + x1[2] * (x2[0] * x3[1] - x3[0] * x2[1]));
}
void particle_to_mesh(size_t nParticles, double *particles, int nElements, double *elements, int *elid, double *Xi)
{
double bbmin[DIMENSION], bbmax[DIMENSION];
int N = DIMENSION + 1;
for (int i = 0; i < DIMENSION; ++i) {
bbmin[i] = elements[i];