Commit 0c732bd2 authored by Jonathan Lambrechts's avatar Jonathan Lambrechts
Browse files

home-made linear system

parent 9d7f2d23
all :
python2 setup.py build_ext --inplace
python setup.py build_ext --inplace
clean :
python setup.py clean
......@@ -18,8 +18,9 @@ cdef extern from "fluid_problem.h":
int n_particles
double *particle_force
pass
int fluid_problem_export(const FluidProblem *problem, const char *output_dir, double t, double it);
int fluid_problem_export(const FluidProblem *problem, const char *output_dir, double t, int it);
void fluid_problem_implicit_euler(FluidProblem *problem, double dt);
void fluid_problem_adapt_mesh(FluidProblem *problem, double gradmin, double gradmax, double lcmin);
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 fluid_problem_set_particles(FluidProblem *problem, int n, double *mass, double *volume, double *position, double *velocity);
......@@ -70,6 +71,8 @@ cdef class fluid_problem :
self._fp.rho = val
elif name == "epsilon" :
self._fp.epsilon = val
def adapt_mesh(self, double gradmin, double gradmax, double lcmin) :
fluid_problem_adapt_mesh(self._fp, gradmin, gradmax, lcmin)
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) :
......
#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <float.h>
#include <math.h>
#include "fluid_problem.h"
#include "mesh_find.h"
......@@ -15,14 +16,24 @@
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};
int fluid_problem_export(const FluidProblem *problem, const char *output_dir, double t, double iter)
int fluid_problem_export(const FluidProblem *problem, const char *output_dir, double t, int iter)
{
if(mesh_write_msh_scalar(problem->mesh, output_dir, "porosity", t, iter, problem->porosity, 1, 0))
char file_name[1024];
sprintf(file_name,"%s/fluid_%05i.msh",output_dir, iter);
FILE *f = fopen(file_name, "w");
if (!f){
printf("Cannot open file \"%s\" for writing.\n", file_name);
return -1;
if(mesh_write_msh_vector(problem->mesh, output_dir, "uv", t, iter, problem->solution, n_fields, 0, 1))
}
if(mesh_write_msh(problem->mesh, f))
return -1;
if(mesh_write_msh_scalar(problem->mesh, f, "porosity", t, iter, problem->porosity, 1, 0))
return -1;
if(mesh_write_msh_scalar(problem->mesh, output_dir, "p", t, iter, problem->solution, n_fields, 2))
if(mesh_write_msh_vector(problem->mesh, f, "uv", t, iter, problem->solution, n_fields, 0, 1))
return -1;
if(mesh_write_msh_scalar(problem->mesh, f, "p", t, iter, problem->solution, n_fields, 2))
return -1;
fclose(f);
return 0;
}
......@@ -345,38 +356,12 @@ void fluid_problem_free(FluidProblem *problem)
mesh_free(problem->mesh);
}
void fluid_problem_set_particles(FluidProblem *problem, int n, double *mass, double *volume, double *position, double *velocity)
static void fluid_problem_compute_porosity(FluidProblem *problem)
{
if(problem->n_particles != n) {
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);
double *volume = problem->particle_volume;
for (int i = 0; i < mesh->n_nodes; ++i){
problem->old_porosity[i] = problem->porosity[i];
problem->porosity[i] = 0;
vertex_volume[i] = 0;
}
......@@ -394,7 +379,7 @@ void fluid_problem_set_particles(FluidProblem *problem, int n, double *mass, dou
vertex_volume[tri[1]] += detj/6;
vertex_volume[tri[2]] += detj/6;
}
for (int i = 0; i < n; ++i) {
for (int i = 0; i < problem->n_particles; ++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];
......@@ -406,3 +391,135 @@ void fluid_problem_set_particles(FluidProblem *problem, int n, double *mass, dou
}
free(vertex_volume);
}
void fluid_problem_adapt_mesh(FluidProblem *problem, double gradmin, double gradmax, double lcmin)
{
Mesh *mesh = problem->mesh;
double *solution = problem->solution;
double *new_mesh_size = malloc(sizeof(double)*mesh->n_nodes);
for (int i = 0; i < mesh->n_nodes; ++i)
new_mesh_size[i] = DBL_MAX;
gradmax = sqrt(gradmax);
gradmin = sqrt(gradmin);
for (int iel = 0; iel < mesh->n_triangles; ++iel) {
int *tri = problem->mesh->triangles+iel*3;
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]};
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 gradU[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 gradV[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 ngrad = pow(gradU[0]*gradU[0]+gradU[1]*gradU[1] + gradV[0]*gradV[0] + gradV[1]*gradV[1], 0.25);
double lc = lcmin /fmin(1,fmax(ngrad/gradmax,gradmin/gradmax));
for (int j = 0; j < 3; ++j)
new_mesh_size[tri[j]] = fmin(new_mesh_size[tri[j]], lc);
}
FILE *f = fopen("adapt/lc.pos", "w");
fprintf(f, "View \"lc\" {\n");
for (int iel = 0; iel < mesh->n_triangles; ++iel) {
int *tri = problem->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]}};
double LC[3] = {new_mesh_size[tri[0]], new_mesh_size[tri[1]], new_mesh_size[tri[2]]};
fprintf(f, "ST(%.8g, %.8g, 0, %.8g, %.8g, 0, %.8g, %.8g, 0){%.8g, %.8g, %.8g};\n",
x[0][0],x[0][1],x[1][0],x[1][1],x[2][0],x[2][1],
LC[0],LC[1],LC[2]);
}
fprintf(f, "};\n");
fclose(f);
free(new_mesh_size);
system("gmsh -2 adapt/mesh.geo");
Mesh *new_mesh = mesh_load_msh("adapt/mesh.msh");
double *new_solution = malloc(sizeof(double)*new_mesh->n_nodes*3);
double *new_xi = malloc(sizeof(double)*new_mesh->n_nodes*2);
int *new_eid = malloc(sizeof(int)*new_mesh->n_nodes);
double *newx = malloc(sizeof(double)*2*new_mesh->n_nodes);
for (int i = 0;i < new_mesh->n_nodes;++i) {
newx[i*2+0] = new_mesh->x[i*3+0];
newx[i*2+1] = new_mesh->x[i*3+1];
}
mesh_particles_to_mesh(mesh, new_mesh->n_nodes, newx, new_eid, new_xi);
for (int i = 0; i < new_mesh->n_nodes; ++i) {
int *tri = problem->mesh->triangles+new_eid[i]*3;
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 *xi = new_xi+i*2;
double phi[] = {1-xi[0]-xi[1],xi[0],xi[1]};
new_solution[i*3+0] = U[0]*phi[0]+U[1]*phi[1]+U[2]*phi[2];
new_solution[i*3+1] = V[0]*phi[0]+V[1]*phi[1]+V[2]*phi[2];
new_solution[i*3+2] = 0.;
}
free(new_eid);
free(new_xi);
free(newx);
free(problem->solution);
problem->solution = new_solution;
free(problem->porosity);
problem->porosity = malloc(sizeof(double)*new_mesh->n_nodes);
free(problem->old_porosity);
problem->old_porosity = malloc(sizeof(double)*new_mesh->n_nodes);
mesh_free(problem->mesh);
problem->mesh = new_mesh;
mesh_particles_to_mesh(problem->mesh, problem->n_particles, problem->particle_position, problem->particle_element_id, problem->particle_uvw);
fluid_problem_compute_porosity(problem);
linear_system_free(problem->linear_system);
problem->linear_system = linear_system_new(problem->mesh, n_fields, problem->n_strong_boundaries, problem->strong_boundaries);
}
void fluid_problem_set_particles(FluidProblem *problem, int n, double *mass, double *volume, double *position, double *velocity)
{
if(problem->n_particles != n) {
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];
}
for (int i = 0; i < problem->mesh->n_nodes; ++i){
problem->old_porosity[i] = problem->porosity[i];
}
fluid_problem_compute_porosity(problem);
}
......@@ -26,9 +26,10 @@ typedef struct {
int *particle_element_id;
} FluidProblem;
int fluid_problem_export(const FluidProblem *problem, const char *output_dir, double t, double iter);
int fluid_problem_export(const FluidProblem *problem, const char *output_dir, double t, int 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, 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);
#endif
#include <math.h>
#include "mesh.h"
#include "linear_system.h"
#include <stdlib.h>
#include <string.h>
#define CONMAX 12
void connectivity_insert(int *connectivity, int i, int j)
{
for (int k = 0; k < CONMAX; ++k) {
int *p = connectivity + CONMAX*i + k;
if (*p == -1)
*p = j;
if (*p == j)
return;
}
printf("ERROR : node %i has more than %i neighbours\n", i, CONMAX);
}
int reverse_cuthill_mckee(Mesh *mesh, int *ordering)
{
int *node_connectivity = malloc(sizeof(int)*mesh->n_nodes*CONMAX);
for (int i = 0; i < mesh->n_nodes*CONMAX; ++i) {
node_connectivity[i] = -1;
}
for (int i = 0; i < mesh->n_triangles; ++i) {
int *tri = mesh->triangles + i*3;
connectivity_insert(node_connectivity, tri[0], tri[1]);
connectivity_insert(node_connectivity, tri[1], tri[0]);
connectivity_insert(node_connectivity, tri[0], tri[2]);
connectivity_insert(node_connectivity, tri[2], tri[0]);
connectivity_insert(node_connectivity, tri[1], tri[2]);
connectivity_insert(node_connectivity, tri[2], tri[1]);
}
int *node_degree = malloc(sizeof(int)*mesh->n_nodes);
for (int i = 0; i < mesh->n_nodes; ++i) {
node_degree[i] = 0;
for (int j = 0; j < CONMAX; ++j) {
if (node_connectivity[CONMAX*i+j] == -1)
break;
node_degree[i] += 1;
}
}
int *queue = malloc(sizeof(int)*mesh->n_nodes);
queue[0] = 0;
for (int i = 0; i < mesh->n_nodes; ++i){
ordering[i] = -1;
if (node_degree[queue[0]] > node_degree[i] )
queue[0] = i;
}
int stage_start = 0;
int stage_end = 1;
int queue_end = 1;
int id = 0;
while(stage_start != stage_end) {
for (int i = stage_start; i < stage_end; ++i) {
int c = queue[i];
ordering[c] = mesh->n_nodes-1 -(id++);
for(int j = 0; j < node_degree[c]; ++j) {
int o = node_connectivity[c*CONMAX+j];
if (ordering[o] == -1) {
ordering[o] = -2;
#if 1
queue[queue_end++] = o;
#else
int k = stage_end;
while(k < queue_end && node_degree[queue[k]] < node_degree[o])
k++;
for (int l = queue_end; l > k; --l)
queue[l] = queue[l-1];
queue[k] = o;
queue_end++;
#endif
}
}
}
stage_start = stage_end;
stage_end = queue_end;
}
int final_bandwidth = 0;
for (int i = 0; i < mesh->n_triangles; ++i) {
int *tri = mesh->triangles + i*3;
int m[3] = {ordering[tri[0]], ordering[tri[1]], ordering[tri[2]]};
int d[3] = {abs(m[0]-m[1]), abs(m[0]-m[2]), abs(m[1]-m[2])};
if (d[0] > final_bandwidth) final_bandwidth = d[0];
if (d[1] > final_bandwidth) final_bandwidth = d[1];
if (d[2] > final_bandwidth) final_bandwidth = d[2];
}
printf("n_nodes : %i final bandwidth : %i\n",mesh->n_nodes, final_bandwidth);
free(queue);
free(node_degree);
free(node_connectivity);
return final_bandwidth;
}
//from wikipedia
//INPUT: A - array of pointers to rows of a square matrix having dimension N
// Tol - small tolerance number to detect failure when the matrix is near degenerate
//OUTPUT: Matrix A is changed, it contains both matrices L-E and U as A=(L-E)+U such that P*A=L*U.
// P - array of N+1 integers containing pivoting of A and P[N] is for determinant computation
// (P[i] is an index of the row in A placed at the i-th row
static int LUPDecompose(double *__restrict__*__restrict A,int N,double Tol,int *P, int bw){
int i,j,k;
for(i=0;i<N;i++)
P[i]=i; //integer vector corresponding to no pivoting
//int c0 = 0, c1 = 0;
for(i=0;i<N;i++){
/*double absA;
double maxA=0.0;
int imax=i;
for(k=i;k<N;k++)
if((absA=fabs(A[k][i]))>maxA){ maxA=absA; imax=k; }
if(maxA<Tol)return(0); //failure, matrix is degenerate
if(imax!=i){
j=P[i]; P[i]=P[imax]; P[imax]=j; //pivoting P
double *ptr=A[i]; A[i]=A[imax]; A[imax]=ptr; //pivoting rows of A
}*/
for(j=i+1;j<i+bw && j<N;j++){
if(A[j][i] == 0.)
continue;
A[j][i]/=A[i][i];
for(k=i+1;k<i+bw && k<N;k++){
/*if(fabs(A[i][k]) < 1e-12){
c1++;
//continue;
}*/
A[j][k]-=A[j][i]*A[i][k];
//c0++;
}
}
}
//printf("%g / %i\n", c1*100./c0, c0);
return(1); //decomposition done
}
//INPUT: A,P filled in LUPDecompose; b - rhs vector; N - dimension
//OUTPUT: x - solution vector of A*x=b
static void LUPSolve(double **A,int *P,double *b,int N,double *x, int bw){
for(int i=0;i<N;i++){
x[i]=b[P[i]];
for(int k=(i-bw<0?0:i-bw);k<i;k++)
x[i]-=A[i][k]*x[k];
}
for(int i=N-1;i>=0;i--){
for(int k=i+1;k<N && k<i+bw;k++)
x[i]-=A[i][k]*x[k];
x[i]=x[i]/A[i][i];
}
}
struct LinearSystemStruct{
double *A;
double **rows;
double *b;
double *x;
int *node_map;
int band_width;
int line_width;
int n_fields;
int n;
Mesh *mesh;
int *isfixed;
};
LinearSystem *linear_system_new(Mesh *mesh, int n_fields, int n_boundaries, const StrongBoundary *boundaries)
{
LinearSystem *system = malloc(sizeof(LinearSystem));
system->n_fields = n_fields;
system->mesh = mesh;
system->n = mesh->n_nodes*system->n_fields;
system->node_map = malloc(sizeof(int)*mesh->n_nodes);
int mbw = reverse_cuthill_mckee(mesh, system->node_map);
system->band_width = system->n_fields*(1+mbw);
system->line_width = system->band_width*2+1;
system->A = malloc(sizeof(double)*system->n*system->line_width);
system->rows = malloc(sizeof(double*)*system->n);
for (int i = 0; i < system->n; ++i)
system->rows[i] = system->A + i*(system->line_width-1) + system->band_width;
system->b = malloc(sizeof(double)*system->n);
system->x = malloc(sizeof(double)*system->n);
system->isfixed = malloc(sizeof(int)*system->n);
for (int i = 0; i < system->n; ++i) {
system->isfixed[i] = 0;
}
for (int ibnd = 0; ibnd < n_boundaries; ++ibnd) {
const StrongBoundary *bnd = 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){
system->isfixed[system->node_map[mesh->phys_nodes[iphys][i]]*system->n_fields+bnd->field] = 1;
}
}
return system;
}
void linear_system_add_to_matrix(LinearSystem *system, int el0, int el1, const double *local_matrix){
int *tri0 = &system->mesh->triangles[el0*3];
int *tri1 = &system->mesh->triangles[el1*3];
int nf = system->n_fields;
for (int i = 0; i < 3; ++i) {
for (int inf = 0; inf < 3; ++inf) {
int ii = system->node_map[tri0[i]]*nf + inf;
for (int j = 0; j < 3; ++j) {
for (int jnf = 0; jnf < nf; ++jnf) {
int jj = system->node_map[tri1[j]]*nf + jnf;
system->rows[ii][jj] += local_matrix[(inf*3+i)*nf*3+jnf*3+j];
}
}
}
}
}
void linear_system_add_to_rhs(LinearSystem *system, int el0, const double *local_vector)
{
const int *tri = system->mesh->triangles + el0*3;
int n_fields = system->n_fields;
for (int i = 0; i < n_fields; ++i) {
for (int j = 0; j < 3; ++j) {
int m = system->node_map[tri[j]]*n_fields+i;
system->b[m] += local_vector[i*3+j];
}
}
}
void linear_system_zero_matrix_and_rhs(LinearSystem *system)
{
for (int i = 0; i < system->n; ++i){
system->b[i] = 0.;
for (int j = 0; j < system->line_width; ++j) {
system->A[i*system->line_width+j] = 0.;
}
}
}
void linear_system_solve(LinearSystem *system, double *solution){
int *P = malloc(sizeof(int)*(system->n));
double **rows = system->rows;
for (int i = 0; i < system->n; ++i){
if (system->isfixed[i]) {
for (int j = i-system->band_width; j <= i+system->band_width; ++j)
rows[i][j] = 0.;
rows[i][i] = 1.;
system->b[i] = 0;
}
}
LUPDecompose(rows, system->n, 1e-8, P, system->band_width);
LUPSolve(rows, P, system->b, system->n, system->x, system->band_width);
for (int i = 0; i < system->mesh->n_nodes; ++i){
int ii = system->node_map[i];
for (int j = 0; j < system->n_fields; ++j)
solution[i*system->n_fields+j] = system->x[ii*system->n_fields+j];
}
}
void linear_system_free(LinearSystem *system)
{
free(system->b);
free(system->x);
free(system->A);
free(system->rows);
free(system->isfixed);
free(system->node_map);
free(system);
}
double linear_system_get_rhs_norm(LinearSystem *system)
{
double norm = 0;
for (int i = 0; i < system->n;++i)
if (!system->isfixed[i])
norm += system->b[i]*system->b[i];
return sqrt(norm);
}
void initialize_linear_solver(int argc, char **argv){}
#include <math.h>
#include "mesh.h"
#include "linear_system.h"
#include <stdlib.h>
#include <string.h>
#define CONMAX 12
void connectivity_insert(int *connectivity, int i, int j)
{
for (int k = 0; k < CONMAX; ++k) {
int *p = connectivity + CONMAX*i + k;
if (*p == -1)
*p = j;
if (*p == j)
return;
}
printf("ERROR : node %i has more than %i neighbours\n", i, CONMAX);
}
int reverse_cuthill_mckee(Mesh *mesh, int *ordering)
{
int *node_connectivity = malloc(sizeof(int)*mesh->n_nodes*CONMAX);
for (int i = 0; i < mesh->n_nodes*CONMAX; ++i) {