Algebraic Multigrid for Compressed VectorH1 space

Ryanleaf · February 19, 2021, 11:39pm

Hello guys,

I met a problem when I was trying to assemble algebraic multigrid preconditioner from a compressed VectorH1 space. Here is the script to show my problem,.

from netgen.geom2d import SplineGeometry from netgen.meshing import MeshingParameters from xfem import * from netgen.csg import CSGeometry, OrthoBrick, Pnt import netgen.gui from netgen.geom2d import unit_square from ngsolve import * from ngsolve.la import EigenValues_Preconditioner from xfem.lsetcurv import * def Setup(): cube = CSGeometry() cube.Add(OrthoBrick(Pnt(-1.5, -1.5, -1.5), Pnt(1.5, 1.5, 1.5))) mesh = Mesh(cube.GenerateMesh(maxh=1, quad_dominated=False)) lsetmeshadap = LevelSetMeshAdaptation(mesh, order=3, threshold=1000, discontinuous_qn=True) phi = Norm(CoefficientFunction((x, y, z))) - 1 deformation = lsetmeshadap.CalcDeformation(phi) lsetp1 = lsetmeshadap.lset_p1 mesh.SetDeformation(deformation) ci = CutInfo(mesh, lsetp1) Vhbase = VectorH1(mesh, order=2, flags={"dgjumps": True}) fes = Compress(Vhbase, GetDofsOfElements(Vhbase, ci.GetElementsOfType(HASNEG))) lset_neg = {"levelset": lsetp1, "domain_type": NEG} lset_pos = {"levelset": lsetp1, "domain_type": POS} lset_doms = [lset_neg, lset_pos] u, v = fes.TnT() a = BilinearForm(fes) a += SymbolicBFI(lset_doms[0], form= u * v) return mesh, fes, a mesh, fes, a = Setup() mg = Preconditioner(a, "h1amg") # register mg to a a.Assemble() a.inner_matrix print('NDOF = ', fes.ndof)
When I run this part of code I immediately got the following segmentation fault,

[code]importing ngs-xfem1.3-2008
optfile ./ng.opt does not exist - using default values
togl-version : 2
OCC module loaded
loading ngsolve library
NGSolve-6.2.2008
Using Lapack
Including sparse direct solver Pardiso
Running parallel using 8 thread(s)
Start Findpoints
Analyze spec points
Find edges
Start Findpoints
Analyze spec points
Find edges
Start Findpoints
Analyze spec points
Find edges
Surface 1 / 6
Optimize Surface
Surface 2 / 6
Optimize Surface
Surface 3 / 6
Optimize Surface
Surface 4 / 6
Optimize Surface
Surface 5 / 6
Optimize Surface
Surface 6 / 6
Optimize Surface

