Commit 2793a8a9 authored by Jonathan Lambrechts's avatar Jonathan Lambrechts
Browse files

isolate linear system python wrapper

parent b70d2fd6
Pipeline #7670 passed with stage
in 2 minutes and 3 seconds
......@@ -62,6 +62,39 @@ def _is_string(s) :
else :
return isinstance(s, basestring)
def _np2c(a,dtype=float) :
tmp = np.require(a,dtype,"C")
r = c_void_p(tmp.ctypes.data)
r.tmp = tmp
return r
class LinearSystem :
def __init__(self, elements, n_fields) :
self._lib = lib2
self._lib.fluid_problem_create_linear_system.restype = c_void_p
self._sysp = c_void_p(self._lib.fluid_problem_create_linear_system(c_int(elements.shape[0]),c_int(elements.shape[1]),_np2c(elements,np.int32),c_int(n_fields)))
self.localsize = n_fields*elements.shape[1]
self._elements = elements
self._lib.fluid_problem_zero_matrix(self._sysp)
self._lib.fluid_problem_get_rhs_norm.restype = c_double
def zero_matrix(self) :
self._lib.fluid_problem_zero_matrix(self._sysp)
def local_to_global(self,localv,localm,rhs):
self._lib.fluid_problem_local_to_global(self._sysp,c_int(self._elements.shape[0]),c_int(self.localsize),_np2c(localv),_np2c(localm),_np2c(rhs))
def rhs_norm(self,rhs) :
return self._lib.fluid_problem_get_rhs_norm(self._sysp,_np2c(rhs));
def solve(self,rhs) :
sol = np.ndarray(rhs.shape)
self._lib.fluid_problem_linear_system_solve(self._sysp,_np2c(rhs),_np2c(sol));
return sol
def __del__(self):
self._lib.fluid_problem_free_linear_system(self._sysp)
class FluidProblem :
"""Create numerical representation of the fluid."""
......@@ -88,11 +121,6 @@ class FluidProblem :
"""
self.strong_cb_ref = []
self.weak_cb_ref = []
def np2c(a) :
tmp = np.require(a,"float","C")
r = c_void_p(tmp.ctypes.data)
r.tmp = tmp
return r
if dim == 2 :
self._lib = lib2
elif dim == 3 :
......@@ -103,7 +131,7 @@ 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), petsc_solver_type.encode(),c_int(drag_in_stab),c_double(drag_coefficient_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), petsc_solver_type.encode(),c_int(drag_in_stab),c_double(drag_coefficient_factor),c_int(temporal), c_int(advection)))
if self._fp == None :
raise NameError("Cannot create fluid problem.")
......@@ -220,12 +248,7 @@ class FluidProblem :
cmax -- optional argument to weigh maximum gradient used in the adaptation formula
cmin -- optional argument to weigh minimum gradient used in the adaptation formula
"""
def np2c(a) :
tmp = np.require(a,"float","C")
r = c_void_p(tmp.ctypes.data)
r.tmp = tmp
return r
self._lib.fluid_problem_adapt_mesh(self._fp, c_double(lcmax), c_double(lcmin), c_double(n_el), c_double(cmax), c_double(cmin), c_int(old_n_particles), np2c(old_particle_position), np2c(old_particle_volume))
self._lib.fluid_problem_adapt_mesh(self._fp, c_double(lcmax), c_double(lcmin), c_double(n_el), c_double(cmax), c_double(cmin), c_int(old_n_particles), _np2c(old_particle_position), _np2c(old_particle_volume))
def _mesh_boundaries(self) :
n = self._lib.fluid_problem_mesh_n_boundaries(self._fp)
......@@ -277,11 +300,6 @@ class FluidProblem :
Keyword arguments:
filename -- name of the file to read
"""
def np2c(a,dtype) :
tmp = np.require(a,dtype,"C")
r = c_void_p(tmp.ctypes.data)
r.tmp = tmp
return r
x,cells,data,cdata,fdata = VTK.read(filename)
mesh_boundaries = {name[9:]:nodes for name,nodes in fdata.items() if name.startswith("Boundary ")}
cbnd_names = (c_char_p*len(mesh_boundaries))(*(i.encode() for i in mesh_boundaries.keys()))
......@@ -291,8 +309,8 @@ class FluidProblem :
bnd_tags = np.repeat(list(range(nbnd)),list([v.shape[0] for v in mesh_boundaries.values()]))
bnd_tags = np.require(bnd_tags,np.int32,"C")
self._lib.fluid_problem_set_elements(self._fp,
c_int(x.shape[0]),np2c(x,np.float64),
c_int(el.shape[0]),np2c(el,np.int32),
c_int(x.shape[0]),_np2c(x,np.float64),
c_int(el.shape[0]),_np2c(el,np.int32),
c_int(bnds.shape[0]),c_void_p(bnds.ctypes.data),c_void_p(bnd_tags.ctypes.data),c_int(len(cbnd_names)),cbnd_names
)
sol = self.solution()
......@@ -325,6 +343,7 @@ class FluidProblem :
self._lib.fluid_problem_compute_node_particle_torque(self._fp, c_double(dt), c_void_p(torques.ctypes.data))
return torques
def implicit_euler(self, dt, check_residual_norm=-1, reduced_gravity=0, stab_param=0.) :
"""Solve the fluid equations.
......@@ -334,43 +353,30 @@ 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
"""
def np2c(a,dtype=float) :
tmp = np.require(a,dtype,"C")
r = c_void_p(tmp.ctypes.data)
r.tmp = tmp
return r
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()
rhs = np.zeros_like(solution_old)
self._lib.fluid_problem_reset_boundary_force(self._fp)
dx = np.ndarray(solution_old.shape)
self._lib.fluid_problem_compute_stab_parameters(self._fp,c_double(dt))
self._lib.fluid_problem_create_linear_system.restype = c_void_p
sysp = c_void_p(self._lib.fluid_problem_create_linear_system(c_int(self.n_elements()),c_int(self.dimension()+1),np2c(self.elements(),np.int32),c_int(self.n_fields())))
self._lib.fluid_problem_zero_matrix(sysp)
localsize = self.n_fields()*(self.dimension()+1)
localv = np.ndarray([localsize*self.n_elements()])
localm = np.ndarray([localsize*localsize*self.n_elements()])
self._lib.fluid_problem_assemble_system(self._fp,np2c(rhs),np2c(solution_old),c_double(dt),np2c(localv),np2c(localm))
self._lib.fluid_problem_local_to_global(sysp,c_int(self.n_elements()),c_int(localsize),np2c(localv),np2c(localm),np2c(rhs))
self._lib.fluid_problem_get_rhs_norm.restype = c_double
sys = LinearSystem(self.elements(),self.n_fields())
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(rhs),_np2c(solution_old),c_double(dt),_np2c(localv),_np2c(localm))
sys.local_to_global(localv,localm,rhs)
if check_residual_norm > 0 :
norm0 = self._lib.fluid_problem_get_rhs_norm(sysp,np2c(rhs));
self._lib.fluid_problem_linear_system_solve(sysp,np2c(rhs),np2c(dx));
self.solution()[:] -= dx
norm0 = sys.rhs_norm(rhs)
self.solution()[:] -= sys.solve(rhs)
if check_residual_norm > 0 :
self._lib.fluid_problem_zero_matrix(sysp)
sys.zero_matrix()
rhs[:] = 0
self._lib.fluid_problem_assemble_system(self._fp,np2c(rhs),np2c(solution_old),c_double(dt),np2c(localv),np2c(localm))
self._lib.fluid_problem_local_to_global(sysp,c_int(self.n_elements()),c_int(localsize),np2c(localv),np2c(localm),np2c(rhs))
norm = self._lib.fluid_problem_get_rhs_norm(sysp,np2c(rhs));
self._lib.fluid_problem_assemble_system(self._fp,_np2c(rhs),_np2c(solution_old),c_double(dt),_np2c(localv),_np2c(localm))
sys.local_to_global(localv,localm,rhs)
norm = sys.rhs_norm(rhs)
if norm > check_residual_norm :
raise ValueError("wrong derivative or linear system precision")
self._lib.fluid_problem_free_linear_system(sysp)
self._lib.fluid_problem_node_force_volume(self._fp,np2c(solution_old),c_double(dt),None,None)
self._lib.fluid_problem_node_force_volume(self._fp,_np2c(solution_old),c_double(dt),None,None)
def compute_cfl(self):
"""Compute the CFL number divided by the time step
......@@ -407,13 +413,8 @@ class FluidProblem :
init -- if True set the old_porosity to zero (has to be used at the first call)
reload -- optional arguments specifying if the simulation is reloaded or if it is the first iteration
"""
def np2c(a) :
tmp = np.require(a,"float","C")
r = c_void_p(tmp.ctypes.data)
r.tmp = tmp
return r
self.n_particles = mass.shape[0]
self._lib.fluid_problem_set_particles(self._fp,c_int(mass.shape[0]),np2c(mass),np2c(volume),np2c(position),np2c(velocity),np2c(contact),None,c_int(init),c_int(reload))
self._lib.fluid_problem_set_particles(self._fp,c_int(mass.shape[0]),_np2c(mass),_np2c(volume),_np2c(position),_np2c(velocity),_np2c(contact),None,c_int(init),c_int(reload))
def _get_matrix(self, f_name, nrow, ncol,typec=c_double) :
......@@ -469,13 +470,8 @@ class FluidProblem :
def set_concentration_cg(self,concentration):
"""Set the concentration at nodes"""
def np2c(a) :
tmp = np.require(a,"float","C")
r = c_void_p(tmp.ctypes.data)
r.tmp = tmp
return r
concentration = concentration.reshape((self.n_nodes(),1))
self._lib.fluid_problem_set_concentration_cg(self._fp,np2c(concentration))
self._lib.fluid_problem_set_concentration_cg(self._fp,_np2c(concentration))
def concentration_dg(self):
"""Return the concentration at discontinuous nodes"""
......
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