Commit 0659028d authored by Jonathan Lambrechts's avatar Jonathan Lambrechts
Browse files

constraints in linear solvers (+ mean pressures) +petsc4pylsys +no

tkinter in show problem

commit 9b1a78fc
Author: Jonathan Lambrechts <jonathan.lambrechts@uclouvain.be>
Date:   Tue Oct 12 11:34:12 2021 +0200

    minor polishing

commit 13c3b511
Author: Jonathan Lambrechts <jonathan.lambrechts@uclouvain.be>
Date:   Mon Oct 11 21:01:43 2021 +0200

    fix validation

commit 0127284b
Author: Jonathan Lambrechts <jonathan.lambrechts@uclouvain.be>
Date:   Mon Oct 11 20:45:34 2021 +0200

    missing files

commit 6aced34c
Author: Jonathan Lambrechts <jonathan.lambrechts@uclouvain.be>
Date:   Mon Oct 11 20:35:37 2021 +0200

    mumps elements with constraints

commit 0357264b
Author: Jonathan Lambrechts <jonathan.lambrechts@uclouvain.be>
Date:   Mon Oct 11 11:52:15 2021 +0200

    wip

commit 2c3e6d43
Author: Jonathan Lambrechts <jonathan.lambrechts@uclouvain.be>
Date:   Mon Oct 11 11:51:59 2021 +0200

    do not use tkinter in show_problem