Meshing subdomain 1 of 1
Delaunay meshing
start tetmeshing
Success !
72 points, 211 elements
Remove Illegal Elements
Volume Optimization
SearchCorrespondingPoint:: did not converge
WARNING: using flags as kwarg is deprecated in <class ‘ngsolve.comp.VectorH1’>, use the flag arguments as kwargs instead!
[Qis-MacBook-Pro:10639] *** Process received signal ***
[Qis-MacBook-Pro:10639] Signal: Segmentation fault: 11 (11)
[Qis-MacBook-Pro:10639] Signal code: Address not mapped (1)
[Qis-MacBook-Pro:10639] Failing at address: 0x7fbcee1e99cc
[Qis-MacBook-Pro:10639] [ 0] 0 libsystem_platform.dylib 0x00007fff6ecf05fd _sigtramp + 29
[Qis-MacBook-Pro:10639] [ 1] 0 ??? 0x0000000a0000000c 0x0 + 42949672972
[Qis-MacBook-Pro:10639] [ 2] 0 libngcore.dylib 0x00000001102f5726 _ZN6ngcore11TaskManager9CreateJobERKNSt3__18functionIFvRNS_8TaskInfoEEEEi + 198
[Qis-MacBook-Pro:10639] [ 3] 0 libngcomp.dylib 0x0000000115b5b1d1 _ZN6ngcomp12H1AMG_MatrixIdEC1ENSt3__110shared_ptrIN4ngla14SparseMatrixTMIdEEEENS3_IN6ngcore8BitArrayEEENS8_9FlatArrayINS8_3INTILi2EiEEmEENSB_IdmEESF_m + 1681
[Qis-MacBook-Pro:10639] [ 4] 0 libngcomp.dylib 0x0000000115b66c8c ZNSt3__110shared_ptrIN6ngcomp12H1AMG_MatrixIdEEE11make_sharedIJRNS0_IN4ngla14SparseMatrixTMIdEEEERNS0_IN6ngcore8BitArrayEEERNSB_5ArrayINSB_3INTILi2EiEEmEERNSF_IdmEESL_iEEES4_DpOT + 236
[Qis-MacBook-Pro:10639] [ 5] 0 libngcomp.dylib 0x0000000115b65540 _ZN6ngcomp20H1AMG_PreconditionerIdE13FinalizeLevelEPKN4ngla10BaseMatrixE + 1120
[Qis-MacBook-Pro:10639] [ 6] 0 libngcomp.dylib 0x0000000115730de4 _ZN6ngcomp14S_BilinearFormIdE10DoAssembleERN6ngcore9LocalHeapE + 12180
[Qis-MacBook-Pro:10639] [ 7] 0 libngcomp.dylib 0x000000011574d044 _ZN6ngcomp12BilinearForm8AssembleERN6ngcore9LocalHeapE + 692
[Qis-MacBook-Pro:10639] [ 8] 0 libngcomp.dylib 0x0000000115ea7504 _ZZN8pybind1112cpp_function10initializeIZ12ExportNgcompRNS_6moduleEE5$142NSt3__110shared_ptrIN6ngcomp12BilinearFormEEEJS9_bEJNS_4nameENS_9is_methodENS_7siblingENS_10call_guardIJNS_18gil_scoped_releaseEEEENS_5arg_vEPKcEEEvOT_PFT0_DpT1_EDpRKT2_ENUlRNS_6detail13function_callEE_8__invokeESW + 148
[Qis-MacBook-Pro:10639] [ 9] 0 libngcomp.dylib 0x0000000115a9f4b9 ZN8pybind1112cpp_function10dispatcherEP7_objectS2_S2 + 4041
[Qis-MacBook-Pro:10639] [10] 0 python 0x000000010f0cb1b8 _PyMethodDef_RawFastCallKeywords + 392
[Qis-MacBook-Pro:10639] [11] 0 python 0x000000010f0cac80 _PyObject_FastCallKeywords + 592
[Qis-MacBook-Pro:10639] [12] 0 python 0x000000010f207ed5 call_function + 453
[Qis-MacBook-Pro:10639] [13] 0 python 0x000000010f205aec _PyEval_EvalFrameDefault + 46092
[Qis-MacBook-Pro:10639] [14] 0 python 0x000000010f1f949e _PyEval_EvalCodeWithName + 414
[Qis-MacBook-Pro:10639] [15] 0 python 0x000000010f25c9a0 PyRun_FileExFlags + 256
[Qis-MacBook-Pro:10639] [16] 0 python 0x000000010f25be17 PyRun_SimpleFileExFlags + 391
[Qis-MacBook-Pro:10639] [17] 0 python 0x000000010f289d3f pymain_main + 9663
[Qis-MacBook-Pro:10639] [18] 0 python 0x000000010f09d66d main + 125
[Qis-MacBook-Pro:10639] [19] 0 libdyld.dylib 0x00007fff6eaf7cc9 start + 1
[Qis-MacBook-Pro:10639] *** End of error message ***
H1AMG: level = 0, num_edges = 19827, nv = 615[/code]
I found when I removed the
" fes = Compress(Vhbase, GetDofsOfElements(Vhbase, ci.GetElementsOfType(HASNEG)))"
And if I use VhBase as FEM space for BilinearForm a and assemble the corresponding algebraic multigrid, everything works. I am wondering if I need to add something extra to make algebraic multigrid work for this compressed VectorH1 space?
Thanks you!

Qi

schruste · February 20, 2021, 11:14am

Dear Qi,

Even with fes = Vhbase it’s not working. The problem is that there are elements where the element matrix is zero and while trying to find proper weights for the H1amg local element matrix inversions fail (the standard use case for h1amg does not assume to work with empty contributions here). What you have to do is to make sure that you don’t compute element matrices for the bilinear form a (and hence for the preconditioner) in the first place on in-active elements. To do so replace your integrator line with

a += SymbolicBFI(lset_doms[0], form= u * v, definedonelements=ci.GetElementsOfType(HASNEG))

Then only on active elements element matrices are computed for the matrix a and the h1amg preconditioner.

Best,
Christoph

Ryanleaf · February 23, 2021, 3:56am

Thank you very much Christoph. This really solved my problem.
Sorry, One thing I just want to make sure.

[quote]Vhbase = VectorH1(mesh, order=order, dirichlet=“bottom|right|back”, flags={“dgjumps”: True})
Vhneg = Compress(Vhbase, GetDofsOfElements(Vhbase, ci.GetElementsOfType(HASNEG)))
Vhpos = Compress(Vhbase, GetDofsOfElements(Vhbase, ci.GetElementsOfType(HASPOS)))
WhGV = FESpace([Vhneg, Vhpos], flags={“dgjumps”: True})[/quote]
Suppose I have a two-phase problem, Can I use Algebraic Multigrid for the matrix from space WhGV?
Or I need to build a block preconditioner for the matrix?
Again, thank you very much.

schruste · February 23, 2021, 7:23am

Dear Qi,

Algorithmically it should work. If it is performing well also depends on the diffusivity contrast, cut position (you are not using stabilization of any kind, right?). An example where it is doing well using h1amg is given in the py_tutorials/mpi directory.

