Commit 7a36fe24 authored by Michel Henry's avatar Michel Henry
Browse files

wip

parent 56557cd2
Pipeline #9971 canceled with stages
......@@ -272,7 +272,7 @@ void fluid_problem_compute_stab_parameters(FluidProblem *problem, double dt) {
amax = fmax(amax,problem->concentration[iel*N_N+i]);
amin = fmin(amin,problem->concentration[iel*N_N+i]);
}
for(int j = 0; j < DIMENSION; ++j) normup += pow2(problem->solution[el[i]*n_fields+j]);
for(int j = 0; j < DIMENSION; ++j) normup += pow2(problem->solution[el[i]*n_fields+j]);//to change
normu = fmax(normu,sqrt(normup));
}
double dxidx[DIMENSION][DIMENSION];
......@@ -414,19 +414,22 @@ static void fluid_problem_f(const FluidProblem *problem, const double *sol, doub
double *g = problem->g;
double rhoreduced = (problem->reduced_gravity ? (rho-problem->rho[0]) : rho);
double drag = problem->volume_drag + problem->quadratic_drag*nold;
// f0[P] = p;
// f00[P*n_fields+P] = 1;
for (int i = 0; i < D; ++i) {
f0[U+i] =
+ c*dp[i]
- g[i]*rhoreduced*c
- bf[i]*c
+ drag*u[i];
//+ 5.3e5*u[i];
+ 5.3e5*u[i];
f00[(U+i)*n_fields+U+i] += drag;//5.3e5;
f01[(U+i)*n_fields+P][i] += c;
if(problem->temporal) {
f0[U+i] += rho*(u[i]-uold[i])/dt;
f0[U+i] += rho*(u[i]-uold[i])/dt;// - rho;
f00[(U+i)*n_fields+U+i] += rho/dt;
}
// #if 0
if (problem->advection) {
// advection :
// dj(uj ui /c) = uj/c dj(ui) + ui/c dj(uj) - uj ui /(c*c) dj(c)
......@@ -469,6 +472,7 @@ static void fluid_problem_f(const FluidProblem *problem, const double *sol, doub
for (int j = 0; j < D; ++j) {
f11[P*n_fields+U+i][i][j] = pspg*f01[(U+i)*n_fields+(U+i)][j];
}
// #endif
}
}
......@@ -726,7 +730,7 @@ static void fluid_problem_volume(FluidProblem *problem, const double *solution_o
const double detj = mesh_dxidx(mesh,iel,dxidx);
double *local_vector = &all_local_vector[fields->local_size*iel];
double *local_matrix = &all_local_matrix[fields->local_size*fields->local_size*iel];
for (int iqp = 0; iqp< N_QP; ++iqp) {
for (int iqp = 0; iqp < N_QP; ++iqp) {
double s[fields->n], sold[fields->n], ds[fields->n][D], dsold[fields->n][D];
double sf[fields->local_size], dsf[fields->local_size][D];
fe_fields_f(fields, QP[iqp], sf);
......@@ -1028,8 +1032,10 @@ void fluid_problem_assemble_system(FluidProblem *problem, const double *solution
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)
// char *forced_row = malloc(sizeof(char)*n_fields*mesh->n_nodes);
int ndof = fluid_problem_get_ndof(problem->fields, problem->mesh);
char *forced_row = malloc(sizeof(char)*ndof);
for (int i = 0; i < ndof; ++i)
forced_row[i] = -1;
for (int ibnd = 0; ibnd < problem->n_strong_boundaries; ++ibnd) {
const StrongBoundary *bnd = problem->strong_boundaries + ibnd;
......@@ -1055,25 +1061,27 @@ void fluid_problem_assemble_system(FluidProblem *problem, const double *solution
fluid_problem_surface_tension(problem, solution_old, problem->grad_a_cg, all_local_vector);
free(grad_a_cg);
}
for (int iel=0; iel < mesh->n_elements; ++iel) {
double *local_vector = &all_local_vector[local_size*iel];
double *local_matrix = &all_local_matrix[local_size*local_size*iel];
for (int inode = 0; inode < element->nlocal; ++inode) {
for(int ifield = 0; ifield < n_fields; ++ifield) {
const int *el = &mesh->elements[iel*N_N];
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< element->nlocal; ++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*element->nlocal] = 0;
}
}
}
// to do
// for (int iel=0; iel < mesh->n_elements; ++iel) {
// double *local_vector = &all_local_vector[local_size*iel];
// double *local_matrix = &all_local_matrix[local_size*local_size*iel];
// for (int inode = 0; inode < element->nlocal; ++inode) {
// for(int ifield = 0; ifield < n_fields; ++ifield) {
// const int *el = &mesh->elements[iel*N_N];
// // char forced_field = forced_row[el[inode]*n_fields+ifield];
// char forced_field = -1;
// if (forced_field == -1) continue;
// int i = inode*n_fields + ifield;
// for (int jnode = 0; jnode< element->nlocal; ++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*element->nlocal] = 0;
// }
// }
// }
free(forced_row);
}
......@@ -1117,9 +1125,7 @@ void fluid_problem_apply_boundary_conditions(FluidProblem *problem)
}
}
bnd->apply(ndof_closure, x, v);
int map[problem->fields->local_size];
fe_fields_local_map(problem->fields, mesh, iel, map);
int *map = &(problem->fields->map[iel*problem->fields->local_size]);
int imap=0;
for(int i=0; problem->fields->n; ++i){
if(i == bnd->field) break;
......@@ -1195,8 +1201,7 @@ static void fluid_problem_compute_node_volume(FluidProblem *problem) {
free(periodic_node_volume);
}
static void compute_porosity(Mesh *mesh, double *node_volume, double *porosity, int n_particles, double *particle_position, double *particle_volume, int *particle_element_id, double *particle_uvw, int *flag)
{
static void compute_porosity(Mesh *mesh, double *node_volume, double *porosity, int n_particles, double *particle_position, double *particle_volume, int *particle_element_id, double *particle_uvw, int *flag) {
#if DIMENSION == 3
FEElement *element = fe_element_new("tetrahedron_p1");
#else
......@@ -1232,21 +1237,25 @@ static void compute_porosity(Mesh *mesh, double *node_volume, double *porosity,
}
}
FluidProblem *fluid_problem_new(double *g, int n_fluids, double *mu, double *rho, double sigma, double coeffStab, double volume_drag, double quadratic_drag, int drag_in_stab, double drag_coeff_factor, double ip_factor, int temporal, int advection)
{
FluidProblem *fluid_problem_new(double *g, int n_fluids, double *mu, double *rho, double sigma, double coeffStab, double volume_drag, double quadratic_drag, int drag_in_stab, double drag_coeff_factor, double ip_factor, int temporal, int advection, int* order_f) {
if (n_fluids != 1 && n_fluids != 2) {
printf("error : only 1 or 2 fluids are supported\n");
}
FluidProblem *problem = (FluidProblem*)malloc(sizeof(FluidProblem));
#if DIMENSION == 3
const char *types[] = {"tetrahedron_p1","tetrahedron_p1","tetrahedron_p1","tetrahedron_p1",NULL};
const int p1p1[4] = {1,1,1,1};
const int isp1 = order_f[0]==p1p1[0];
const char *types_fp1[] = {"tetrahedron_p1","tetrahedron_p1","tetrahedron_p1","tetrahedron_p1",NULL};
const char *types_fp2[] = {"tetrahedron_p2","tetrahedron_p2","tetrahedron_p2","tetrahedron_p1",NULL};
const char *type_p[] = {"tetrahedron_p1",NULL};
#else
const char *types[] = {"triangle_p1","triangle_p1", "triangle_p1",NULL};
// const char *types[] = {"triangle_p2","triangle_p2", "triangle_p1",NULL};
const int p1p1[3] = {1,1,1};
int isp1 = order_f[0]==p1p1[0];
const char *types_fp1[] = {"triangle_p1","triangle_p1", "triangle_p1",NULL};
const char *types_fp2[] = {"triangle_p2","triangle_p2", "triangle_p1",NULL};
const char *type_p[] = {"triangle_p1",NULL};
#endif
problem->fields = fe_fields_new(types);
problem->fields = fe_fields_new(isp1?types_fp1:types_fp2);
problem->field_porosity = fe_fields_new(type_p);
const char *type_a[] = {"triangle_p1dg",NULL};
problem->field_concentration = fe_fields_new(type_a);
......@@ -1585,6 +1594,37 @@ void fluid_problem_adapt_gen_mesh(FluidProblem *problem, double lcmax, double lc
#endif
}
int fluid_problem_get_ndof(const FEFields *fields, const Mesh *mesh){
int c = -1;
for(int iel=0; iel<mesh->n_elements; ++iel){
int *map = &(fields->map[iel*fields->local_size]);
for(int j = 0; j < fields->local_size; ++j){
c = map[j] > c ? map[j]:c;
}
}
if (c == -1){
printf("Error : ndof = -1\n");
exit(1);
}
return ++c;
}
void fluid_problem_set_problem(FluidProblem *problem) {
free(problem->porosity);
int ndof = fluid_problem_get_ndof(problem->field_porosity, problem->mesh);
problem->porosity = malloc(ndof*sizeof(double));
free(problem->oldporosity);
problem->oldporosity = malloc(ndof*sizeof(double));
for(int i = 0; i < ndof; ++i){
problem->porosity[i] = 1.;
problem->oldporosity[i] = 1.;
}
free(problem->solution);
int n_fields = fluid_problem_n_fields(problem);
ndof = fluid_problem_get_ndof(problem->fields, problem->mesh);
problem->solution = calloc(ndof, sizeof(double));
}
// allocate all mesh-dependant structures
void fluid_problem_set_mesh(FluidProblem *problem, Mesh *mesh) {
if(problem->mesh)
......@@ -1597,27 +1637,29 @@ void fluid_problem_set_mesh(FluidProblem *problem, Mesh *mesh) {
problem->bulk_force = malloc(sizeof(double)*mesh->n_nodes*DIMENSION);
for (int i=0; i<mesh->n_nodes*DIMENSION; ++i)
problem->bulk_force[i] = 0;
free(problem->porosity);
problem->porosity = (double*)malloc(mesh->n_nodes*sizeof(double));
// free(problem->porosity);
// problem->porosity = (double*)malloc(mesh->n_nodes*sizeof(double));
free(problem->grad_a_cg);
problem->grad_a_cg = malloc(sizeof(double)*mesh->n_nodes*D);
problem->all_kappa = malloc(sizeof(double)*mesh->n_elements);
free(problem->oldporosity);
problem->oldporosity = (double*)malloc(mesh->n_nodes*sizeof(double));
for(int i = 0; i < mesh->n_nodes; ++i){
problem->porosity[i] = 1.;
problem->oldporosity[i] = 1.;
}
// free(problem->oldporosity);
// problem->oldporosity = (double*)malloc(mesh->n_nodes*sizeof(double));
// for(int i = 0; i < mesh->n_nodes; ++i){
// problem->porosity[i] = 1.;
// problem->oldporosity[i] = 1.;
// }
free(problem->boundary_force);
problem->boundary_force = (double*)malloc(DIMENSION*mesh->n_boundaries*sizeof(double));
for(int i = 0; i < mesh->n_boundaries*DIMENSION; i++) problem->boundary_force[i] = 0;
free(problem->solution);
int n_fields = fluid_problem_n_fields(problem);
problem->solution = (double*)malloc(mesh->n_nodes*n_fields*sizeof(double));
// free(problem->solution);
// int n_fields = fluid_problem_n_fields(problem);
// int ndof = fluid_problem_get_ndof(problem->fields, problem->mesh);
// problem->solution = calloc(ndof, sizeof(double));
// problem->solution = (double*)malloc(mesh->n_nodes*n_fields*sizeof(double));// to change
if (problem->n_fluids == 2) {
free(problem->concentration);
problem->concentration = (double*)malloc(mesh->n_elements*N_N*sizeof(double));
......@@ -1627,9 +1669,9 @@ void fluid_problem_set_mesh(FluidProblem *problem, Mesh *mesh) {
// for(int i = 0; i < mesh->n_nodes*D; ++i) problem->mesh_velocity[i] = 0;
problem->taup = malloc(sizeof(double)*mesh->n_elements);
for (int i = 0; i < problem->mesh->n_elements; ++i) problem->taup[i] = 0;
for (int i = 0; i < problem->mesh->n_nodes*n_fields; ++i){
problem->solution[i] = 0.;
}
// for (int i = 0; i < problem->mesh->n_nodes*n_fields; ++i){ // to change
// problem->solution[i] = 0.;
// }
fluid_problem_compute_node_volume(problem);
if (problem->boundaries) mesh_boundaries_free(problem->boundaries,problem->mesh->n_boundaries);
problem->boundaries = mesh_boundaries_create(problem->mesh);
......@@ -1765,7 +1807,7 @@ void fluid_problem_set_particles(FluidProblem *problem, int n, double *mass, dou
for (int i = 0; i < problem->mesh->n_nodes; ++i) {
problem->oldporosity[i] += problem->porosity[i]-c[i];
for (int d = 0; d < DIMENSION; ++d) {
problem->solution[n_fields*i+U+d] *= problem->porosity[i]/c[i];
problem->solution[n_fields*i+U+d] *= problem->porosity[i]/c[i]; // to change
}
}
free(c);
......@@ -1804,7 +1846,7 @@ void fluid_problem_move_particles(FluidProblem *problem, int n, double *position
int nf = fluid_problem_n_fields(problem);
for (int i = 0; i < problem->mesh->n_nodes; i++) {
for (int d = 0; d < DIMENSION; ++d) {
problem->solution[i*nf+U+d] *= problem->porosity[i]/(problem->porosity[i] - particle_enter_porosity[i] + particle_leave_porosity[i]);
problem->solution[i*nf+U+d] *= problem->porosity[i]/(problem->porosity[i] - particle_enter_porosity[i] + particle_leave_porosity[i]); // to change
}
}
free(particle_enter_porosity);
......
......@@ -111,7 +111,7 @@ void fluid_problem_move_particles(FluidProblem *problem, int n, double *position
void fluid_problem_set_bulk_force(FluidProblem *problem, double *force);
void fluid_problem_set_concentration_cg(FluidProblem *problem, double *concentration);
void fluid_problem_free(FluidProblem *problem);
FluidProblem *fluid_problem_new(double *g, int n_fluids, double *mu, double *rho, double sigma, double coeffStab, double volume_drag, double quadratic_drag, int drag_in_stab, double drag_coeff_factor, double ip_factor, int temporal, int advection);
FluidProblem *fluid_problem_new(double *g, int n_fluids, double *mu, double *rho, double sigma, double coeffStab, double volume_drag, double quadratic_drag, int drag_in_stab, double drag_coeff_factor, double ip_factor, int temporal, int advection, int* order_f);
void fluid_problem_adapt_mesh(FluidProblem *problem, Mesh *new_mesh, int old_n_particles, double *old_particle_position, double *old_particle_volume);
void fluid_problem_set_mesh(FluidProblem *problem, Mesh *mesh);
int *fluid_problem_particle_element_id(FluidProblem *problem);
......@@ -133,7 +133,8 @@ void fluid_problem_mesh_local_map(FluidProblem *p, int *idx);
void fluid_problem_concentration_local_map(FluidProblem *p, int *idx);
void fluid_problem_porosity_local_map(FluidProblem *p, int *idx);
void fluid_problem_local_map(FluidProblem *p, int *idx);
int fluid_problem_get_ndof(const FEFields *fields, const Mesh *mesh);
void fluid_problem_set_problem(FluidProblem *problem);
// mesh data (for i/o)
int fluid_problem_n_nodes(const FluidProblem *p);
......
......@@ -36,13 +36,13 @@ void fe_fields_eval(const FEFields *fields, const double *xi, const Mesh *mesh,
void fe_fields_local_map(const FEFields *fields, const Mesh *mesh, int iel, int map[]);
void fe_fields_add_to_local_vector(const FEFields *fields, const double *f0, const double f1[][D], const double *sf, const double dsf[][D], double jw, double *local_vector);
static void fe_fields_add_to_local_matrix(const FEFields *fields, const double *f00, const double f01[][D], const double f10[][D], const double f11[][D][D], const double *sf, const double dsf[][D], double jw, double *local_matrix) {
int jc = 0;
int n_fields = fields->n;
for (int jfield = 0; jfield < fields->n; ++jfield) {
for (int jphi = 0; jphi < fields->element[jfield]->nlocal; ++jphi){
int ic = 0;
for (int ifield = 0; ifield < fields->n; ++ifield){
for (int iphi = 0; iphi < fields->element[ifield]->nlocal; ++iphi){
int c = 0;
int ic = 0;
for(int ifield=0; ifield < fields->n; ++ifield){
for(int iphi = 0; iphi < fields->element[ifield]->nlocal; ++iphi){
int jc = 0;
for(int jfield=0; jfield < fields->n; ++jfield){
for(int jphi=0; jphi < fields->element[jfield]->nlocal; ++jphi){
double m = 0;
if (f00) {
m += jw*sf[ic]*sf[jc]*f00[ifield*fields->n+jfield];
......@@ -64,24 +64,66 @@ static void fe_fields_add_to_local_matrix(const FEFields *fields, const double *
m += jw*sf[ic]*dsf[jc][jd]*f01[ifield*fields->n+jfield][jd];
}
}
local_matrix[(ifield*n_fields+iphi)*n_fields*n_fields+jfield*n_fields+jphi] += m;
ic ++;
local_matrix[c++] += m;
jc++;
}
}
jc ++;
ic++;
}
}
// int c = 0;
// int jc = 0;
// int n_fields = fields->n;
// for (int jfield = 0; jfield < fields->n; ++jfield) {
// for (int jphi = 0; jphi < fields->element[jfield]->nlocal; ++jphi){
// int ic = 0;
// for (int ifield = 0; ifield < fields->n; ++ifield){
// for (int iphi = 0; iphi < fields->element[ifield]->nlocal; ++iphi){
// double m = 0;
// if (f00) {
// m += jw*sf[ic]*sf[jc]*f00[ifield*fields->n+jfield];
// }
// if (f10) {
// for (int id = 0; id < D; ++id) {
// m += jw*dsf[ic][id]*sf[jc]*f10[ifield*fields->n+jfield][id];
// }
// }
// if (f11) {
// for (int id = 0; id < D; ++id) {
// for (int jd = 0; jd < D; ++jd) {
// m += jw*dsf[ic][id]*dsf[jc][jd]*f11[ifield*fields->n+jfield][id][jd];
// }
// }
// }
// if (f01) {
// for (int jd = 0; jd < D; ++jd) {
// m += jw*sf[ic]*dsf[jc][jd]*f01[ifield*fields->n+jfield][jd];
// }
// }
// local_matrix[c++] += m;
// ic ++;
// }
// }
// jc ++;
// }
// }
}
void fe_fields_free(FEFields *fields);
#if D==2
#define N_QP 3
// #define N_QP 3
#define N_QP 13
#define N_LQP 2
static double LQP[][1] = {{0.21132486540518713}, {0.7886751345948129}};
static double LQW[] = {0.5,0.5};
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};
// 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};
static double QW[] = {-0.0747850222338409948,0.087807628716603997,0.087807628716603997,0.087807628716603997,0.0266736178044190003,0.0266736178044190003,0.0266736178044190003,0.0385568804451279987,0.0385568804451279987,0.0385568804451279987,0.0385568804451279987,0.0385568804451279987,0.0385568804451279987};
static double QP[][2] = {{0.333333333333332982,0.333333333333332982},{0.260345966079039981,0.260345966079039981},{0.260345966079039981,0.479308067841919982},{0.479308067841919982,0.260345966079039981},{0.0651301029022160055,0.0651301029022160055},{0.0651301029022160055,0.869739794195568017},{0.869739794195568017,0.0651301029022160055},{0.312865496004874011,0.638444188569810001},{0.638444188569810001,0.0486903154253160025},{0.0486903154253160025,0.312865496004874011},{0.312865496004874011,0.0486903154253160025},{0.638444188569810001,0.312865496004874011},{0.0486903154253160025,0.638444188569810001}};
static double detDD(const double m[2][2])
{
......
......@@ -28,6 +28,7 @@
The continuous phase (mixture) is computed by solving averaged Navier-Stokes equations on unstructured meshes with the continuous finite element method.
"""
from builtins import print
from ctypes import *
from time import time
import numpy as np
......@@ -127,7 +128,7 @@ def _load_msh(mesh_file_name, lib, dim) :
class FluidProblem :
"""Creates the numerical representation of the fluid."""
def __init__(self, dim, g, mu, rho, sigma=0, coeff_stab=0.01, volume_drag=0., quadratic_drag=0., solver=None, solver_options="", drag_in_stab=1, drag_coefficient_factor=1, ip_factor=1, temporal=True, advection=True):
def __init__(self, dim, g, mu, rho, sigma=0, coeff_stab=0.01, volume_drag=0., quadratic_drag=0., solver=None, solver_options="", drag_in_stab=1, drag_coefficient_factor=1, ip_factor=1, temporal=True, advection=True, order=[1,1,1]):
"""Builds the fluid structure.
Keyword arguments:
......@@ -155,6 +156,9 @@ class FluidProblem :
self.weak_cb_ref = []
self.sys = None
self.mean_pressure = None
self.order_fields = order
self.order_porosity = [1]
self.order_concentration = self.order_porosity
if dim == 2 :
self._lib = lib2
elif dim == 3 :
......@@ -165,7 +169,11 @@ class FluidProblem :
n_fluids = np.require(rho,"float","C").reshape([-1]).shape[0]
self._n_fluids = n_fluids
self._dim = dim
self._fp = c_void_p(self._lib.fluid_problem_new(_np2c(g), n_fluids, _np2c(mu), _np2c(rho), c_double(sigma), c_double(coeff_stab), c_double(volume_drag), c_double(quadratic_drag), c_int(drag_in_stab),c_double(drag_coefficient_factor),c_double(ip_factor),c_int(temporal), c_int(advection)))
self._fp = c_void_p(self._lib.fluid_problem_new(
_np2c(g), n_fluids, _np2c(mu), _np2c(rho), c_double(sigma), c_double(coeff_stab),
c_double(volume_drag), c_double(quadratic_drag), c_int(drag_in_stab),c_double(drag_coefficient_factor),
c_double(ip_factor), c_int(temporal), c_int(advection),_np2c(np.array(self.order_fields, dtype=np.int32), dtype=np.int32)
))
if self._fp == None :
raise NameError("Cannot create fluid problem.")
......@@ -186,8 +194,15 @@ class FluidProblem :
"""
mesh = _load_msh(mesh_file_name, self._lib, self.dimension())
self._lib.fluid_problem_set_mesh(self._fp, mesh)
idx = self.set_numbering([1,1])
idx = self.set_numbering([1]*self.dimension())
self._lib.fluid_problem_mesh_local_map(self._fp, _np2c(idx, dtype=np.int32))
idx = self.set_numbering(self.order_fields)
idx_concentration = self.set_numbering(self.order_concentration)
idx_porosity = self.set_numbering(self.order_porosity)
self._lib.fluid_problem_fields_local_map(self._fp, _np2c(idx, dtype=np.int32))
self._lib.fluid_problem_porosity_local_map(self._fp, _np2c(idx_porosity, dtype=np.int32))
self._lib.fluid_problem_concentration_local_map(self._fp, _np2c(idx_concentration, dtype=np.int32))
self._lib.fluid_problem_set_problem(self._fp)
self.sys = None
gmsh.model.remove()
......@@ -366,7 +381,7 @@ class FluidProblem :
offsets = np.cumsum(np.repeat([self._dim+1],connectivity.shape[0])).astype(np.int32)
VTK.write(output_dir+"/fluid",it,t,(types,connectivity,offsets),self.coordinates(),data,field_data,cell_data)
def read_vtk(self, dirname, i):
def read_vtk(self, dirname, i): # to update for p2p1
"""Reads output file to reload computation.
Keyword arguments:
......@@ -435,18 +450,18 @@ class FluidProblem :
self._lib.fluid_problem_apply_boundary_conditions(self._fp)
self._lib.fluid_problem_reset_boundary_force(self._fp)
self._lib.fluid_problem_compute_stab_parameters(self._fp,c_double(dt))
idx = self.set_numbering([1,1,1])
idx_p1 = self.set_numbering([1])
idx = self.set_numbering(self.order_fields)
idx_concentration = self.set_numbering(self.order_concentration)
idx_porosity = self.set_numbering(self.order_porosity)
self._lib.fluid_problem_fields_local_map(self._fp, _np2c(idx, dtype=np.int32))
self._lib.fluid_problem_porosity_local_map(self._fp, _np2c(idx_p1, dtype=np.int32))
self._lib.fluid_problem_concentration_local_map(self._fp, _np2c(idx_p1, dtype=np.int32))
self._lib.fluid_problem_porosity_local_map(self._fp, _np2c(idx_porosity, dtype=np.int32))
self._lib.fluid_problem_concentration_local_map(self._fp, _np2c(idx_concentration, dtype=np.int32))
# print(time()-timer)
# idx = np.zeros((self.n_elements()*self.local_size()), dtype=np.int32)
# self.local_map(idx)
idx = idx.reshape([self.n_elements(), -1])
periodic_map = self._get_matrix("periodic_mapping", self.n_nodes(), 1, c_int).flatten()
solution_old = self.solution().copy()
solution_old = self.solution_flat().copy()
if self.sys is None :
constraints = []
if self.mean_pressure is not None:
......@@ -461,13 +476,19 @@ class FluidProblem :
localm = np.ndarray([sys.localsize**2*self.n_elements()])
self._lib.fluid_problem_assemble_system(self._fp,_np2c(solution_old),c_double(dt),_np2c(localv),_np2c(localm))
# print(time()-timer)
# np.set_printoptions(precision=1,threshold=1e-8)
# mass_matrix = localm[:sys.localsize**2].reshape(sys.localsize, -1)
# print(360*mass_matrix[:6,:6])
constraint_value = []
if self.mean_pressure is not None:
constraint_value.append((self.node_volume().reshape([-1]), -self.mean_pressure*np.sum(self.node_volume())))
rhs = sys.local_to_global(localv,localm,self.solution().reshape([-1]),constraint_value)
self.solution()[:] -= sys.solve(rhs).reshape([-1, self.n_fields()])[periodic_map]
rhs = sys.local_to_global(localv,localm,self.solution_flat(),constraint_value)
# self.solution()[:] -= sys.solve(rhs).reshape([-1, self.n_fields()])[periodic_map]
# self.solution()[:] -= sys.solve(rhs).reshape([-1, self.n_fields()])
self.solution_flat()[:] -= sys.solve(rhs)
# exit(0)
# print(time()-timer)
if check_residual_norm > 0 :
if check_residual_norm > 0 and False :
self._lib.fluid_problem_assemble_system(self._fp,_np2c(solution_old),c_double(dt),_np2c(localv),_np2c(localm))
rhs = sys.local_to_global(localv,localm,self.solution().reshape([-1]),constraint_value)
norm = np.linalg.norm(rhs)
......@@ -477,7 +498,6 @@ class FluidProblem :
self._lib.fluid_problem_node_force_volume(self._fp,_np2c(solution_old),c_double(dt),None,None)
# print(time()-timer)
def compute_cfl(self):
"""Computes the CFL number divided by the time step
......@@ -532,7 +552,14 @@ class FluidProblem :
def solution(self) :
"""Gives access to the fluid field value at the mesh nodes."""
return self._get_matrix("solution", self.n_nodes(), self.n_fields())
sol = self.solution_flat()[:self.n_fields()*self.n_nodes()].reshape(-1, self.n_fields())
return sol
# return self._get_matrix("solution", self.n_nodes(), self.n_fields())
def solution_flat(self) :
"""Gives access to the fluid field value at each degree of freedom"""
ndof = np.max(self.set_numbering(self.order_fields))+1
return self._get_matrix("solution", ndof, 1).reshape([-1])
def velocity(self) :
"""Gives access to the fluid velocity value at the mesh nodes."""
......@@ -654,6 +681,32 @@ class FluidProblem :
self._lib.fluid_problem_local_map(self._fp, _np2c(idx, dtype=np.int32))
@timeit
# def set_numbering(self, order=[]):
# mesh_cl = MeshCl(self.elements())
# maps = [LagrangeSpace(mesh_cl, o).element_nodes for o in order]
# ndofs = [lmap.shape[1] for lmap in maps]
# ndofs_unique = np.unique(ndofs)
# n_order = len(order)
# # map en u0,v0,p0, u1,v1,p1, ...
# idx_uvp = np.full([self.n_elements(), np.sum(ndofs)], -1, dtype=np.int32)
# for i,ndof in enumerate(ndofs_unique):
# ndof_by_order = ndof if i == 0 else ndof-ndofs_unique[i-1]
# offset_dof = 0 if i == 0 else idx_uvp.max()+1
# offset_ind = np.where(idx_uvp[0]==-1)[0][0]
# order_to_use = np.where(ndofs >= ndof)[0]
# offset_order = 0 if i == 0 else ndof-ndofs_unique[i-1]
# n_order = len(order_to_use)
# for j in range(ndof_by_order):
# for k in range(n_order):
# lmap = maps[order_to_use[k]][:,offset_order:] - maps[order_to_use[k]][:,offset_order:].min()
# idx_uvp[:,offset_ind + j*n_order+k] = n_order*lmap[:,j]+k+offset_dof
# # map en u0,u1, ..., v0, v1, ..., p0, p1, ...
# idx = np.full([self.n_elements(), np.sum(ndofs)], -1, dtype=np.int32)
# n = ndofs_unique[0]
# for i in range(n_order):
# idx[:,i*n:(i+1)*n] = idx_uvp[:,i::n_order]
# return idx.reshape([-1])
def set_numbering(self, order=[]):
mesh_cl = MeshCl(self.elements())
maps = [LagrangeSpace(mesh_cl, o).element_nodes for o in order]
......@@ -674,10 +727,14 @@ class FluidProblem :
lmap = maps[order_to_use[k]][:,offset_order:] - maps[order_to_use[k]][:,offset_order:].min()
idx_uvp[:,offset_ind + j*n_order+k] = n_order*lmap[:,j]+k+offset_dof
# map en u0,u1, ..., v0, v1, ..., p0, p1, ...
idx = np.full([self.n_elements(), np.sum(ndofs)], -1, dtype=np.int32)
n = ndofs_unique[0]
for i in range(n_order):
idx[:,i*n:(i+1)*n] = idx_uvp[:,i::n_order]
idx = np.full([self.n_elements(),np.sum(ndofs)],-1,dtype=np.int32)
if np.all(np.array(order)==1):
permutation = np.repeat(n_order*np.arange(ndofs[0])[None,:], n_order, axis=0).reshape([-1]) \
+ np.repeat(np.arange(n_order), ndofs[0])
idx = idx_uvp[:,permutation]
if order==[2,2,1]: # 2D
permutation = np.array([0,3,6,9,11,13,1,4,7,10,12,14,2,5,8], dtype=np.int32)
idx = idx_uvp[:,permutation]
return idx.reshape([-1])
def _n_physicals(self):
......
......@@ -102,7 +102,7 @@ ii = 0
#
# FLUID PROBLEM
#
fluid = fluid.FluidProblem(2,g,[nu*rho],[rho])
fluid = fluid.FluidProblem(2,g,[nu*rho],[rho], solver="petsc4py")
# Set the mesh geometry for the fluid computation
fluid.load_msh("mesh.msh")
fluid.set_wall_boundary("Bottom")
......