parent 21676316
Pipeline #9837 passed with stages
in 7 minutes and 45 seconds
......@@ -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;
......
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include "mpi.h"
#include "dmumps_c.h"
#include "string.h"
......@@ -14,6 +15,57 @@ void mumps_set_options(DMUMPS_STRUC_C *id, int *options, int n){
}
}
DMUMPS_STRUC_C *mumps_new(int n, int nnz, int *row, int *col) {
DMUMPS_STRUC_C *id = malloc(sizeof(DMUMPS_STRUC_C));
int myid, ierr;
ierr = MPI_Init(NULL, NULL);
//ierr = MPI_Comm_rank(MPI_COMM_WORLD, &myid);
/*Initialize a MUMPS instance. Use MPI_COMM_WORLD.*/
id->job=JOB_INIT;
id->par=1;
id->sym=0;
id->comm_fortran=USE_COMM_WORLD;
dmumps_c(id);
/*Define the problem on the host*/
id->n = n;
id->nz = nnz;
id->irn = malloc(nnz*sizeof(int));
memcpy(id->irn, row, sizeof(int)*nnz);
id->jcn = malloc(nnz*sizeof(int));
memcpy(id->jcn, col, sizeof(int)*nnz);
/*No outputs*/
/*id->ICNTL(1)=-1;
id->ICNTL(2)=-1;
id->ICNTL(3)=-1;*/
/* output level */
id->ICNTL(4)=-1;
/* element format */
id->ICNTL(5)=0;
/* permutes the matrix to a zero-free diagonal */
id->ICNTL(6)=2;
/* reordering metis */
id->ICNTL(7)=5;
/* transpose matrix (local matrices c->f) */
id->ICNTL(9)=1;
/* no scaling */
id->ICNTL(8)=0;
/* When significant extra fill-in is caused by numerical pivoting, increasing ICNTL(14) may help */
/* 50 makes all the validation functional */
id->ICNTL(14)=50;
/* num of openmp threads */
id->ICNTL(16)=4;
return id;
}
void mumps_solve(DMUMPS_STRUC_C *id, double *a, double *rhs) {
int init = id->a == NULL;
id->a = a;
id->rhs = rhs;
id->job = (init ? 6 : 5);
dmumps_c(id);
}
DMUMPS_STRUC_C *mumps_element_new(int n, int nelt, int *eltptr, int *eltvar, double *a_elt, double *rhs) {
DMUMPS_STRUC_C *id = malloc(sizeof(DMUMPS_STRUC_C));
int myid, ierr;
......@@ -58,6 +110,7 @@ DMUMPS_STRUC_C *mumps_element_new(int n, int nelt, int *eltptr, int *eltvar, dou
id->ICNTL(16)=4;
return id;
}
void mumps_element_solve(DMUMPS_STRUC_C *id, double *a_elt, double *rhs) {
int init = id->a_elt == NULL;
id->a_elt = a_elt;
......
from os.path import isfile
import tkinter as tk
from tkinter.filedialog import askdirectory
# Create the box to select the output directory
boxroot = tk.Tk()
boxroot.withdraw()
dirname = askdirectory(parent=boxroot)
boxroot.destroy()
import os.path
from paraview.simple import *
try :
dirname = os.path.dirname(GetActiveSource().FileName)
except:
raise ValueError("Cannot get the directory of the active source.")
Delete()
# Show the fluid velocity
def showFluid(fluid, renderView, animationScene):
fluidDisplay = Show(fluid, renderView)
......@@ -65,7 +62,7 @@ def showContacts(particleProblem, renderView):
RenameSource('Contacts', contacts)
return contacts
except NameError:
print("MigFlow extract contacts Filter has not been charged !\nPlease see the wiki to add it in your Paraview filter.")
print("MigFlow extract contacts Filter has not been charged !\nPlease see the wiki to add it in your Paraview filter.")
return None
def showFluidVelocity(fluid, renderView):
......@@ -85,7 +82,7 @@ def showFluidVelocity(fluid, renderView):
animationScene = GetAnimationScene()
animationScene.PlayMode = 'Snap To TimeSteps'
renderView = GetActiveViewOrCreate('RenderView')
renderView.InteractionMode = '3D'
renderView.InteractionMode = '2D'
# Read the pvd files
particleProblemFile = str(dirname)+"/particle_problem.pvd"
......@@ -114,4 +111,4 @@ if isfile(fluidFile):
# Show the first time step
animationScene.GoToFirst()
renderView.Update()
renderView.ResetCamera()
\ No newline at end of file
renderView.ResetCamera()
......@@ -24,9 +24,11 @@ set(SRC
fluid.py
_tools.py
lmgc90Interface.py
csr.py
scontact.py
time_integration.py
petsclsys.py
petsc4pylsys.py
scipylsys.py
VTK.py
)
......
......@@ -37,6 +37,11 @@ try :
have_petsc = True
except :
have_petsc = False
try :
from . import petsc4pylsys
have_petsc4py = True
except :
have_petsc4py = False
try :
from . import scipylsys
have_scipy = True
......@@ -58,5 +63,8 @@ def get_linear_system_package(choices=None):
if have_scipy and choice=="scipy":
print("using scipy linear system")
return scipylsys.LinearSystem
if have_petsc4py and choice=="petsc4py":
print("using petsc4py linear system")
return petsc4pylsys.LinearSystem
raise ValueError("linear system not available")
import numpy as np
class CSR :
def __init__(self, idx, rhsidx,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,:]
pairs = pairs.reshape([-1])
allpairs = [pairs.reshape([-1])]
num = np.max(idx)
self.ndof = num+1
self.constraints = constraints
for c in constraints :
num += 1
pairs = np.ndarray([c.size*2+1],dtype=([('i0',np.int32),('i1',np.int32)]))
pairs['i0'][:c.size] = c
pairs['i1'][:c.size] = num
pairs['i0'][c.size:c.size*2] = num
pairs['i1'][c.size:c.size*2] = c
pairs['i0'][c.size*2] = num
pairs['i1'][c.size*2] = num
allpairs.append(pairs)
pairs = np.concatenate(allpairs)
pairs, pmap = np.unique(pairs,return_inverse=True)
self.map = []
count = 0
for p in allpairs :
self.map.append(pmap[count:count+p.size])
count += p.size
self.row = np.hstack([np.array([0],dtype=np.int32),
np.cumsum(np.bincount(pairs["i0"]),
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)
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
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) :
v += np.bincount(cmap, np.concatenate((cv[0],cv[0],[0])), self.col.size)
......@@ -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 :
......
......@@ -25,6 +25,8 @@ import numpy as np
import ctypes
import ctypes.util
import os
from ._tools import timeit
from .csr import CSR
dir_path = os.path.dirname(os.path.realpath(__file__))
lib2 = np.ctypeslib.load_library("libmbfluid2",dir_path)
......@@ -37,15 +39,62 @@ def _np2c(a,dtype=float,order="C") :
r.tmp = tmp
return r
class LinearSystem2 :
@timeit
def __init__(self, elements, n_fields, options=None, constraints=[]) :
self.localsize = elements.shape[1]*n_fields
self.constraints = list(c[:,0]*n_fields + c[:,1] for c in constraints)
idx = ((elements*n_fields)[:,:,None]+np.arange(n_fields)[None,None,:]).reshape([elements.shape[0],-1])
self.idx = (elements[:,None,:]*n_fields+np.arange(n_fields)[None,:,None]).reshape([-1])
self.csr = CSR(idx, self.idx,self.constraints)
self.val = np.zeros(self.csr.col.size)
row = np.repeat(np.arange(self.csr.size,dtype=np.int32),self.csr.row[1:]-self.csr.row[:-1])
lib2.mumps_new.restype = ctypes.c_void_p
self.mumps_ptr = ctypes.c_void_p(lib2.mumps_new(ctypes.c_int(self.csr.size),ctypes.c_int(self.csr.row[-1]),_np2c(row+1,np.int32),_np2c(self.csr.col+1,np.int32)))
options = None if options=="" else options
if options is not None:
options_vec = []
for key, val in options.items():
if isinstance(key,str):
key = options_key[key]
options_vec.extend([key,val])
options_vec = np.array(options_vec, dtype=np.int32)
lib2.mumps_set_options(self.mumps_ptr, _np2c(options_vec, np.int32), ctypes.c_int(np.size(options_vec)//2))
@timeit
def local_to_global(self,localv,localm,u, constraints_value):
self.csr.assemble_mat(localm, constraints_value, self.val)
return self.csr.assemble_rhs(localv, u, constraints_value)
@timeit
def solve(self,rhs) :
r = np.copy(rhs)
lib2.mumps_solve(self.mumps_ptr,_np2c(self.val),_np2c(r))
return r[:self.csr.ndof]
def __delete__(self):
lib2.mumps_element_delete(self.mumps_ptr)
class LinearSystem :
def __init__(self, elements, n_fields, options=None) :
@timeit
def __init__(self, elements, n_fields, options=None, constraints=[]) :
n_nodes = np.max(elements)+1
self.globalsize = n_nodes*n_fields
self.globalsize = n_nodes*n_fields+len(constraints)
self.localsize = elements.shape[1]*n_fields
self.idx = (elements[:,None,:]*n_fields+np.arange(n_fields)[None,:,None]).reshape([-1])
idxf = ((elements[:,:,None]*n_fields+np.arange(n_fields)[None,None,:]).reshape([-1])+1).astype(np.int32)
eltptr = np.arange(1,elements.size*n_fields+2,self.localsize)
self.idx = (elements[:,None,:]*n_fields+np.arange(n_fields)[None,:,None])
idxf = [self.idx.swapaxes(1,2)+1]
eltptr = [np.arange(1,elements.size*n_fields+2,self.localsize)]
self.ndof = n_nodes*n_fields
self.n_fields = n_fields
self.constraints = constraints
for i,constraint in enumerate(constraints):
idxc = np.c_[np.full(constraint.shape[0], self.ndof+i+1), constraint[:,0]*n_fields+constraint[:,1]+1].flatten()
idxf.append(idxc)
eltptr.append(eltptr[-1][-1]+np.arange(2,constraint.shape[0]*2+2,2))
lib2.mumps_element_new.restype = ctypes.c_void_p
eltptr = np.concatenate([i.flat for i in eltptr])
idxf = np.concatenate([i.flat for i in idxf])
self.mumps_ptr = ctypes.c_void_p(lib2.mumps_element_new(ctypes.c_int(self.globalsize),ctypes.c_int(eltptr.size-1),_np2c(eltptr,np.int32),_np2c(idxf,np.int32)))
options = None if options=="" else options
if options is not None:
......@@ -54,19 +103,30 @@ class LinearSystem :
if isinstance(key,str):
key = options_key[key]
options_vec.extend([key,val])
print(options_vec)
options_vec = np.array(options_vec, dtype=np.int32)
lib2.mumps_set_options(self.mumps_ptr, _np2c(options_vec, np.int32), ctypes.c_int(np.size(options_vec)//2))
@timeit
def local_to_global(self,localv,localm,u, constraints_values):
rhs = np.bincount(self.idx.flatten(),localv,self.globalsize)
for i, (c, cv) in enumerate(zip(self.constraints,constraints_values)):
rhs[i+self.ndof] = cv[1] + np.sum(u[c[:,0]*self.n_fields+c[:,1]]*cv[0])
self.v = np.empty(localm.size+sum(c.shape[0]*4 for c in self.constraints))
p = localm.size
self.v[:p] = localm.flat
for constraint, cv in constraints_values:
localc = np.zeros((constraint.shape[0],4))
localc[:,1] = constraint
localc[:,2] = constraint
self.v[p:localc.size+p] = localc.flat
p += localc
return rhs
def local_to_global(self,localv,localm,rhs):
rhs[:] = np.bincount(self.idx,localv,self.globalsize)
self.v = localm.reshape([-1])
@timeit
def solve(self,rhs) :
r = np.copy(rhs)
lib2.mumps_element_solve(self.mumps_ptr,_np2c(self.v),_np2c(r))
return r
return r[:self.ndof]
def __delete__(self):
lib2.mumps_element_delete(self.mumps_ptr)
......
import petsc4py
import sys
petsc4py.init(sys.argv)
from petsc4py import PETSc
import numpy as np
from ._tools import timeit
from .csr import CSR
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])
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])
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)
PETSc.Options().prefixPop()
self.ksp = PETSc.KSP().create()
self.ksp.setOptionsPrefix(b"fluid_")
self.ksp.setFromOptions()
self.csr = CSR(idx, 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)
self.mat.assemble()
return self.csr.assemble_rhs(localv, u, constraints_value)
@timeit
def solve(self,rhs) :
x = np.ndarray(rhs.shape)
prhs = PETSc.Vec().createWithArray(rhs.reshape([-1]))