Best,
Christoph

Ryanleaf · February 23, 2021, 3:44pm

Dear Christoph,

Thank you very much for your help, I was asking because when I tried to use Algebraic multigrid for the space WhGV, again I am having segmentation fault.

[code]from netgen.geom2d import SplineGeometry
from netgen.meshing import MeshingParameters
from xfem import *
from netgen.csg import CSGeometry, OrthoBrick, Pnt
import netgen.gui
from netgen.geom2d import unit_square
from ngsolve import *
from ngsolve.la import EigenValues_Preconditioner
from xfem.lsetcurv import *

def eps(u):
return 0.5 * (Grad(u) + Grad(u).trans)
def Setup():
cube = CSGeometry()
cube.Add(OrthoBrick(Pnt(-1.5, -1.5, -1.5), Pnt(1.5, 1.5, 1.5)))
mesh = Mesh(cube.GenerateMesh(maxh=1, quad_dominated=False))
lsetmeshadap = LevelSetMeshAdaptation(mesh, order=3, threshold=1000, discontinuous_qn=True)
phi = Norm(CoefficientFunction((x, y, z))) - 1
deformation = lsetmeshadap.CalcDeformation(phi)
lsetp1 = lsetmeshadap.lset_p1
mesh.SetDeformation(deformation)
ci = CutInfo(mesh, lsetp1)
Vhbase = VectorH1(mesh, order=2, dirichlet=[1,2,3,4,5,6])
Vhneg = Compress(Vhbase, GetDofsOfElements(Vhbase, ci.GetElementsOfType(HASNEG)))
Vhpos = Compress(Vhbase, GetDofsOfElements(Vhbase, ci.GetElementsOfType(HASPOS)))
WhGV = FESpace([Vhneg, Vhpos], flags={“dgjumps”: True})
coef_normal = CoefficientFunction((x, y, z))
lset_neg = {“levelset”: lsetp1, “domain_type”: NEG}
lset_pos = {“levelset”: lsetp1, “domain_type”: POS}
lset_doms = [lset_neg, lset_pos]
elements_doms = [HASNEG, HASPOS]
u, v = WhGV.TnT()
a = BilinearForm(WhGV)
for i in [0, 1]: # loop over domains
a += SymbolicBFI(lset_doms[i], form=2 * InnerProduct(eps(u[i]), eps(v[i])),
definedonelements=ci.GetElementsOfType(elements_doms[i]))
a += SymbolicBFI(lset_doms[i], form=u[i] * v[i], definedonelements=ci.GetElementsOfType(elements_doms[i]))
return mesh, WhGV, a
mesh, fes, a = Setup()
mg = Preconditioner(a, “h1amg”) # register mg to a
a.Assemble()
print('NDOF = ', fes.ndof)
preI = Projector(mask=fes.FreeDofs(), range=True)
lams = EigenValues_Preconditioner(mat=a.mat, pre=mg.mat)
print(“smallest eign value”, max(lams)/min(lams))
[/code]

