ngsPETSc Docs
Home Page: https://ngspetsc.readthedocs.io/en/latest/ngsPETSc.html
License: MIT License
Multi-purpose finite element library
Find the Open Source Community on https://ngsolve.org
Support & Services: https://cerbsim.com
When a curved firedrake-netgen mesh is created in serial, if you reload it in parallel via Checkpointing, then the boundary IDs are not labelled correctly. See the MFE below.
from firedrake import *
from firedrake.output import VTKFile
# Load the mesh and apply DirichletBCs (run this in parallel)
import os.path
if os.path.exists("out/mesh.h5"):
with CheckpointFile("out/mesh.h5", 'r') as file:
mesh = file.load_mesh("firedrake_default")
V = FunctionSpace(mesh, "CG", 1)
U = Function(V)
DirichletBC(V, -10.0, 1).apply(U)
DirichletBC(V, 10.0, 2).apply(U)
VTKFile("out/func.pvd").write(U)
# Build the mesh and save it (run this in serial)
else:
from netgen.occ import *
box = Box(Pnt(0,0,0), Pnt(1,1,1))
body = box.MakeFillet(box.edges, 0.1)
body.faces.Min(X).bc("left")
body.faces.Max(X).bc("right")
geo = OCCGeometry(body)
ngmesh = geo.GenerateMesh(maxh=0.15)
print(ngmesh.GetRegionNames(dim=2))
label_flags = {"labels" : {"left" : 1, "right" : 2, "default" : 3, "fillet" : 4}}
mesh = Mesh(Mesh(ngmesh, netgen_flags=label_flags).curve_field(2))
# Important: If one instead generates the mesh via the code commented below, no problems arise!
'''
mesh = UnitCubeMesh(10, 10, 10)
V = VectorFunctionSpace(mesh, "CG", 4)
newcoords = assemble(interpolate(mesh.coordinates + as_vector([0.25 * sin(np.pi*mesh.coordinates[1]), 0, 0]), V))
mesh = Mesh(newcoords)
'''
with CheckpointFile("out/mesh.h5", "w") as file:
file.save_mesh(mesh)
Running the above code twice (the first time in serial, and the second time in parallel), results in the boundary conditions being applied incorrectly, as shown below:
On the other hand, if you run it the second time in serial (i.e. you load the mesh in serial), or you don't curve the mesh
(i.e. you do mesh = Mesh(ngmesh, netgen_flags=label_flags) ) then the boundary conditions are correctly applied.
Furthermore, if you instead make a curved mesh in vanilla Firedrake (see the commented code above), then the boundary conditions are again correctly applied, even when loading the mesh in parallel.
Note also that, if you instead make a "curved" linear mesh (i.e. you do mesh = Mesh(Mesh(ngmesh, netgen_flags=label_flags).curve_field(1)) ) then you still get the problematic behaviour with boundary conditions being applied incorrectly.
Minimal working example, using the Firedrake-Netgen interface:
import math
from firedrake import *
import netgen.occ as ngocc
from firedrake.output import VTKFile
Cy1 = ngocc.Cylinder(ngocc.Pnt(-0.5, 0.5, 0.0), ngocc.X, r=0.25, h=0.5)
Cy2 = ngocc.Cylinder(ngocc.Pnt(-0.5, -0.5, 0.0), ngocc.X, r=0.25, h=0.5)
Cy3 = ngocc.Cylinder(ngocc.Pnt(0.0, 0.0, 0.0), ngocc.X, r=1.0, h=2.5)
Cn4 = ngocc.Cone(ngocc.gp_Ax2(ngocc.Pnt(2.5, 0.0, 0.0), ngocc.X), r1=1.0, r2=0.25, h=1.0, angle=2*math.pi)
Cy5 = ngocc.Cylinder(ngocc.Pnt(3.5, 0.0, 0.0), ngocc.X, r=0.25, h=1.0)
body = Cy1 + Cy2 + Cy3 + Cn4 + Cy5
body = body.MakeFillet((body.edges[-0.4 < ngocc.X])[ngocc.X < 4.4], 0.075)
geo = ngocc.OCCGeometry(body)
ngmsh = geo.GenerateMesh(maxh=0.5)
msh = Mesh(Mesh(ngmsh).curve_field(2))
VTKFile("mesh.pvd").write(msh)
The resulting fillet along the interface of the big cylinder and cone is not generated correctly. See the image below where I am visualizing the mesh in Paraview. The fillet has somehow been deformed.
Note that, if one does not curve the mesh, this issue does not seem to occur.
Upstream has changed the way MPI is wrapped
Running a simple notebook such as https://github.com/fem-on-colab/fem-on-colab/blob/main/ngsolve/test-ngspetsc.ipynb results in
RuntimeError Traceback (most recent call last)
<ipython-input-1-3a9ca44d9240> in <cell line: 1>()
----> 1 mesh = Mesh(unit_square.GenerateMesh(maxh=0.5).Distribute(COMM_WORLD))
RuntimeError: ng_openmpi.so: cannot open shared object file: No such file or directory
cc @UZerbinati: it's unclear to me whether the issue is in ngsPETSc, or in netgen itself.
I tried to run the testsuite with a complex PETSc installation, and got
FAILED test_eps.py::test_eps_ghep_eigfuncs - ZeroDivisionError: float division by zero
FAILED test_ksp.py::test_ksp_preonly_lu - assert 0.5333333333319396 < 0.0001
FAILED test_ksp.py::test_ksp_cg_gamg - assert 0.5333333333319396 < 0.0001
FAILED test_pc.py::test_pc_gamg - assert 0.5333333333319396 < 0.0001
FAILED test_pc.py::test_pc_hiptmaier_xu_sor - netgen.libngpy._meshing.NgException: BaseVector::operator/=: division by zero
FAILED test_pc.py::test_pc_hiptmaier_xu_bjacobi - netgen.libngpy._meshing.NgException: BaseVector::operator/=: division by zero
FAILED test_snes.py::test_snes_toy_newtonls - assert -6 in [4, 3, 2]
plus one of the SNES tests getting stuck at
7822 SNES Function norm 8.626434292221e-01
7823 SNES Function norm 8.626434292221e-01
7824 SNES Function norm 8.626434292221e-01
7825 SNES Function norm 8.626434292221e-01
7826 SNES Function norm 8.626434292221e-01
7827 SNES Function norm 8.626434292221e-01
7828 SNES Function norm 8.626434292221e-01
7829 SNES Function norm 8.626434292221e-01
7830 SNES Function norm 8.626434292221e-01
7831 SNES Function norm 8.626434292221e-01
7832 SNES Function norm 8.626434292221e-01
7833 SNES Function norm 8.626434292221e-01
7834 SNES Function norm 8.626434292221e-01
and forcing me to interrupt the test session.
I do realize that I may sound like a fool trying to use a complex PETSc since (to my limited understanding) ngsolve only uses real numbers. However, I am doing that because ngsPETSc is also employed as a dependency of other libraries (firedrake and FEniCSx), which may legitimately employ PETSc in complex mode.
Possible course of actions:
PETSc, orngsPETSc are supposed to work with complex PETSc (I guess, all parts interfacing firedrake and FEniCSx), and which parts are not supposed to work (I guess, those interfacing ngsolve).MFE below:
from firedrake import *
from firedrake.output import VTKFile
from netgen.occ import *
box = Box(Pnt(0, 0, 0), Pnt(1, 1, 1))
body = box.MakeFillet(box.edges, 0.1)
geo = OCCGeometry(body)
ngmesh = geo.GenerateMesh(maxh=0.1)
mesh = Mesh(Mesh(ngmesh).curve_field(2))
VTKFile("mesh.pvd").write(mesh)
Running this code in serial, the mesh is created with no issues.
But running it in parallel (e.g. mpiexec -n 2 python3 test.py) results in the following warning:
firedrake:WARNING Not able to curve element 4724, residual is: 0.18171646033208677
So it seems some of the elements are not being curved in this case.
React
A declarative, efficient, and flexible JavaScript library for building user interfaces.
🖖 Vue.js is a progressive, incrementally-adoptable JavaScript framework for building UI on the web.
Typescript
TypeScript is a superset of JavaScript that compiles to clean JavaScript output.
An Open Source Machine Learning Framework for Everyone
Django
The Web framework for perfectionists with deadlines.
Laravel
A PHP framework for web artisans
Bring data to life with SVG, Canvas and HTML. 📊📈🎉
JavaScript (JS) is a lightweight interpreted programming language with first-class functions.
Some thing interesting about web. New door for the world.
A server is a program made to process requests and deliver data to clients.
Machine learning is a way of modeling and interpreting data that allows a piece of software to respond intelligently.
Some thing interesting about visualization, use data art
Some thing interesting about game, make everyone happy.
Facebook
We are working to build community through open source technology. NB: members must have two-factor auth.
Microsoft
Open source projects and samples from Microsoft.
Google
Google ❤️ Open Source for everyone.
Alibaba
Alibaba Open Source for everyone
D3
Data-Driven Documents codes.
Tencent
China tencent open source team.