Skip to content
GitLab
Commit 9fe1e3db authored by Michel Henry's avatar Michel Henry
Browse files

wip, pressure wrong at factor -2 in global matrix

parent 7810bf04
Pipeline #9873 failed with stages
in 2 minutes and 41 seconds
Hide whitespace changes
Inline Side-by-side
fluid/fluid_problem.c
View file @ 9fe1e3db
......@@ -441,10 +441,8 @@ static void fluid_problem_f(const FluidProblem *problem, const double *sol, doub
}
for (int j = 0; j < D; ++j) {
f1[U+i][j] = mu*(du[i][j]-u[i]*dc[j]/c+ du[j][i]-u[j]*dc[i]/c);
// ca ne correspond pas aux bons indices
f10[(U+i)*n_fields+U+j][j] += -mu*dc[i]/c;
f10[(U+i)*n_fields+U+i][i] += -mu*dc[j]/c;
// equivalent a [(U+i)*n_fields+U+j]*D + j != (U+i)*D+j)*n_fields + U+j
f11[(U+i)*n_fields+U+j][j][i] += mu;
f11[(U+i)*n_fields+U+i][j][j] += mu;
......@@ -462,7 +460,6 @@ static void fluid_problem_f(const FluidProblem *problem, const double *sol, doub
f11[(U+i)*n_fields+(U+j)][i][j] += lsic;
}
f1[P][i] = -u[i];
f10[P*n_fields+U+i][i] += -1;
// PSPG
double pspg = taup/rho;
......@@ -1106,8 +1103,11 @@ void fluid_problem_apply_boundary_conditions(FluidProblem *problem)
x[i*D+k] = mesh->x[j*3+k];
}
bnd->apply(mbnd->n_nodes, x, v);
// int map[problem->fields->local_size];
// fe_fields_local_map(problem->fields, problem->mesh, iel, map);
for (int i = 0; i < mbnd->n_nodes; ++i)
problem->solution[mbnd->nodes[i]*n_fields+bnd->field] = v[i];
problem->solution[mbnd->nodes[i]*n_fields+bnd->field] = v[i]; // pas correcte
free(x);
free(v);
}
......@@ -1610,7 +1610,7 @@ void fluid_problem_adapt_mesh(FluidProblem *problem, Mesh *new_mesh, int old_n_p
for (int j=0; j<n_fields; ++j) {
new_solution[i*n_fields+j] = 0;
for (int k = 0; k < element->nlocal; ++k)
new_solution[i*n_fields+j] += problem->solution[el[k]*n_fields+j]*phi[k];
new_solution[i*n_fields+j] += problem->solution[el[k]*n_fields+j]*phi[k]; // pas juste
}
if (problem->n_fluids == 2) {
concentration_cg[i] = 0;
......@@ -1972,6 +1972,7 @@ void fluid_problem_reset_boundary_force(FluidProblem *problem) {
int fluid_problem_local_size(FluidProblem *problem) {
return problem->fields->local_size;
}
// int fluid_problem_global_size(FluidProblem *problem) {
// int n_nodes = problem->mesh->n_nodes;
// int n_edges = problem->mesh->n_edges;
......@@ -1979,8 +1980,14 @@ int fluid_problem_local_size(FluidProblem *problem) {
// int n_elt = problem->mesh->n_elements;
// int global_size = problem->fields->nd[0]*n_nodes
// + problem->fields->nd[1]*n_edges
// + problem->fields->nd[2]*n_edges;
// + problem->fields->nd[2]*n_faces;
// #if DIMENSION == 3
// global_size += problem->fields->nd[3]*n_elt;
// #endif
// }
void fluid_problem_local_map(FluidProblem *p, int map[]){
for(int iel = 0; iel < p->mesh->n_elements; ++iel){
fe_fields_local_map(p->fields, p->mesh, iel, &(map[iel*p->fields->local_size]));
}
}
\ No newline at end of file
fluid/fluid_problem.h
View file @ 9fe1e3db
......@@ -127,8 +127,8 @@ double *fluid_problem_porosity(const FluidProblem *p);
double *fluid_problem_concentration_dg(const FluidProblem *p);
double *fluid_problem_element_size(const FluidProblem *p);
void fluid_problem_advance_concentration(FluidProblem *problem, double dt);
int fluid_problem_local_size(FluidProblem *problem);
void fluid_problem_local_map(FluidProblem *p, int map[]);
int fluid_problem_local_size(FluidProblem *p);
// mesh data (for i/o)
......
fluid/fluid_problem_fem.c
View file @ 9fe1e3db
......@@ -290,4 +290,4 @@ void fe_fields_eval(const FEFields *fields, const double *xi, const Mesh *mesh,
double sf[fields->local_size];
fe_fields_f(fields, xi, sf);
fe_fields_eval_sf(fields, mesh, iel, sf, v, s);
}
}
\ No newline at end of file
fluid/fluid_problem_fem.h
View file @ 9fe1e3db
......@@ -29,9 +29,9 @@ void fe_fields_eval_grad_sf(const FEFields *fields, const Mesh *mesh, int iel, c
void fe_fields_eval_sf(const FEFields *fields, const Mesh *mesh, int iel, const double sf[], const double *v, double s[]);
void fe_fields_eval_grad(const FEFields *fields, const double *xi, const double dxidx[D][D], const Mesh *mesh, int iel, const double *v, double s[], double ds[][DIMENSION]);
void fe_fields_eval(const FEFields *fields, const double *xi, const Mesh *mesh, int iel, const double *v, double s[]);
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) {
printf("ndof_local = %6d\n", fields->local_size); // dof_local
int jc = 0;
int n_fields = fields->n;
for (int jfield = 0; jfield < fields->n; ++jfield) {
......@@ -40,7 +40,6 @@ static void fe_fields_add_to_local_matrix(const FEFields *fields, const double *
for (int ifield = 0; ifield < fields->n; ++ifield){
for (int iphi = 0; iphi < fields->element[ifield]->nlocal; ++iphi){
double m = 0;
// int index = ifield*fields->n+jfield;
if (f00) {
m += jw*sf[ic]*sf[jc]*f00[ifield*fields->n+jfield];
}
......@@ -61,7 +60,6 @@ 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[] += m; // elt*dof_local
local_matrix[(ifield*n_fields+iphi)*n_fields*n_fields+jfield*n_fields+jphi] += m;
ic ++;
}
......@@ -69,15 +67,6 @@ static void fe_fields_add_to_local_matrix(const FEFields *fields, const double *
jc ++;
}
}
// for(int jfield=0; jfields < fields->n; ++jfield){
// for(int jphi=0; jphi < fields->element[jfield]->nlocal; ++jphi){
// for(int ifields=0; ifields < fields->n; ++ifields){
// for(int iphi=0; iphi < fields->element[ifields]->nlocal; ++iphi){
// local_matrix[]
// }
// }
// }
// }
}
void fe_fields_free(FEFields *fields);
......
python/csr.py
View file @ 9fe1e3db
......@@ -7,7 +7,7 @@ class CSR :
num = np.max(idx)
self.ndof = num+1
idx = idx.astype('uint64')
pairs = (idx[:,:,None]*shift+idx[:,None,:]).flatten()
pairs = (idx[:,None,:]*shift+idx[:,:,None]).flatten()
allpairs = [pairs]
self.constraints = constraints
self.idx = idx.reshape([-1]).astype('int64')
......@@ -32,30 +32,12 @@ class CSR :
self.size = self.row.size-1
def assemble_rhs(self, localv, u, constraints_value) :
localv_print = localv.reshape([4,-1])
ii = 0
for elt in range(4):
print("-------- elt : ------------ ", elt)
for i in range(9):
print("\t:", "%3d"%self.idx[ii], "->%2.4f" % localv_print[elt,i], end="")
ii += 1
print("")
rhs = np.bincount(self.idx,localv,self.size)
for i, (c, cv) in enumerate(zip(self.constraints, constraints_value)):
rhs[i+self.ndof] = cv[1] + np.sum(u[c]*cv[0])
return rhs
def assemble_mat(self, localm, constraints_value, v) :
localm_print = localm.reshape([4,9,-1])
ii = 0
for elt in range(4):
print("-------- elt : ------------ ", elt)
for i in range(9):
for j in range(9):
print("\t:", "%4d"%self.map[0][ii], "->%2.5f" % localm_print[elt,i,j], end="")
ii += 1
print("")
v[:] = np.bincount(self.map[0],localm,self.col.size)
for cmap, cv in zip (self.map[1:], constraints_value) :
v += np.bincount(cmap, np.concatenate((cv[0],cv[0],[0])), self.col.size)
......
python/fluid.py
View file @ 9fe1e3db
......@@ -157,11 +157,6 @@ class FluidProblem :
self._lib = lib3
else :
raise ValueError("dim should be 2 or 3.")
# if etypes is None and dim == 2:
# etypes = ["triangle_p1","triangle_p1", "triangle_p1"]
# if etypes is None and dim == 3:
# etypes = ["tetrahedron_p1","tetrahedron_p1","tetrahedron_p1","tetrahedron_p1"]
self._lib.fluid_problem_new.restype = c_void_p
n_fluids = np.require(rho,"float","C").reshape([-1]).shape[0]
self._n_fluids = n_fluids
......@@ -432,6 +427,10 @@ class FluidProblem :
solution_old = self.solution().copy()
self._lib.fluid_problem_reset_boundary_force(self._fp)
self._lib.fluid_problem_compute_stab_parameters(self._fp,c_double(dt))
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()
if self.sys is None :
constraints = []
......@@ -440,7 +439,8 @@ class FluidProblem :
# todo periodic count !!!
constraint = np.column_stack((np.arange(self.n_nodes()), np.repeat(self._dim, self.n_nodes())))
constraints.append(constraint)
self.sys = self.linear_system_package(periodic_map[self.elements()],self.n_fields(),self.solver_options, constraints)
self.sys = self.linear_system_package(idx, self.n_fields(), self.solver_options, constraints)
# self.sys = self.linear_system_package(periodic_map[self.elements()], self.n_fields(),self.solver_options, constraints)
sys = self.sys
localv = np.ndarray([sys.localsize*self.n_elements()])
localm = np.ndarray([sys.localsize**2*self.n_elements()])
......@@ -626,6 +626,13 @@ class FluidProblem :
def node_volume(self) :
"""Returns the volume associated with each node"""
return self._get_matrix("node_volume",self.n_nodes(),1)
def local_size(self):
self._lib.fluid_problem_local_size.restype = c_int
return self._lib.fluid_problem_local_size(self._fp)
def local_map(self, map):
self._lib.fluid_problem_local_map(self._fp, _np2c(map, dtype=np.int32))
def _n_physicals(self):
f = self._lib.fluid_problem_private_n_physical
......
python/petsc4pylsys.py
View file @ 9fe1e3db
......@@ -11,22 +11,25 @@ except :
print("dof_numbering not available !!!")
class LinearSystemAIJ:
def __init__(self, elements, n_fields, options, constraints) :
def __init__(self, idx, n_fields, options, constraints) :
# idx = ((elements*n_fields)[:,:,None]+np.arange(n_fields)[None,None,:]).reshape([elements.shape[0],-1])
self.localsize = elements.shape[1]*n_fields
# self.localsize = elements.shape[1]*n_fields
self.idx = idx
self.localsize = idx.shape[1]
self.constraints = list(c[:,0]*n_fields + c[:,1] for c in constraints)
# self.idx = (elements[:,None,:]*n_fields+np.arange(n_fields)[None,:,None]).reshape([-1])
self.space_v = LagrangeSpace(MeshCl(elements), order=1).element_nodes
self.space_p = LagrangeSpace(MeshCl(elements), order=1).element_nodes
dimension = n_fields-1
nmax_v = self.space_v.max()+1
self.localsize = dimension*self.space_v.shape[1] + self.space_p.shape[1]
self.globalsize = nmax_v*dimension + self.space_p.max()+1
# self.space_v = LagrangeSpace(MeshCl(elements), order=1).element_nodes
# self.space_p = LagrangeSpace(MeshCl(elements), order=1).element_nodes
# dimension = n_fields-1
# nmax_v = self.space_v.max()+1
# self.localsize = dimension*self.space_v.shape[1] + self.space_p.shape[1]
# self.globalsize = nmax_v*dimension + self.space_p.max()+1
# format u0,u1,...,v0,v1,...,p0,p1,...
self.idx = np.tile(self.space_v, dimension) \
+ (np.repeat(np.arange(dimension), self.space_v.shape[1])*nmax_v)[None,:]
pmap = dimension*nmax_v + self.space_p
self.idx = np.concatenate([self.idx, pmap], axis=1)
# self.idx = np.tile(self.space_v, dimension) \
# + (np.repeat(np.arange(dimension), self.space_v.shape[1])*nmax_v)[None,:]
# pmap = dimension*nmax_v + self.space_p
# self.idx = np.concatenate([self.idx, pmap], axis=1)
PETSc.Options().prefixPush("fluid_")
PETSc.Options().insertString("-pc_type ilu -pc_factor_levels 5 -pc_factor_shift_type INBLOCKS -pc_factor_shift_amount 1e-16")
PETSc.Options().insertString(options)
......@@ -44,6 +47,7 @@ class LinearSystemAIJ:
def local_to_global(self, localv, localm, u, constraints_value):
self.csr.assemble_mat(localm, constraints_value, self.val)
self.mat.assemble()
self.mat.view()
return self.csr.assemble_rhs(localv, u, constraints_value)
@timeit
......
validation/poiseuille/poiseuille.py
View file @ 9fe1e3db
......@@ -89,7 +89,7 @@ class Poiseuille(unittest.TestCase) :
tic = time.perf_counter()
while ii < 100 :
#Fluid solver
time_integration.iterate(fluid,None,dt,check_residual_norm=1e-6)
time_integration.iterate(fluid,None,dt)
t += dt
#Output files writting
if ii %outf == 0 :
......@@ -97,6 +97,7 @@ class Poiseuille(unittest.TestCase) :
fluid.write_vtk(outputdir, ioutput, t)
ii += 1
print("%i : %.2g/%.2g (cpu %.6g)" % (ii, t, tEnd, time.perf_counter() - tic))
exit(0)
s = fluid.solution()
x = fluid.coordinates()
......
Supports Markdown
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