[quote]H1AMG: level = 0, num_edges = 102939, [Qis-MacBook-Pro:34374] *** Process received [Qis-MacBook-Pro:34374] Signal: Segmentation [Qis-MacBook-Pro:34374] Signal code: (0)
[Qis-MacBook-Pro:34374] Failing at address: [Qis-MacBook-Pro:34374] [ 0] 0 libsystem_platform.dylib [Qis-MacBook-Pro:34374] [ 1] 0 libsystem_c.dylib [Qis-MacBook-Pro:34374] [ 2] 0 libngcomp.dylib [Qis-MacBook-Pro:34374] [ 3] 0 libngcore.dylib [Qis-MacBook-Pro:34374] [ 4] 0 libngcomp.dylib [Qis-MacBook-Pro:34374] [ 5] 0 libngcomp.dylib [Qis-MacBook-Pro:34374] [ 6] 0 libngcomp.dylib [Qis-MacBook-Pro:34374] [ 7] 0 libngcomp.dylib [Qis-MacBook-Pro:34374] [ 8] 0 libngcomp.dylib [Qis-MacBook-Pro:34374] [ 9] 0 libngcomp.dylib [Qis-MacBook-Pro:34374] [10] 0 libngcomp.dylib [Qis-MacBook-Pro:34374] [11] 0 python [Qis-MacBook-Pro:34374] [12] 0 python [Qis-MacBook-Pro:34374] [13] 0 python [Qis-MacBook-Pro:34374] [14] 0 python [Qis-MacBook-Pro:34374] [15] 0 python [Qis-MacBook-Pro:34374] [16] 0 python [Qis-MacBook-Pro:34374] [17] 0 python [Qis-MacBook-Pro:34374] [18] 0 python [Qis-MacBook-Pro:34374] [19] 0 python [Qis-MacBook-Pro:34374] [20] 0 libdyld.dylib [Qis-MacBook-Pro:34374] *** End of error Is this somehow the same problem with the nv = 1692
signal ***
fault: 11 (11)
0x0
0x00007fff6ecf05fd _sigtramp + 29
0x00007fff95242d80 __c_locale + 0
0x000000010784016a ZNSt3__110__function6__funcIZN6ngcore11ParallelForImZN6ngcomp12H1AMG_MatrixINS_7complexIdEEEC1ENS_10shared_ptrIN4ngla14SparseMatrixTMIS7_EEEENS9_INS2_8BitArrayEEENS2_9FlatArrayINS2_3INTILi2EiEEmEENSG_IdmEESK_mEUliE1_EEvNS2_7T_RangeIT_EET0_iNS2_10TotalCostsEEUlRNS2_8TaskInfoEE_NS_9allocatorIST_EEFvSS_EEclESS + 90
0x00000001020cd726 _ZN6ngcore11TaskManager9CreateJobERKNSt3__18functionIFvRNS_8TaskInfoEEEEi + 198
0x000000010782ed14 _ZN6ngcomp12H1AMG_MatrixIdEC1ENSt3__110shared_ptrIN4ngla14SparseMatrixTMIdEEEENS3_IN6ngcore8BitArrayEEENS8_9FlatArrayINS8_3INTILi2EiEEmEENSB_IdmEESF_m + 4564
0x0000000107839c8c ZNSt3__110shared_ptrIN6ngcomp12H1AMG_MatrixIdEEE11make_sharedIJRNS0_IN4ngla14SparseMatrixTMIdEEEERNS0_IN6ngcore8BitArrayEEERNSB_5ArrayINSB_3INTILi2EiEEmEERNSF_IdmEESL_iEEES4_DpOT + 236
0x0000000107838540 _ZN6ngcomp20H1AMG_PreconditionerIdE13FinalizeLevelEPKN4ngla10BaseMatrixE + 1120
0x0000000107403de4 _ZN6ngcomp14S_BilinearFormIdE10DoAssembleERN6ngcore9LocalHeapE + 12180
0x0000000107420044 _ZN6ngcomp12BilinearForm8AssembleERN6ngcore9LocalHeapE + 692
0x0000000107b7a504 _ZZN8pybind1112cpp_function10initializeIZ12ExportNgcompRNS_6moduleEE5$142NSt3__110shared_ptrIN6ngcomp12BilinearFormEEEJS9_bEJNS_4nameENS_9is_methodENS_7siblingENS_10call_guardIJNS_18gil_scoped_releaseEEEENS_5arg_vEPKcEEEvOT_PFT0_DpT1_EDpRKT2_ENUlRNS_6detail13function_callEE_8__invokeESW + 148
0x00000001077724b9 ZN8pybind1112cpp_function10dispatcherEP7_objectS2_S2 + 4041
0x0000000100ea31b8 _PyMethodDef_RawFastCallKeywords + 392
0x0000000100ea2c80 _PyObject_FastCallKeywords + 592
0x0000000100fdfed5 call_function + 453
0x0000000100fddaec _PyEval_EvalFrameDefault + 46092
0x0000000100fd149e _PyEval_EvalCodeWithName + 414
0x00000001010349a0 PyRun_FileExFlags + 256
0x0000000101033e17 PyRun_SimpleFileExFlags + 391
0x0000000101061d3f pymain_main + 9663
0x0000000100e7566d main + 125
0x00007fff6eaf7cc9 start + 1
message ***[/quote]
previous one?

Best,
Qi

schruste · February 23, 2021, 10:34pm

Dear Qi,

Ok, my last reply was a bit premature. In the demo-file that I refered to we use a XFEM representation which harmonizes better with the h1amg in this case. The problem with the version that you tried is that on some elements there are some active dofs which however have no contribution to some local element matrices (e.g. on an element neighboring to cut element). This is in contrast to the XFEM formulation where dofs only occur on elements where they also have a contribution. This is not a flaw by the CutFEM formulation but rather of the implementation via the Compress-space. I see three possible solutions:

switch to an XFEM formulation
add a very small regularization which is independent of the cut configuration, e.g.

a += SymbolicBFI(form=1e-12 * u[i] * v[i])

Implement a version of the CompressedFESpace that also removes the link from elements to dofs where desired.

Best,
Christoph

Ryanleaf · February 24, 2021, 3:52am

Dear Christoph,

Thank you very much for your kind help.
I’ll discuss these possible solutions with my advisors.

Best,
Qi

schruste · February 24, 2021, 8:49am

OK, I just did the third approach. The new variant that comes with ngsxfem (on github, current master branch) now has a RestrictedFESpace that you can obtain by replacing “Compress” with “Restrict”. Instead of a BitArray of marked dofs you provide a BitArray of marked elements. The resulting dofs will be the same. However, when asked for element degrees of freedoms the RestrictedFESpace knows where no dofs should be placed at all.
With the change from

