Commit 60391df2 authored by Michel Henry's avatar Michel Henry
Browse files

wip

parent 3a0417b7
......@@ -26,6 +26,7 @@ def timeprint(timers) :
print(timers)
atexit.register(timeprint, timers)
from . import petsc4pylsys
try :
from . import mumpslsys
......
import numpy as np
class CSR :
def __init__(self, idx, rhsidx,constraints):
def __init__(self, idx, constraints):
pairs = np.ndarray([idx.shape[0],idx.shape[1],idx.shape[1]],dtype=([('i0',np.int32),('i1',np.int32)]))
pairs['i0'][:,:,:] = idx[:,:,None]
pairs['i1'][:,:,:] = idx[:,None,:]
......@@ -10,6 +10,7 @@ class CSR :
num = np.max(idx)
self.ndof = num+1
self.constraints = constraints
self.idx = idx.reshape([-1])
for c in constraints :
num += 1
pairs = np.ndarray([c.size*2+1],dtype=([('i0',np.int32),('i1',np.int32)]))
......@@ -32,10 +33,9 @@ class CSR :
dtype=np.int32)])
self.col = pairs['i1'].copy()
self.size = self.row.size-1
self.rhsidx = rhsidx
def assemble_rhs(self, localv, u, constraints_value) :
rhs = np.bincount(self.rhsidx,localv,self.size)
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
......
......@@ -5,14 +5,28 @@ from petsc4py import PETSc
import numpy as np
from ._tools import timeit
from .csr import CSR
try :
from .dof_numbering import MeshCl, LagrangeSpace
except :
print("dof_numbering not available !!!")
class LinearSystemAIJ:
def __init__(self, elements, n_fields, options, constraints) :
idx = ((elements*n_fields)[:,:,None]+np.arange(n_fields)[None,None,:]).reshape([elements.shape[0],-1])
# idx = ((elements*n_fields)[:,:,None]+np.arange(n_fields)[None,None,:]).reshape([elements.shape[0],-1])
self.localsize = elements.shape[1]*n_fields
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.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
# 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)
PETSc.Options().prefixPush("fluid_")
PETSc.Options().insertString("-pc_type lu -pc_factor_shift_type NONZERO -pc_factor_shift_amount 1e-16")
PETSc.Options().insertString(options)
......@@ -20,10 +34,11 @@ class LinearSystemAIJ:
self.ksp = PETSc.KSP().create()
self.ksp.setOptionsPrefix(b"fluid_")
self.ksp.setFromOptions()
self.csr = CSR(idx, self.idx,self.constraints)
self.csr = CSR(self.idx,self.constraints)
self.val = np.zeros(self.csr.col.size)
self.mat = PETSc.Mat().createAIJWithArrays(self.csr.size,(self.csr.row,self.csr.col,self.val),bsize=1)
@timeit
def local_to_global(self, localv, localm, u, constraints_value):
self.csr.assemble_mat(localm, constraints_value, self.val)
......
......@@ -64,7 +64,7 @@ r = 8e-3 # particles radius
rin = 0.01905 # central disk radius
rhop = 2500 # particles density
# Numerical parameters
outf = 1000 # number of iterations between output files
outf = 50 # number of iterations between output files
dt = 1e-3 # time step
tEnd = 10 # final time
# Geometrical parameters
......
......@@ -68,13 +68,13 @@ 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)
#Format: strong_boundaries = [(Boundary tag, Fluid field, Value)
fluid = mbfluid.FluidProblem(2,g,nu*rho,rho,solver="mumps")
fluid.set_mean_pressure(-100000)
fluid = mbfluid.FluidProblem(2,g,nu*rho,rho,solver="petsc4py")
# fluid.set_mean_pressure(-100000)
fluid.load_msh("mesh.msh")
fluid.set_wall_boundary("Left",velocity=[0,0])
fluid.set_wall_boundary("Bottom",velocity=[0,0])
fluid.set_wall_boundary("Top",velocity=[vUp,0])
fluid.set_wall_boundary("Top",velocity=[vUp,0],pressure=0)
fluid.set_wall_boundary("Right",velocity=[0,0])
ii = 0
t = 0
......
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