Commit bec26ebb authored by Michel Henry's avatar Michel Henry
Browse files

wip

parent 6ed38c0f
Pipeline #9843 failed with stages
in 1 minute and 21 seconds
......@@ -29,6 +29,7 @@ set(SRC
petsclsys.py
scipylsys.py
VTK.py
dof_numbering.py
)
foreach(f ${SRC})
......
......@@ -143,8 +143,7 @@ def test(dim, order):
gmsh.view.add_list_data(view, "SP", xdof.shape[0], data.flat)
gmsh.merge("view_dof.opt")
if __name__ == "__main__":
test(2,2)
gmsh.fltk.run()
gmsh.finalize()
# test(2,2)
# gmsh.fltk.run()
# gmsh.finalize()
......@@ -34,6 +34,7 @@ import os
import sys
from . import VTK
from ._tools import gmsh,timeit
try :
from .petsclsys import LinearSystem
except :
......@@ -467,10 +468,9 @@ 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))
localsize = self._lib.fluid_problem_local_size(self._fp)
periodic_map = self._get_matrix("periodic_mapping", self.n_nodes(), 1, c_int).flatten()
if self.sys is None :
self.sys = LinearSystem(periodic_map[self.elements()],self.n_fields(),self.solver_options, localsize)
self.sys = LinearSystem(periodic_map[self.elements()],self.n_fields(),self.solver_options)
sys = self.sys
rhs = np.zeros(sys.globalsize)
localv = np.ndarray([sys.localsize*self.n_elements()])
......
......@@ -21,41 +21,68 @@
import petsc4py
petsc4py.init()
from petsc4py import PETSc
from .dof_numbering import MeshCl, LagrangeSpace
import numpy as np
try :
from .dof_numbering import MeshCl, LagrangeSpace
except :
print("dof_numbering not available !!!")
class LinearSystem :
def __init__(self, elements, n_fields, options, localsize) :
def __init__(self, elements, n_fields, options) :
nnodes = np.max(elements)+1
nn = elements.shape[1]
pairs = np.ndarray((elements.shape[0],nn*(nn-1)),dtype=([('i0',np.int32),('i1',np.int32)]))
pairs['i0'] = elements[:,np.repeat(range(nn),nn-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.elements = elements
# format u0,u1,...,v0,v1,...,p0,p1,...
dofmap = np.tile(self.space_v, dimension)
dofmap = dofmap \
+ (np.repeat(np.arange(dimension), self.space_v.shape[1])*nmax_v)[None,:]
pmap = dimension*nmax_v + self.space_p
dofmap = np.concatenate([dofmap, pmap], axis=1)
# nn = elements.shape[1]
# pairs = np.ndarray((elements.shape[0],nn*(nn-1)),dtype=([('i0',np.int32),('i1',np.int32)]))
# pairs['i0'] = elements[:,np.repeat(range(nn),nn-1)]
# idx = np.hstack([[i for i in range(nn) if i!=j] for j in range(nn)])
# pairs['i1'] = elements[:,idx]
# pairs = np.unique(pairs)
# nnz = (np.bincount(pairs["i0"])+1).astype(np.int32)
nn = self.space_v.shape[1]
# nn = dofmap.shape[1]
pairs = np.ndarray((elements.shape[0],nn*(nn-1)), dtype=([('i0', np.int32), ('i1', np.int32)]))
pairs['i0'] = dofmap[:, np.repeat(range(nn), nn-1)]
idx = np.hstack([[i for i in range(nn) if i!=j] for j in range(nn)])
pairs['i1'] = elements[:,idx]
pairs['i1'] = dofmap[:,idx]
pairs = np.unique(pairs)
nnz = (np.bincount(pairs["i0"])+1).astype(np.int32)
PETSc.Options().prefixPush("fluid_")
PETSc.Options().insertString(options)
PETSc.Options().prefixPop()
self.mat = PETSc.Mat().createBAIJ(nnodes*n_fields,n_fields,nnz)
# self.mat = PETSc.Mat().createBAIJ(self.globalsize, n_fields, nnz)
self.ksp = PETSc.KSP().create()
self.ksp.setOptionsPrefix(b"fluid_")
self.ksp.setFromOptions()
self.mat.zeroEntries()
###
self.localsize = localsize
self.globalsize = n_fields*nnodes
self.elements = elements
# self.idx = dofmap.reshape([-1])
self.idx = (self.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
self.idx = np.asarray(self.idx, np.int32)
def local_to_global(self,localv,localm,rhs):
self.mat.zeroEntries()
np.add.at(rhs.reshape([-1]),self.idx,localv)
n_fields = self.localsize//self.elements.shape[1]
# print(self.space_v)
vmap = np.tile(self.elements*n_fields,[1,n_fields]) + np.repeat(np.arange(n_fields,dtype=np.int32),n_fields)[None,:]
# n_fields = self.localsize//self.elements.shape[1]
# vmap = np.tile(self.elements*n_fields,[1,n_fields]) + np.repeat(np.arange(n_fields,dtype=np.int32),n_fields)[None,:]
vmap = self.idx.reshape([self.elements.shape[0],-1])
np.add.at(rhs.reshape([-1]),self.idx,localv) # ok
for e,m in zip(vmap,localm.reshape([-1,self.localsize**2])) :
self.mat.setValues(e,e,m,PETSc.InsertMode.ADD)
#for e,m in zip(self.elements,localm.reshape([-1,self.localsize**2])) :
......
......@@ -24,8 +24,9 @@ import scipy.sparse.linalg
class LinearSystem :
_initialized = False
def __init__(self, elements, n_fields,options,localsize) :
def __init__(self, elements, n_fields,options) :
self.el = elements
localsize = n_fields*elements.shape[1]
idx = (elements[:,:,None]*n_fields+np.arange(n_fields)[None,None,:]).reshape([-1,localsize])
self.rows = np.repeat(idx[:,:,None],localsize,axis=2)
self.rows = self.rows.reshape([-1])
......
......@@ -65,12 +65,15 @@ class Poiseuille(unittest.TestCase) :
#Object fluid creation + Boundary condition of the fluid (field 0 is horizontal velocity; field 1 is vertical velocity; field 2 is pressure)
fluid = mbfluid.FluidProblem(2,g,nu*rho,rho)
fluid.load_msh("mesh.msh")
fluid.set_open_boundary("LeftUp",velocity=[lambda x : 1/(20*mu)*x[:,1]*(1-x[:, 1]),0])
fluid.set_open_boundary("LeftDown",velocity=[lambda x : 1/(20*mu)*x[:,1]*(1-x[:, 1]),0])
fluid.set_open_boundary("Left",velocity=[lambda x : 1/(20*mu)*x[:,1]*(1-x[:, 1]),0])
fluid.set_open_boundary("Right",pressure=0)
# fluid.set_open_boundary("LeftUp",velocity=[lambda x : 1/(20*mu)*x[:,1]*(1-x[:, 1]),0])
# fluid.set_open_boundary("LeftDown",velocity=[lambda x : 1/(20*mu)*x[:,1]*(1-x[:, 1]),0])
fluid.set_wall_boundary("Bottom",velocity=[0,0])
fluid.set_wall_boundary("Top",velocity=[0,0])
fluid.set_open_boundary("RightUp",pressure=0)
fluid.set_open_boundary("RightDown",pressure=0)
# fluid.set_open_boundary("RightUp",pressure=0)
# fluid.set_open_boundary("RightDown",pressure=0)
ii = 0
......
Markdown is supported
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