Vhneg = Compress(Vhbase, GetDofsOfElements(Vhbase, ci.GetElementsOfType(HASNEG)))
Vhpos = Compress(Vhbase, GetDofsOfElements(Vhbase, ci.GetElementsOfType(HASPOS)))

to

Vhneg = Restrict(Vhbase, ci.GetElementsOfType(HASNEG))
Vhpos = Restrict(Vhbase, ci.GetElementsOfType(HASPOS))

your example should go through again.

Best,
Christoph

Ryanleaf · February 25, 2021, 12:20am

Dear Christoph,

Thank you very much for your kind input, I will try to working on this.

Best,

Qi

Ryanleaf · March 3, 2021, 5:24pm

Dear Christoph,

I met another segmentation fault when I tried to add Dirichlet boundary conditions to this Restrict space. Here is part of script to demonstrate my problem.

[code]from netgen.geom2d import SplineGeometry
from netgen.meshing import MeshingParameters
from xfem import *
from netgen.csg import CSGeometry, OrthoBrick, Pnt
import netgen.gui
from netgen.geom2d import unit_square
from ngsolve import *
from ngsolve.la import EigenValues_Preconditioner
from xfem.lsetcurv import *

def eps(u):
return 0.5 * (Grad(u) + Grad(u).trans)
def Setup():
cube = CSGeometry()
cube.Add(OrthoBrick(Pnt(-1.5, -1.5, -1.5), Pnt(1.5, 1.5, 1.5)))
mesh = Mesh(cube.GenerateMesh(maxh=1, quad_dominated=False))
lsetmeshadap = LevelSetMeshAdaptation(mesh, order=3, threshold=1000, discontinuous_qn=True)
phi = Norm(CoefficientFunction((x, y, z))) - 1
deformation = lsetmeshadap.CalcDeformation(phi)
lsetp1 = lsetmeshadap.lset_p1
mesh.SetDeformation(deformation)
ci = CutInfo(mesh, lsetp1)
Vhbase = VectorH1(mesh, order=2, dirichlet=[1,2,3,4,5,6])
Vhneg = Restrict(Vhbase, ci.GetElementsOfType(HASNEG))
Vhpos = Restrict(Vhbase, ci.GetElementsOfType(HASPOS))
WhGV = FESpace([Vhneg, Vhpos], flags={“dgjumps”: True})
coef_normal = CoefficientFunction((x, y, z))
lset_neg = {“levelset”: lsetp1, “domain_type”: NEG}
lset_pos = {“levelset”: lsetp1, “domain_type”: POS}
lset_doms = [lset_neg, lset_pos]
elements_doms = [HASNEG, HASPOS]
u, v = WhGV.TnT()
a = BilinearForm(WhGV)
for i in [0, 1]: # loop over domains
a += SymbolicBFI(lset_doms[i], form=2* InnerProduct(eps(u[i]), eps(v[i])),
definedonelements=ci.GetElementsOfType(elements_doms[i]))
a += SymbolicBFI(lset_doms[i], form=u[i] * v[i], definedonelements=ci.GetElementsOfType(elements_doms[i]))
return mesh, WhGV, a
mesh, WhGV, a = Setup()
gfu = GridFunction(WhGV)
gfu.components[1].Set(CoefficientFunction((x,y,z)),BND)
mg = Preconditioner(a, “h1amg”) # register mg to a
a.Assemble()
print('NDOF = ', WhGV.ndof)
preI = Projector(mask=WhGV.FreeDofs(), range=True)
lams = EigenValues_Preconditioner(mat=a.mat, pre=mg.mat)
print(“condition number”, max(lams)/min(lams))[/code]

