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

wip

parent 6b85afeb
Pipeline #9900 failed with stages
in 2 minutes and 45 seconds
......@@ -56,6 +56,9 @@ try :
except :
have_scipy = False
# only use petsc4py for p2p1 now
have_mumps, have_petsc, have_scipy = False, False, False
def get_linear_system_package(choices=None):
if choices is None:
choices = ["mumps","petsc","petsc4py","scipy"]
......
import numpy as np
import time
from ._tools import timeit
class CSR :
@timeit
def __init__(self, rhsidx,constraints):
shift = 2**32
num = np.max(rhsidx)
......@@ -31,12 +32,14 @@ class CSR :
self.size = self.row.size-1
self.rhsidx = rhsidx.reshape([-1])
@timeit
def assemble_rhs(self, localv, u, constraints_value) :
rhs = np.bincount(self.rhsidx,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
@timeit
def assemble_mat(self, localm, constraints_value, v) :
v[:] = np.bincount(self.map[0],localm,self.col.size)
for cmap, cv in zip (self.map[1:], constraints_value) :
......
......@@ -29,11 +29,12 @@
"""
from ctypes import *
from time import time
import numpy as np
import os
import sys
from . import VTK
from ._tools import gmsh, get_linear_system_package, _np2c
from ._tools import gmsh, get_linear_system_package, _np2c, timeit
dir_path = os.path.dirname(os.path.realpath(__file__))
lib2 = np.ctypeslib.load_library("libmbfluid2",dir_path)
......@@ -133,7 +134,7 @@ class FluidProblem :
density -- List of fluid phases densities (this should be a vector whom dimension specify if there is one or two fluid)
sigma -- Surface tension (only when there are two fluids)
coeff_stab -- Optional argument used as coefficient for extra diffusivity added to stabilise the advection of concentration (only for two fluid problem)
solver -- possible solver are "mumps", "petsc", "scipy"
solver -- possible solver are "mumps", "petsc", "scipy", "petsc4py"
solver_options -- Optional argument to specify solver option (only when petsc is used)
drag_in_stab -- States if the drag force is in the stabilisation term
drag_coefficient_factor -- Factor multiplicating the drag coefficient
......@@ -412,6 +413,7 @@ class FluidProblem :
self._lib.fluid_problem_compute_node_particle_torque(self._fp, c_double(dt), c_void_p(torques.ctypes.data))
return torques
@timeit
def implicit_euler(self, dt, check_residual_norm=-1, reduced_gravity=0, stab_param=0.) :
"""Solves the fluid equations.
......@@ -421,13 +423,14 @@ class FluidProblem :
reduced_graviy -- If True simulations are run with a reduced gravity (not to use in two fluids simulations)
stab_param -- If zero pspg/supg/lsic stabilisation is computed otherwise the value is used as a coefficient for a pressure laplacian stabilisation term
"""
# timer = time()
self._lib.fluid_problem_set_reduced_gravity(self._fp,c_int(reduced_gravity))
self._lib.fluid_problem_set_stab_param(self._fp,c_double(stab_param))
self._lib.fluid_problem_apply_boundary_conditions(self._fp)
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))
# 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])
......@@ -445,11 +448,13 @@ class FluidProblem :
localv = np.ndarray([sys.localsize*self.n_elements()])
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)
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]
# print(time()-timer)
if check_residual_norm > 0 :
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)
......@@ -458,6 +463,8 @@ class FluidProblem :
if norm > check_residual_norm :
raise ValueError("wrong derivative or linear system precision")
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
......@@ -631,6 +638,7 @@ class FluidProblem :
self._lib.fluid_problem_local_size.restype = c_int
return self._lib.fluid_problem_local_size(self._fp)
@timeit
def local_map(self, map):
self._lib.fluid_problem_local_map(self._fp, _np2c(map, dtype=np.int32))
......
......@@ -23,7 +23,6 @@ class LinearSystemAIJ:
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)
......
......@@ -41,7 +41,6 @@ lib.PetscInitialize(None, None, None)
COMM_WORLD = ctypes.c_void_p.in_dll(lib,"PETSC_COMM_WORLD")
#atexit.register(lib.PetscFinalize)
class PetscObj:
def __init__(self, constructor, destructor, *args):
......@@ -104,11 +103,11 @@ class LinearSystem:
lib.MatSeqBAIJSetPreallocationCSR(self.mat,
ctypes.c_int(n_fields), _np2c(self.csr.row, np.int32), _np2c(self.csr.col, np.int32), None)
else:
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.csr = CSR(idx, self.idx,self.constraints)
self.csr = CSR(self.idx,self.constraints)
self.val = np.zeros(self.csr.col.size)
self.mat = PetscObj("MatCreateSeqAIJWithArrays", "MatDestroy", COMM_WORLD,
ctypes.c_int(self.csr.size), ctypes.c_int(self.csr.size),
......
......@@ -69,12 +69,12 @@ 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="petsc4py")
# fluid.set_mean_pressure(-100000)
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],pressure=0)
fluid.set_wall_boundary("Top",velocity=[vUp,0])
fluid.set_wall_boundary("Right",velocity=[0,0])
ii = 0
t = 0
......
......@@ -64,7 +64,7 @@ class Poiseuille(unittest.TestCase) :
outf = 1 # number of iterations between output files
#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,solver="petsc4py")
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])
......
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