Commit 0357264b authored by Jonathan Lambrechts's avatar Jonathan Lambrechts
Browse files

wip

parent 2c3e6d43
......@@ -1140,7 +1140,7 @@ static void fluid_problem_dg_to_cg_grad(FluidProblem *problem, const double *dg,
}
}
void fluid_problem_assemble_system(FluidProblem *problem, double *rhs, const double *solution_old, double dt, double *all_local_vector, double *all_local_matrix)
void fluid_problem_assemble_system(FluidProblem *problem, const double *solution_old, double dt, double *all_local_vector, double *all_local_matrix)
{
const Mesh *mesh = problem->mesh;
......
......@@ -67,63 +67,63 @@ def _np2c(a,dtype=float) :
return r
def _load_msh(mesh_file_name, lib, dim) :
"""
mesh_file_name -- Name of the mesh.msh file containing information about the domain
"""
etype = 2 if dim == 2 else 4
betype = 1 if dim == 2 else 2
gmsh.model.add("tmp")
gmsh.open(mesh_file_name)
gmsh.model.mesh.renumberNodes()
ntags, x, _ = gmsh.model.mesh.getNodes()
x = x.reshape([-1,3])
x[ntags-1] = x
_, el = gmsh.model.mesh.getElementsByType(etype)
el = (el-1).reshape([-1,dim+1])
bnd = []
btag = []
bname = []
periodic_nodes = []
for edim, etag in gmsh.model.getEntities() :
ptag, cnodes, pnodes, _ = gmsh.model.mesh.getPeriodicNodes(edim, etag)
if ptag == etag or len(cnodes) == 0: continue
periodic_nodes.extend(zip(cnodes,pnodes))
for ip,(pdim, ptag) in enumerate(gmsh.model.getPhysicalGroups(dim-1)) :
bname.append(gmsh.model.getPhysicalName(pdim,ptag))
for etag in gmsh.model.getEntitiesForPhysicalGroup(pdim, ptag):
for eltype, _, elnodes in zip(*gmsh.model.mesh.getElements(pdim,etag)):
if eltype != betype:
continue
bnd.append(elnodes-1)
btag.append(np.full((elnodes.shape[0]//dim), ip))
periodic_nodes = np.array(periodic_nodes,dtype=np.int32).reshape([-1,2])-1
bnd = np.hstack(bnd).reshape([-1,dim])
btag = np.hstack(btag)
cbname = (c_char_p*len(bname))(*(i.encode() for i in bname))
# remove nodes and boundaries not connected to elements
keepnodes = np.full(x.shape[0], False, bool)
keepnodes[el] = True
newidx = np.cumsum(keepnodes)-1
el = newidx[el]
x = x[keepnodes,:]
periodic_nodes = newidx[periodic_nodes]
keepbnd = np.full(bnd.shape[0], True, bool)
for i in range(dim) :
keepbnd = np.logical_and(keepbnd, keepnodes[bnd[:,i]])
bnd = newidx[bnd][keepbnd,:]
btag = btag[keepbnd]
is_parent = np.full([x.shape[0]],True,bool)
is_parent[periodic_nodes[:,0]] = False
periodic_parent = np.cumsum(is_parent)-1
periodic_parent[periodic_nodes[:,0]] = periodic_parent[periodic_nodes[:,1]]
lib.mesh_new_from_elements.restype = c_void_p
return c_void_p(lib.mesh_new_from_elements(
c_int(x.shape[0]),_np2c(x,np.float64),
c_int(el.shape[0]),_np2c(el,np.int32),
c_int(bnd.shape[0]),_np2c(bnd,np.int32),
_np2c(btag,np.int32),c_int(len(cbname)),cbname,
_np2c(periodic_parent,np.int32)))
"""
mesh_file_name -- Name of the mesh.msh file containing information about the domain
"""
etype = 2 if dim == 2 else 4
betype = 1 if dim == 2 else 2
gmsh.model.add("tmp")
gmsh.open(mesh_file_name)
gmsh.model.mesh.renumberNodes()
ntags, x, _ = gmsh.model.mesh.getNodes()
x = x.reshape([-1,3])
x[ntags-1] = x
_, el = gmsh.model.mesh.getElementsByType(etype)
el = (el-1).reshape([-1,dim+1])
bnd = []
btag = []
bname = []
periodic_nodes = []
for edim, etag in gmsh.model.getEntities() :
ptag, cnodes, pnodes, _ = gmsh.model.mesh.getPeriodicNodes(edim, etag)
if ptag == etag or len(cnodes) == 0: continue
periodic_nodes.extend(zip(cnodes,pnodes))
for ip,(pdim, ptag) in enumerate(gmsh.model.getPhysicalGroups(dim-1)) :
bname.append(gmsh.model.getPhysicalName(pdim,ptag))
for etag in gmsh.model.getEntitiesForPhysicalGroup(pdim, ptag):
for eltype, _, elnodes in zip(*gmsh.model.mesh.getElements(pdim,etag)):
if eltype != betype:
continue
bnd.append(elnodes-1)
btag.append(np.full((elnodes.shape[0]//dim), ip))
periodic_nodes = np.array(periodic_nodes,dtype=np.int32).reshape([-1,2])-1
bnd = np.hstack(bnd).reshape([-1,dim])
btag = np.hstack(btag)
cbname = (c_char_p*len(bname))(*(i.encode() for i in bname))
# remove nodes and boundaries not connected to elements
keepnodes = np.full(x.shape[0], False, bool)
keepnodes[el] = True
newidx = np.cumsum(keepnodes)-1
el = newidx[el]
x = x[keepnodes,:]
periodic_nodes = newidx[periodic_nodes]
keepbnd = np.full(bnd.shape[0], True, bool)
for i in range(dim) :
keepbnd = np.logical_and(keepbnd, keepnodes[bnd[:,i]])
bnd = newidx[bnd][keepbnd,:]
btag = btag[keepbnd]
is_parent = np.full([x.shape[0]],True, bool)
is_parent[periodic_nodes[:,0]] = False
periodic_parent = np.cumsum(is_parent)-1
periodic_parent[periodic_nodes[:,0]] = periodic_parent[periodic_nodes[:,1]]
lib.mesh_new_from_elements.restype = c_void_p
return c_void_p(lib.mesh_new_from_elements(
c_int(x.shape[0]),_np2c(x,np.float64),
c_int(el.shape[0]),_np2c(el,np.int32),
c_int(bnd.shape[0]),_np2c(bnd,np.int32),
_np2c(btag,np.int32),c_int(len(cbname)),cbname,
_np2c(periodic_parent,np.int32)))
class FluidProblem :
......@@ -156,6 +156,7 @@ class FluidProblem :
self.strong_cb_ref = []
self.weak_cb_ref = []
self.sys = None
self.mean_pressure = None
if dim == 2 :
self._lib = lib2
elif dim == 3 :
......@@ -189,6 +190,10 @@ class FluidProblem :
self.sys = None
gmsh.model.remove()
def set_mean_pressure(self, mean_pressure):
self.mean_pressure = mean_pressure
self.sys = None
def set_wall_boundary(self, tag, pressure=None, velocity=None, compute_viscous_term=1) :
"""Sets the weak boundaries (=normal fluxes) for the fluid problem.
......@@ -360,10 +365,6 @@ class FluidProblem :
offsets = np.cumsum(np.repeat([self._dim+1],connectivity.shape[0])).astype(np.int32)
VTK.write(output_dir+"/fluid",it,t,(types,connectivity,offsets),self.coordinates(),data,field_data,cell_data)
def export_vtk(self, output_dir, t, it, stab=False) :
print("Careful this function is deprecated... \n\tUse write_vtk instead !")
self.write_vtk(output_dir, it, t, stab)
def read_vtk(self, dirname, i):
"""Reads output file to reload computation.
......@@ -396,34 +397,6 @@ class FluidProblem :
self.old_porosity()[:] = data["old_porosity"]
self.sys = None
def import_vtk(self, filename) :
print("Careful this function is deprecated... \n\tUse read_vtk instead !")
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()))
el = cells["connectivity"].reshape([-1,self._dim+1])
nbnd = len(mesh_boundaries)
bnds = np.vstack(list(mesh_boundaries.values()))
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.mesh_new_from_elements.restype = c_void_p
_mesh = c_void_p(self._lib.mesh_new_from_elements(
c_int(x.shape[0]),_np2c(x,np.float64),
c_int(el.shape[0]),_np2c(el,np.int32),
c_int(bnds.shape[0]),_np2c(bnds,np.int32),
_np2c(bnd_tags,np.int32),c_int(len(cbnd_names)),cbnd_names,
_np2c(data["parent_node_id"],np.int32) if "parent_node_id" in data else None))
self._lib.fluid_problem_set_mesh(self._fp, _mesh)
sol = self.solution()
sol[:,:self._dim] = data["velocity"][:,:self._dim]
sol[:,[self._dim]] = data["pressure"]
if self._n_fluids == 2 :
self.concentration_dg()[:] = cdata["concentration"]
self.porosity()[:] = data["porosity"]
self.old_porosity()[:] = data["old_porosity"]
self.sys = None
def compute_node_force(self, dt) :
"""Computes the forces to apply on each grain resulting from the fluid motion.
......@@ -463,18 +436,25 @@ class FluidProblem :
periodic_map = self._get_matrix("periodic_mapping", self.n_nodes(), 1, c_int).flatten()
if self.sys is None :
self.sys = self.linear_system_package(periodic_map[self.elements()],self.n_fields(),self.solver_options)
constraints = []
if self.mean_pressure is not None:
node_volume = self.node_volume()
# 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)
sys = self.sys
rhs = np.zeros(sys.globalsize)
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)
self.solution()[:] -= sys.solve(rhs).reshape([-1,self.n_fields()])[periodic_map]
self._lib.fluid_problem_assemble_system(self._fp,_np2c(solution_old),c_double(dt),_np2c(localv),_np2c(localm))
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]
if check_residual_norm > 0 :
rhs[:] = 0
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)
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)
norm = np.linalg.norm(rhs)
print("norm",norm)
if norm > check_residual_norm :
......
......@@ -69,12 +69,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)
fluid = mbfluid.FluidProblem(2,g,nu*rho,rho,solver="petsc")
fluid.set_mean_pressure(0)
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",pressure=0,velocity=[vUp,0])
fluid.set_wall_boundary("Top",velocity=[vUp,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