[quote]
Remove Illegal Elements
Volume Optimization
SearchCorrespondingPoint:: did not converge
WARNING: using flags as kwarg is deprecated in <class ‘ngsolve.comp.FESpace’>, use the flag arguments as kwargs instead!
[qi-Alienware-17-R4:08254] *** Process received signal ***
[qi-Alienware-17-R4:08254] Signal: Segmentation fault (11)
[qi-Alienware-17-R4:08254] Signal code: Address not mapped (1)
[qi-Alienware-17-R4:08254] Failing at address: 0x3e0000003f
[qi-Alienware-17-R4:08254] [ 0] /lib/x86_64-linux-gnu/libc.so.6(+0x46210)[0x7fbef6ab3210]
[qi-Alienware-17-R4:08254] [ 1] /home/qi/ngsuite/ngsolve-install/lib/python3/dist-packages/ngsolve/…/…/…/libngfem.so(_ZNK5ngfem22T_DifferentialOperatorINS_16DiffOpIdVectorH1ILi3ELNS_4VorBE1EEEE10CalcMatrixERKNS_13FiniteElementERKNS_30SIMD_BaseMappedIntegrationRuleEN5ngbla15BareSliceMatrixIN6ngcore4SIMDIdLi4EEELNSB_8ORDERINGE1EEE+0x11e)[0x7fbee8e2886e]
[qi-Alienware-17-R4:08254] [ 2] /home/qi/ngsuite/ngsolve-install/lib/python3/dist-packages/ngsolve/…/…/…/libngfem.so(_ZNK5ngfem30SymbolicBilinearFormIntegrator22T_CalcElementMatrixAddIdddEEvRKNS_13FiniteElementERKNS_21ElementTransformationEN5ngbla10FlatMatrixIT1_LNS8_8ORDERINGE1EmEERbRN6ngcore9LocalHeapE+0x8d6)[0x7fbee9b5d7a6]
[qi-Alienware-17-R4:08254] [ 3] /home/qi/ngsuite/ngsolve-install/lib/python3/dist-packages/ngsolve/…/…/…/libngfem.so(_ZNK5ngfem30SymbolicBilinearFormIntegrator17CalcElementMatrixERKNS_13FiniteElementERKNS_21ElementTransformationEN5ngbla10FlatMatrixIdLNS7_8ORDERINGE1EmEERN6ngcore9LocalHeapE+0x8b)[0x7fbee9b3e50b]
[qi-Alienware-17-R4:08254] [ 4] /home/qi/ngsuite/ngsolve-install/lib/python3/dist-packages/ngsolve/…/…/…/libngcomp.so(ZZN6ngcomp9SetValuesIdEEvSt10shared_ptrIN5ngfem19CoefficientFunctionEERNS_12GridFunctionENS2_4VorBEPKNS_6RegionEPNS2_20DifferentialOperatorERN6ngcore9LocalHeapEbbiENKUlNS_7FESpace7ElementESF_E0_clESH_SF+0x13f6)[0x7fbeeb2e23c6]
[qi-Alienware-17-R4:08254] [ 5] /home/qi/ngsuite/ngsolve-install/lib/python3/dist-packages/ngsolve/…/…/…/libngcomp.so(ZNSt17_Function_handlerIFvN6ngcomp7FESpace7ElementERN6ngcore9LocalHeapEEZNS0_9SetValuesIdEEvSt10shared_ptrIN5ngfem19CoefficientFunctionEERNS0_12GridFunctionENS9_4VorBEPKNS0_6RegionEPNS9_20DifferentialOperatorES5_bbiEUlS2_S5_E0_E9_M_invokeERKSt9_Any_dataOS2_S5+0x5a)[0x7fbeeb2e277a]
[qi-Alienware-17-R4:08254] [ 6] /home/qi/ngsuite/ngsolve-install/lib/python3/dist-packages/ngsolve/…/…/…/libngcomp.so(+0x750a6b)[0x7fbeeb042a6b]
[qi-Alienware-17-R4:08254] [ 7] /home/qi/ngsuite/ngsolve-install/lib/python3/dist-packages/netgen/…/…/…/libngcore.so(_ZN6ngcore11TaskManager9CreateJobERKSt8functionIFvRNS_8TaskInfoEEEi+0x10d)[0x7fbef458e95d]
[qi-Alienware-17-R4:08254] [ 8] /home/qi/ngsuite/ngsolve-install/lib/python3/dist-packages/ngsolve/…/…/…/libngcomp.so(_ZN6ngcomp15IterateElementsERKNS_7FESpaceEN5ngfem4VorBERN6ngcore9LocalHeapERKSt8functionIFvNS0_7ElementES7_EE+0x322)[0x7fbeeb042742]
[qi-Alienware-17-R4:08254] [ 9] /home/qi/ngsuite/ngsolve-install/lib/python3/dist-packages/ngsolve/…/…/…/libngcomp.so(_ZN6ngcomp9SetValuesIdEEvSt10shared_ptrIN5ngfem19CoefficientFunctionEERNS_12GridFunctionENS2_4VorBEPKNS_6RegionEPNS2_20DifferentialOperatorERN6ngcore9LocalHeapEbbi+0xdf4)[0x7fbeeb2d96e4]
[qi-Alienware-17-R4:08254] [10] /home/qi/ngsuite/ngsolve-install/lib/python3/dist-packages/ngsolve/…/…/…/libngcomp.so(_ZN6ngcomp9SetValuesESt10shared_ptrIN5ngfem19CoefficientFunctionEERNS_12GridFunctionENS1_4VorBEPNS1_20DifferentialOperatorERN6ngcore9LocalHeapEbbi+0x224)[0x7fbeeb2b6724]
[qi-Alienware-17-R4:08254] [11] /home/qi/ngsuite/ngsolve-install/lib/python3/dist-packages/ngsolve/…/…/…/libngcomp.so(+0xd21527)[0x7fbeeb613527]
[qi-Alienware-17-R4:08254] [12] /home/qi/ngsuite/ngsolve-install/lib/python3/dist-packages/ngsolve/…/…/…/libngcomp.so(+0x9b2aba)[0x7fbeeb2a4aba]
[qi-Alienware-17-R4:08254] [13] /usr/bin/python3.8[0x6b224d]
[qi-Alienware-17-R4:08254] [14] /usr/bin/python3.8(PyObject_Call+0x27e)[0x5f2c0e]
[qi-Alienware-17-R4:08254] [15] /usr/bin/python3.8(_PyEval_EvalFrameDefault+0x1f70)[0x56b7b0]
[qi-Alienware-17-R4:08254] [16] /usr/bin/python3.8(_PyEval_EvalCodeWithName+0x26a)[0x56822a]
[qi-Alienware-17-R4:08254] [17] /usr/bin/python3.8(_PyFunction_Vectorcall+0x393)[0x5f6033]
[qi-Alienware-17-R4:08254] [18] /usr/bin/python3.8(_PyEval_EvalFrameDefault+0x8f6)[0x56a136]
[qi-Alienware-17-R4:08254] [19] /usr/bin/python3.8(_PyEval_EvalCodeWithName+0x26a)[0x56822a]
[qi-Alienware-17-R4:08254] [20] /usr/bin/python3.8(PyEval_EvalCode+0x27)[0x68c1e7]
[qi-Alienware-17-R4:08254] [21] /usr/bin/python3.8[0x67d5a1]
[qi-Alienware-17-R4:08254] [22] /usr/bin/python3.8[0x67d61f]
[qi-Alienware-17-R4:08254] [23] /usr/bin/python3.8(PyRun_FileExFlags+0x9b)[0x67d6db]
[qi-Alienware-17-R4:08254] [24] /usr/bin/python3.8(PyRun_SimpleFileExFlags+0x17e)[0x67da6e]
[qi-Alienware-17-R4:08254] [25] /usr/bin/python3.8(Py_RunMain+0x212)[0x6b6132]
[qi-Alienware-17-R4:08254] [26] /usr/bin/python3.8(Py_BytesMain+0x2d)[0x6b64bd]
[qi-Alienware-17-R4:08254] [27] /lib/x86_64-linux-gnu/libc.so.6(__libc_start_main+0xf3)[0x7fbef6a940b3]
[qi-Alienware-17-R4:08254] [28] /usr/bin/python3.8(_start+0x2e)[0x5f927e]
[qi-Alienware-17-R4:08254] *** End of error message ***[/quote]

