- Thank you received: 0
Algebraic Multigrid for Compressed VectorH1 space
3 years 9 months ago - 3 years 9 months ago #3538
by Ryanleaf
Algebraic Multigrid for Compressed VectorH1 space was created by Ryanleaf
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,.
When I run this part of code I immediately got the following segmentation fault,
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
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,.
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 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)
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
" 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
Last edit: 3 years 9 months ago by Ryanleaf.
3 years 9 months ago #3539
by schruste
Replied by schruste on topic Algebraic Multigrid for Compressed VectorH1 space
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
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
3 years 9 months ago - 3 years 9 months ago #3540
by Ryanleaf
Replied by Ryanleaf on topic Algebraic Multigrid for Compressed VectorH1 space
Thank you very much Christoph. This really solved my problem.
Sorry, One thing I just want to make sure.
Or I need to build a block preconditioner for the matrix?
Again, thank you very much.
Sorry, One thing I just want to make sure.
Suppose I have a two-phase problem, Can I use Algebraic Multigrid for the matrix from space WhGV?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})
Or I need to build a block preconditioner for the matrix?
Again, thank you very much.
Last edit: 3 years 9 months ago by Ryanleaf.
3 years 9 months ago #3541
by schruste
Replied by schruste on topic Algebraic Multigrid for Compressed VectorH1 space
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
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
3 years 9 months ago - 3 years 9 months ago #3544
by Ryanleaf
Replied by Ryanleaf on topic Algebraic Multigrid for Compressed VectorH1 space
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.
Best,
Qi
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))
Is this somehow the same problem with the previous one?H1AMG: level = 0, num_edges = 102939, nv = 1692
[Qis-MacBook-Pro:34374] *** Process received signal ***
[Qis-MacBook-Pro:34374] Signal: Segmentation fault: 11 (11)
[Qis-MacBook-Pro:34374] Signal code: (0)
[Qis-MacBook-Pro:34374] Failing at address: 0x0
[Qis-MacBook-Pro:34374] [ 0] 0 libsystem_platform.dylib 0x00007fff6ecf05fd _sigtramp + 29
[Qis-MacBook-Pro:34374] [ 1] 0 libsystem_c.dylib 0x00007fff95242d80 __c_locale + 0
[Qis-MacBook-Pro:34374] [ 2] 0 libngcomp.dylib 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
[Qis-MacBook-Pro:34374] [ 3] 0 libngcore.dylib 0x00000001020cd726 _ZN6ngcore11TaskManager9CreateJobERKNSt3__18functionIFvRNS_8TaskInfoEEEEi + 198
[Qis-MacBook-Pro:34374] [ 4] 0 libngcomp.dylib 0x000000010782ed14 _ZN6ngcomp12H1AMG_MatrixIdEC1ENSt3__110shared_ptrIN4ngla14SparseMatrixTMIdEEEENS3_IN6ngcore8BitArrayEEENS8_9FlatArrayINS8_3INTILi2EiEEmEENSB_IdmEESF_m + 4564
[Qis-MacBook-Pro:34374] [ 5] 0 libngcomp.dylib 0x0000000107839c8c _ZNSt3__110shared_ptrIN6ngcomp12H1AMG_MatrixIdEEE11make_sharedIJRNS0_IN4ngla14SparseMatrixTMIdEEEERNS0_IN6ngcore8BitArrayEEERNSB_5ArrayINSB_3INTILi2EiEEmEERNSF_IdmEESL_iEEES4_DpOT_ + 236
[Qis-MacBook-Pro:34374] [ 6] 0 libngcomp.dylib 0x0000000107838540 _ZN6ngcomp20H1AMG_PreconditionerIdE13FinalizeLevelEPKN4ngla10BaseMatrixE + 1120
[Qis-MacBook-Pro:34374] [ 7] 0 libngcomp.dylib 0x0000000107403de4 _ZN6ngcomp14S_BilinearFormIdE10DoAssembleERN6ngcore9LocalHeapE + 12180
[Qis-MacBook-Pro:34374] [ 8] 0 libngcomp.dylib 0x0000000107420044 _ZN6ngcomp12BilinearForm8AssembleERN6ngcore9LocalHeapE + 692
[Qis-MacBook-Pro:34374] [ 9] 0 libngcomp.dylib 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
[Qis-MacBook-Pro:34374] [10] 0 libngcomp.dylib 0x00000001077724b9 _ZN8pybind1112cpp_function10dispatcherEP7_objectS2_S2_ + 4041
[Qis-MacBook-Pro:34374] [11] 0 python 0x0000000100ea31b8 _PyMethodDef_RawFastCallKeywords + 392
[Qis-MacBook-Pro:34374] [12] 0 python 0x0000000100ea2c80 _PyObject_FastCallKeywords + 592
[Qis-MacBook-Pro:34374] [13] 0 python 0x0000000100fdfed5 call_function + 453
[Qis-MacBook-Pro:34374] [14] 0 python 0x0000000100fddaec _PyEval_EvalFrameDefault + 46092
[Qis-MacBook-Pro:34374] [15] 0 python 0x0000000100fd149e _PyEval_EvalCodeWithName + 414
[Qis-MacBook-Pro:34374] [16] 0 python 0x00000001010349a0 PyRun_FileExFlags + 256
[Qis-MacBook-Pro:34374] [17] 0 python 0x0000000101033e17 PyRun_SimpleFileExFlags + 391
[Qis-MacBook-Pro:34374] [18] 0 python 0x0000000101061d3f pymain_main + 9663
[Qis-MacBook-Pro:34374] [19] 0 python 0x0000000100e7566d main + 125
[Qis-MacBook-Pro:34374] [20] 0 libdyld.dylib 0x00007fff6eaf7cc9 start + 1
[Qis-MacBook-Pro:34374] *** End of error message ***
Best,
Qi
Last edit: 3 years 9 months ago by Ryanleaf.
3 years 9 months ago #3548
by schruste
Replied by schruste on topic Algebraic Multigrid for Compressed VectorH1 space
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:
1. switch to an XFEM formulation
2. add a very small regularization which is independent of the cut configuration, e.g.
3. Implement a version of the CompressedFESpace that also removes the link from elements to dofs where desired.
Best,
Christoph
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:
1. switch to an XFEM formulation
2. add a very small regularization which is independent of the cut configuration, e.g.
Code:
a += SymbolicBFI(form=1e-12 * u[i] * v[i])
Best,
Christoph
The following user(s) said Thank You: Ryanleaf
Time to create page: 0.113 seconds