Problem projecting with scaled meshes

dolfinx
1

I am trying to calculate electric field and charge density from poisson equation using a projection function. The code gives NaN for meshes with scale =1E-6 (this is the actual scale of the problem), but works for scale =1. I am able to solve poission equation for both meshes, but projection does not seem to work.

from mpi4py import MPI
from dolfinx import mesh
import ufo
import numpy,gmsh
from dolfinx.io import gmshio, XDMFFile
from dolfinx import fem
from dolfinx.fem import functionspace
from dolfinx import default_scalar_type
from dolfinx.fem.petsc import LinearProblem
import pyvista


scale =1E-6
L = 600.0*scale #Length 
W = 200.0*scale #Width 
H=50*scale  # Thickness
model_rank = 0
mesh_comm = MPI.COMM_WORLD

gdim=3


def simpleproject(e,V,bc=[]):
    """Simple Projection function"""
    dx = ufl.dx(V.mesh)
    # Define variational problem for projection
    v = ufl.TrialFunction(V)
    w = ufl.TestFunction(V)
    sol=fem.Function(V)
    a = ufl.dot(v, w) * dx
    L = ufl.dot(ufl.grad(e), w) * dx
    problem = LinearProblem(a,L,bc,petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
    sol=problem.solve()
    return sol


gmsh.initialize()
gmsh.model.add("Structure")

box=gmsh.model.occ.addBox(0,0,0,L,W,H,1)
gmsh.model.occ.synchronize()
gmsh.model.addPhysicalGroup(gdim,[3],tag=1)
boundaries = gmsh.model.getBoundary([(gdim,3)], oriented=False)
boundary_ids = [b[1] for b in boundaries]
gmsh.model.addPhysicalGroup(gdim-1,boundary_ids,tag=20)
gmsh.model.mesh.generate(gdim)

domain,_ ,_  = gmshio.model_to_mesh(gmsh.model, mesh_comm, model_rank, gdim=gdim)
gmsh.finalize()

tdim=domain.topology.dim
fdim=tdim-1
p=1
V = functionspace(domain, ("Lagrange", p))
V2= functionspace(domain, ("Lagrange", p,(tdim,)))
uD = fem.Function(V)
uD.interpolate(lambda x: 1 + x[0]**2 + 2 * x[1]**2)
Efield=simpleproject(-uD,V2)
Ez=Efield.sub(gdim-1).collapse()

from dolfinx import plot
u_topology, u_cell_types, u_geometry = plot.vtk_mesh(V)
u_grid = pyvista.UnstructuredGrid(u_topology, u_cell_types, u_geometry)
u_grid.point_data["u"] = Ez.x.array.realu_grid.set_active_scalars("u") 
u_plotter = pyvista.Plotter() 
u_plotter.add_mesh(u_grid, show_edges=True) 
u_plotter.view_xy() 
if not pyvista.OFF_SCREEN: 
    u_plotter.show()

Screenshot 2024-09-06 at 4.19.08 PM
Screenshot 2024-09-06 at 4.19.08 PM962×722 45.2 KB

Screenshot 2024-09-06 at 4.19.47 PM
Screenshot 2024-09-06 at 4.19.47 PM956×680 72.5 KB

2

I would not advice you to scale the mesh to such a small order of magnitude as you introduce floating point error in your code.
For instance, for any integral, you need to compute the determinant of the Jacobian of a single cell.
This is equivalent to the volume of a cell. If each side of your Domain is 1e-6, and say you have 1 element in each direction, then the volume of a hexahedron is 1e-18, which is below double precision. The volume of a tetrahedron is of course even smaller.

Please rather rescale your equation, if you aim to work in micrometers.

1 Like
3

For further questions, please also put some effort into your minimal reproducible example. Currently, you have not given a code that is runnable and it has tons of unused imports.
I cannot even reproduce your nan example after trying to fix your code:

from mpi4py import MPI
import gmsh
from dolfinx.io import gmshio
from dolfinx import fem
from dolfinx.fem import functionspace
from dolfinx.fem.petsc import LinearProblem
import ufl

scale = 1e-6
L = 600.0 * scale  # Length
W = 200.0 * scale  # Width
H = 50 * scale  # Thickness
model_rank = 0
mesh_comm = MPI.COMM_WORLD

gdim = 3


def simpleproject(e, V, bc=[]):
    """Simple Projection function"""
    dx = ufl.dx(V.mesh)
    # Define variational problem for projection
    v = ufl.TrialFunction(V)
    w = ufl.TestFunction(V)
    sol = fem.Function(V)
    a = ufl.dot(v, w) * dx
    L = ufl.dot(ufl.grad(e), w) * dx
    problem = LinearProblem(
        a, L, bc, petsc_options={"ksp_type": "preonly", "pc_type": "lu"}
    )
    sol = problem.solve()
    return sol


gmsh.initialize()
gmsh.model.add("Structure")

box = gmsh.model.occ.addBox(0, 0, 0, L, W, H, 1)
gmsh.model.occ.synchronize()
gmsh.model.addPhysicalGroup(gdim, [3], tag=1)
boundaries = gmsh.model.getBoundary([(gdim, box)], oriented=False)
boundary_ids = [b[1] for b in boundaries]
gmsh.model.addPhysicalGroup(gdim - 1, boundary_ids, tag=20)
gmsh.model.mesh.generate(gdim)

domain, _, _ = gmshio.model_to_mesh(gmsh.model, mesh_comm, model_rank, gdim=gdim)
gmsh.finalize()

tdim = domain.topology.dim
fdim = tdim - 1
p = 1
V = functionspace(domain, ("Lagrange", p))
V2 = functionspace(domain, ("Lagrange", p, (gdim,)))
uD = fem.Function(V)
uD.interpolate(lambda x: 1 + x[0] ** 2 + 2 * x[1] ** 2)
Efield = simpleproject(-uD, V2)
print(Efield.x.array)
4

Sorry for the typos in the code…

from mpi4py import MPI
from dolfinx import mesh
import ufl
import numpy,gmsh
from dolfinx.io import gmshio
from dolfinx import fem
from dolfinx.fem import functionspace
from dolfinx import default_scalar_type
from dolfinx.fem.petsc import LinearProblem
import pyvista


scale =1E-6
L = 600.0*scale #Length 
W = 200.0*scale #Width 
H=50*scale  # Thickness
model_rank = 0
mesh_comm = MPI.COMM_WORLD

gdim=3


def simpleproject(e,V,bc=[]):
    """Simple Projection function"""
    dx = ufl.dx(V.mesh)
    # Define variational problem for projection
    v = ufl.TrialFunction(V)
    w = ufl.TestFunction(V)
    sol=fem.Function(V)
    a = ufl.dot(v, w) * dx
    L = ufl.dot(ufl.grad(e), w) * dx
    problem = LinearProblem(a,L,bc,petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
    sol=problem.solve()
    return sol


gmsh.initialize()
gmsh.model.add("Structure")

box=gmsh.model.occ.addBox(0,0,0,L,W,H,1)
gmsh.model.occ.synchronize()
gmsh.model.addPhysicalGroup(gdim,[1],tag=1)
boundaries = gmsh.model.getBoundary([(gdim,1)], oriented=False)
boundary_ids = [b[1] for b in boundaries]
gmsh.model.addPhysicalGroup(gdim-1,boundary_ids,tag=20)
gmsh.model.mesh.generate(gdim)

domain,_ ,_  = gmshio.model_to_mesh(gmsh.model, mesh_comm, model_rank, gdim=gdim)
gmsh.finalize()

tdim=domain.topology.dim
fdim=tdim-1
p=1
V = functionspace(domain, ("Lagrange", p))
V2= functionspace(domain, ("Lagrange", p,(tdim,)))
uD = fem.Function(V)
uD.interpolate(lambda x: 1 + x[0]**2 + 2 * x[1]**2)
Efield=simpleproject(-uD,V2)
Ez=Efield.sub(gdim-1).collapse()

from dolfinx import plot
u_topology, u_cell_types, u_geometry = plot.vtk_mesh(V)
u_grid = pyvista.UnstructuredGrid(u_topology, u_cell_types, u_geometry)
u_grid.point_data["u"] = Ez.x.array.real
u_grid.set_active_scalars("u") 
u_plotter = pyvista.Plotter() 
u_plotter.add_mesh(u_grid, show_edges=True) 
u_plotter.view_xy() 
if not pyvista.OFF_SCREEN: 
    u_plotter.show()
5

Scaling the equation did the trick. The problem was indeed the jacobian blowing up.

Thanks a lot!.