Is there something I should modify to enforce Dirichlet boundary conditions?
Best,
Qi

schruste · March 5, 2021, 3:17pm

Thanks for the report. This was a bug on my end. I only implemented the Restricted spaces for volume terms and access. The Set went into a trap then. This should be fixed now, see current master on github and have a try with that.

Best,
Christoph

Ryanleaf · March 6, 2021, 2:50am

Dear Christoph,

Thank you very much for fixing the bug. This update fixed the boundary issue.
I am sorry I think may meet another problem about multigrid.:pinch:
Same script:

[code]from netgen.geom2d import SplineGeometry
from netgen.meshing import MeshingParameters
from xfem import *
from netgen.csg import CSGeometry, OrthoBrick, Pnt
import netgen.gui
from netgen.geom2d import unit_square
from ngsolve import *
from ngsolve.la import EigenValues_Preconditioner
from xfem.lsetcurv import *

def eps(u):
return 0.5 * (Grad(u) + Grad(u).trans)
def Setup():
cube = CSGeometry()
cube.Add(OrthoBrick(Pnt(-1.5, -1.5, -1.5), Pnt(1.5, 1.5, 1.5)))
mesh = Mesh(cube.GenerateMesh(maxh=1/8, quad_dominated=False))
lsetmeshadap = LevelSetMeshAdaptation(mesh, order=3, threshold=1000, discontinuous_qn=True)
phi = Norm(CoefficientFunction((x, y, z))) - 1
deformation = lsetmeshadap.CalcDeformation(phi)
lsetp1 = lsetmeshadap.lset_p1
mesh.SetDeformation(deformation)
ci = CutInfo(mesh, lsetp1)
Vhbase = VectorH1(mesh, order=2, dirichlet=[1,2,3,4,5,6])
Vhneg = Restrict(Vhbase, ci.GetElementsOfType(HASNEG))
Vhpos = Restrict(Vhbase, ci.GetElementsOfType(HASPOS))
WhGV = FESpace([Vhneg, Vhpos], flags={“dgjumps”: True})
lset_neg = {“levelset”: lsetp1, “domain_type”: NEG}
lset_pos = {“levelset”: lsetp1, “domain_type”: POS}
lset_doms = [lset_neg, lset_pos]
elements_doms = [HASNEG, HASPOS]
u, v = WhGV.TnT()
a = BilinearForm(WhGV)
for i in [0, 1]: # loop over domains
a += SymbolicBFI(lset_doms[i], form=2* InnerProduct(eps(u[i]), eps(v[i])),
definedonelements=ci.GetElementsOfType(elements_doms[i]))
a += SymbolicBFI(lset_doms[i], form=u[i] * v[i], definedonelements=ci.GetElementsOfType(elements_doms[i]))
return mesh, WhGV, a
mesh, WhGV, a = Setup()
gfu = GridFunction(WhGV)
gfu.components[1].Set(CoefficientFunction((x,y,z)),BND)
mg = Preconditioner(a, “h1amg”) # register mg to a
a.Assemble()
print('NDOF = ', WhGV.ndof)
preI = Projector(mask=WhGV.FreeDofs(), range=True)
lams = EigenValues_Preconditioner(mat=a.mat, pre=mg.mat)
print(“condition number”, max(lams)/min(lams))[/code]
I noticed when I refine maxh to 1/8, I begin to get the following error:

