Hi Netgen/NGSolve developers,
recently I switched from the old *.geo files to the python capabilities of netgen and ngsolve for meshing complex geometries.
However, I noticed, that importing Mesh from ngsolve.comp prevents other python libraries from executing code in parallel (especially the pymumps library).
An example code would look like the following:
import numpy as np
from scipy.sparse import coo_matrix,triu,identity
from ngsolve.comp import Mesh
import mumps
# GENERATE RANDOM SPARSE BANDED SYMMETRIC POSITIVE DEFINITE MATRIX
n = 50_000
n_bandwidth = int(n*0.005)
n_data = int(n**2*0.002)
row_ids = np.random.randint(0,high=n_bandwidth,size=n_data)+np.repeat(np.arange(n),int(n_data/n))
col_ids = np.random.randint(0,high=n_bandwidth,size=n_data)+np.repeat(np.arange(n),int(n_data/n))
tmp = np.zeros_like(col_ids)
tmp[::2] = n_bandwidth*20
col_ids+= tmp
row_ids = np.clip(row_ids,a_min=0,a_max=n-1)
col_ids = np.clip(col_ids,a_min=0,a_max=n-1)
K = coo_matrix( (np.random.rand(n_data), (row_ids,col_ids) ), shape=(n, n) )
K = K +K.T+ identity(n)
F = np.random.rand(n)
# SOLVE WITH MUMPS
ctx = mumps.DMumpsContext(sym=1, par=1) # set matrix to symmetric,
if ctx.myid == 0:
ctx.set_shape( K.shape[0] )
ctx.set_centralized_sparse( triu(K) )
u = F.copy().flatten()
ctx.set_rhs(u)
ctx.run(job=6) # Analysis + Factorization + Solve
ctx.destroy() # Free memory
print( 'residum: Ku -F = {}'.format(np.linalg.norm((K@u -F))) )
When ever I put the from ngsolve.comp import Mesh statment before import mumps no parallel execution is done.
Why is that? Is that a wanted behavior?