[quote]Traceback (most recent call last):
File “/tmp/iterative_solver_krylov_hybrid.py/main.py”, line 42, in
a.Assemble()
netgen.libngpy._meshing.NgException: Inverse matrix: Matrix singularin Assemble BilinearForm ‘biform_from_py’
[/quote]
Am I having inappropriate settings in this script to causing this problem?
Would you please let me know where to find a description about which version of algebraic multigrid is currently implemented?
Again thank you very much for your kind help.

Best,
Qi

schruste · March 6, 2021, 9:02pm

Dear Qi,

I can not reproduce the singular matrix problem.

However, the performance of the preconditioner is indeed not satisfactory. I am not an expert on the implemented h1amg precond, but it is not specifically tailored for unfitted FEM situations. In your examples you have stabilization for bad cuts so that it seems not surprising that an AMG tailored for more standard H1-conforming fitted discretizations suffers from bad cut situations.

Best,
Christoph

Ryanleaf · March 8, 2021, 5:04am

Dear Christoph,

Thank you very much for the quick reply.
It’s a little bit strange. Previously I got this singular matrix problem on my personal PC.
Today I updated the ngsxfem&ngsolve in cluster and run the exact the same script for maxh=1/8 & maxh=1/12.

[code]from netgen.geom2d import SplineGeometry
from netgen.meshing import MeshingParameters
from xfem import *
from netgen.csg import CSGeometry, OrthoBrick, Pnt
import netgen.gui
from netgen.geom2d import unit_square
from ngsolve import *
from ngsolve.la import EigenValues_Preconditioner
from xfem.lsetcurv import *

def eps(u):
return 0.5 * (Grad(u) + Grad(u).trans)
def Setup():
cube = CSGeometry()
cube.Add(OrthoBrick(Pnt(-1.5, -1.5, -1.5), Pnt(1.5, 1.5, 1.5)))
mesh = Mesh(cube.GenerateMesh(maxh=1/12, quad_dominated=False))
lsetmeshadap = LevelSetMeshAdaptation(mesh, order=3, threshold=1000, discontinuous_qn=True)
phi = Norm(CoefficientFunction((x, y, z))) - 1
deformation = lsetmeshadap.CalcDeformation(phi)
lsetp1 = lsetmeshadap.lset_p1
mesh.SetDeformation(deformation)
ci = CutInfo(mesh, lsetp1)
Vhbase = VectorH1(mesh, order=2, dirichlet=[1,2,3,4,5,6])
Vhneg = Restrict(Vhbase, ci.GetElementsOfType(HASNEG))
Vhpos = Restrict(Vhbase, ci.GetElementsOfType(HASPOS))
WhGV = FESpace([Vhneg, Vhpos], flags={“dgjumps”: True})
lset_neg = {“levelset”: lsetp1, “domain_type”: NEG}
lset_pos = {“levelset”: lsetp1, “domain_type”: POS}
lset_doms = [lset_neg, lset_pos]
elements_doms = [HASNEG, HASPOS]
u, v = WhGV.TnT()
a = BilinearForm(WhGV)
for i in [0, 1]: # loop over domains
a += SymbolicBFI(lset_doms[i], form=2* InnerProduct(eps(u[i]), eps(v[i])),
definedonelements=ci.GetElementsOfType(elements_doms[i]))
a += SymbolicBFI(lset_doms[i], form=u[i] * v[i], definedonelements=ci.GetElementsOfType(elements_doms[i]))
return mesh, WhGV, a
mesh, WhGV, a = Setup()
gfu = GridFunction(WhGV)
gfu.components[1].Set(CoefficientFunction((x,y,z)),BND)
mg = Preconditioner(a, “h1amg”) # register mg to a
a.Assemble()
print('NDOF = ', WhGV.ndof)
preI = Projector(mask=WhGV.FreeDofs(), range=True)
lams = EigenValues_Preconditioner(mat=a.mat, pre=mg.mat)
print(“condition number”, max(lams)/min(lams))
[/code]
I gott this singular matrix problem in both experiments on cluster too.

[quote]
catch in AssembleBilinearform 2: Inverse matrix: Matrix singular
Traceback (most recent call last):
File “”, line 42, in netgen.libngpy._meshing.NgException: Inverse matrix: Matrix singularin Assemble BilinearForm ‘biform_from_py’[/quote]
Do you know any possible settings that may cause this problem?

Best,
Qi