Cannot interpolate a function on a unit interval with SEGV error

Discoverer · March 6, 2023, 7:21am

I want to test the SUPG method with a one dimensional problem. And when I interpolate the exact solution on a very fine grid with 500 cells(actually no matter how many cells, I still get the errors), as the codes show,

x=SpatialCoordinate(domain)
domain_e=mesh.create_interval(MPI.COMM_WORLD,500,[0,1])
V_e= fem.FunctionSpace(domain_e,("CG",1))
u_exact=fem.Function(V_e)
u_exact.interpolate(fem.Expression((exp(a[0]*x[0]/kappa)-1)/(exp(a[0]/kappa)-1),V_e.element.interpolation_points()))

I receive the error

[0]PETSC ERROR: ------------------------------------------------------------------------
[0]PETSC ERROR: Caught signal number 11 SEGV: Segmentation Violation, probably memory access out of range
[0]PETSC ERROR: Try option -start_in_debugger or -on_error_attach_debugger
[0]PETSC ERROR: or see https://petsc.org/release/faq/#valgrind
[0]PETSC ERROR: or try http://valgrind.org on GNU/linux and Apple Mac OS X to find memory corruption errors
[0]PETSC ERROR: configure using --with-debugging=yes, recompile, link, and run 
[0]PETSC ERROR: to get more information on the crash.
[0]PETSC ERROR: Run with -malloc_debug to check if memory corruption is causing the crash.
malloc(): invalid size (unsorted)
[LAPTOP-8EKSL8AG:05599] *** Process received signal ***
[LAPTOP-8EKSL8AG:05599] Signal: Aborted (6)
[LAPTOP-8EKSL8AG:05599] Signal code:  (-6)

So I dont know the reason. please help,thanks in advance!
The complete code is as follows:

from dolfinx import mesh,fem
from ufl import (inner,dot,div,grad,
                 TrialFunction,TestFunction,dx,
                 CellDiameter,Min,sqrt,lhs,rhs,SpatialCoordinate,exp
                    )
from mpi4py import MPI
from petsc4py import PETSc
import numpy as np

use_SUPG=True
N=16
k=1

domain=mesh.create_interval(MPI.COMM_WORLD,N,[0,1])

kappa=fem.Constant(domain,PETSc.ScalarType(5e-4))
a=fem.Constant(domain,np.array((1,),dtype=PETSc.ScalarType))

V=fem.FunctionSpace(domain,("CG",k))
u,v=TrialFunction(V),TestFunction(V)
f=fem.Constant(domain,PETSc.ScalarType(0))

res_Galerkin=(kappa*inner(grad(u),grad(v))+dot(a,grad(u))*v-f*v)*dx

h=CellDiameter(domain)
Cinv=fem.Constant(domain,PETSc.ScalarType(6*k*k))

tau=Min(h*h/(Cinv*kappa),h/(2*sqrt(dot(a,a))))

res_SUPG=tau*(-div(kappa*grad(u))+dot(a,grad(u))-f)*dot(a,grad(v))*dx

res= res_Galerkin
if use_SUPG:
    res +=res_SUPG

uh=fem.Function(V)

point0=mesh.locate_entities_boundary(domain,0,lambda x:np.isclose(x[0],0))
point1=mesh.locate_entities_boundary(domain,0,lambda x:np.isclose(x[0],1))
bc=fem.Constant(domain,PETSc.ScalarType(0))
bc0=fem.dirichletbc(bc,fem.locate_dofs_topological(V,0,point0),V)
bc=fem.Constant(domain,PETSc.ScalarType(1))
bc1=fem.dirichletbc(bc,fem.locate_dofs_topological(V,0,point1),V)

bd=[bc0,bc1]
A=lhs(res)
L=rhs(res)

problem=fem.petsc.LinearProblem(A,L,bd)
uh=problem.solve()

x=SpatialCoordinate(domain)

domain_e=mesh.create_interval(MPI.COMM_WORLD,500,[0,1])
V_e= fem.FunctionSpace(domain_e,("CG",1))
u_exact=fem.Function(V_e)


u_exact.interpolate(fem.Expression((exp(a[0]*x[0]/kappa)-1)/(exp(a[0]/kappa)-1),V_e.element.interpolation_points()))


e=uh-u_exact

print("H1 seminorm error = "
      +str(sqrt(fem.assemble_scalar(fem.form(dot(grad(e),grad(e))*dx)))))
print("L2 error = "
      +str(sqrt(fem.assemble_scalar(fem.form(e*e*dx)))))

x = V.tabulate_dof_coordinates()
x_order = np.argsort(x[:,0])

import matplotlib.pyplot as plt
plt.plot(x[x_order, 0], u_exact.x.array[x_order],label='exact solution')
plt.plot(x[x_order, 0],uh.x.array[x_order],label='numerical solution')

plt.legend()
plt.autoscale()
plt.show()

dokken · March 6, 2023, 8:42am

I am not sure what ufl-version you are using, as I cannot find Min as part of ufl.

Discoverer · March 6, 2023, 8:49am

Thank you for reply.
Maybe I mix up the old dolfin with dolfinx.The Min is in dolfin and I take it for granted that dolfinx also includes it.I try to use min_value instead but the error still occurs.

dokken · March 6, 2023, 8:55am

You need to step back a bit and go into detail on what you are trying to do.

It seems like you are trying to interpolate an expression defined on a coarse mesh onto a finer mesh.
This is not supported. You need to interpolate the function from the coarse to the fine mesh, and then in turn interpolate an expression containing this interpolated function.

Discoverer · March 6, 2023, 9:09am

Thank you very much.I know where I am wrong. I want to plot the exact solution with a higher resolution than the numerical resut.And I want to put them both in the same figure.It seems that I need to redefine the functionspace and parameter involved in exact solution and then plot them.

dokken · March 6, 2023, 10:15am

There are multiple things in your code that you need to be careful about.

Trying to assemble solutions from different grids in the same form:

Here you should rather use the exact solution on the coarse mesh.

Next up you have various overflow issues due to the fact that kappa is so small, so your exponential term goes to infinity.

Increasing kappa and implementing the various fixes I’ve suggested, you have a working code:

from IPython import embed
import matplotlib.pyplot as plt
from dolfinx import mesh, fem
from functools import partial
from ufl import (inner, dot, div, grad,
                 TrialFunction, TestFunction, dx,
                 CellDiameter, min_value, sqrt, lhs, rhs, SpatialCoordinate, exp
                 )
from mpi4py import MPI
from petsc4py import PETSc
import numpy as np

use_SUPG = True
N = 16
k = 1

domain = mesh.create_interval(MPI.COMM_WORLD, N, [0, 1])

kappa = fem.Constant(domain, PETSc.ScalarType(5e-1))
a = fem.Constant(domain, np.array((1,), dtype=PETSc.ScalarType))

V = fem.FunctionSpace(domain, ("CG", k))
u, v = TrialFunction(V), TestFunction(V)
f = fem.Constant(domain, PETSc.ScalarType(0))

res_Galerkin = (kappa*inner(grad(u), grad(v))+dot(a, grad(u))*v-f*v)*dx

h = CellDiameter(domain)
Cinv = fem.Constant(domain, PETSc.ScalarType(6*k*k))

tau = min_value(h*h/(Cinv*kappa), h/(2*sqrt(dot(a, a))))

res_SUPG = tau*(-div(kappa*grad(u))+dot(a, grad(u))-f)*dot(a, grad(v))*dx

res = res_Galerkin
if use_SUPG:
    res += res_SUPG

uh = fem.Function(V)

point0 = mesh.locate_entities_boundary(
    domain, 0, lambda x: np.isclose(x[0], 0))
point1 = mesh.locate_entities_boundary(
    domain, 0, lambda x: np.isclose(x[0], 1))
bc = fem.Constant(domain, PETSc.ScalarType(0))
bc0 = fem.dirichletbc(bc, fem.locate_dofs_topological(V, 0, point0), V)
bc = fem.Constant(domain, PETSc.ScalarType(1))
bc1 = fem.dirichletbc(bc, fem.locate_dofs_topological(V, 0, point1), V)

bd = [bc0, bc1]
A = lhs(res)
L = rhs(res)

problem = fem.petsc.LinearProblem(A, L, bd)
uh = problem.solve()

x = SpatialCoordinate(domain)

domain_e = mesh.create_interval(MPI.COMM_WORLD, 500, [0, 1])


def u_exact(x, a, kappa, exponential_func):
    return (exponential_func(a[0]*x[0]/kappa)-1)/(exponential_func(a[0]/float(kappa))-1)


V_e = fem.FunctionSpace(domain_e, ("Lagrange", 1))
u_ex = fem.Function(V_e)
u_ex.interpolate(partial(u_exact, a=a.value,
                 kappa=float(kappa), exponential_func=np.exp))
x_fine = V_e.tabulate_dof_coordinates()
x_e_order = np.argsort(x_fine[:, 0])

e = uh-u_exact(x, a, kappa, exp)
print("H1 seminorm error = "
      + str(sqrt(fem.assemble_scalar(fem.form(dot(grad(e), grad(e))*dx)))))
print("L2 error = "
      + str(sqrt(fem.assemble_scalar(fem.form(e*e*dx)))))

x = V.tabulate_dof_coordinates()
x_order = np.argsort(x[:, 0])

plt.plot(x_fine[x_e_order, 0], u_ex.x.array[x_e_order],
         "-ro", label='exact solution')
plt.plot(x[x_order, 0], uh.x.array[x_order], "-bs", label='numerical solution')

plt.legend()
plt.autoscale()
plt.savefig("test_fig.png")

For further questions, please carefully review the code I’ve posted and make sure you understand it.

Discoverer · March 6, 2023, 1:48pm

Thank you for your patience and really great advices!
So if I understand right, for this line,

u_ex.interpolate(partial(u_exact, a=a.value,
                 kappa=float(kappa), exponential_func=np.exp))

this time the x which the function interpolate supply to the function u_exact is equal to x=SpatialCoordinate(domain_e)?
Then I want to ask that when defining
kappa = fem.Constant(domain, PETSc.ScalarType(5e-1)) a = fem.Constant(domain, np.array((1,), dtype=PETSc.ScalarType))
the domain on which these parameters are defined has no effect on interpolate()?Since as in your codes, the a,kappa is on coarse grid domain when used in u_ex.interpolate(partial(u_exact, a=a.value, kappa=float(kappa), exponential_func=np.exp)) while the u_ex is on fine grid domain_e.

Further when I use your codes on my computer , when running this line

u_ex.interpolate(partial(u_exact, a=a.value,
                 kappa=float(kappa), exponential_func=np.exp))

it sends error about float(),something like this,

 f = float(self)
  File "/usr/lib/python3/dist-packages/ufl/core/expr.py", line 318, in __float__
    v = float(self._ufl_evaluate_scalar_())
  File "/usr/lib/python3/dist-packages/ufl/core/expr.py", line 313, in _ufl_evaluate_scalar_
    return self(())  # No known x

Do you know how to fix it?

dokken · March 6, 2023, 2:48pm

Discoverer:

it sends error about float(),something like this,

 f = float(self)
  File "/usr/lib/python3/dist-packages/ufl/core/expr.py", line 318, in __float__
    v = float(self._ufl_evaluate_scalar_())
  File "/usr/lib/python3/dist-packages/ufl/core/expr.py", line 313, in _ufl_evaluate_scalar_
    return self(())  # No known x

Do you know how to fix it?

What version of DOLFINx are you using?
This was fixed in: Fix infinite recursion for `dolfinx.fem.Constant.__float__()` calls by mikics · Pull Request #2482 · FEniCS/dolfinx · GitHub
as far as I can tell.

Discoverer:

So if I understand right, for this line,
u_ex.interpolate(partial(u_exact, a=a.value,
                 kappa=float(kappa), exponential_func=np.exp))
this time the x which the function interpolate supply to the function u_exact is equal to x=SpatialCoordinate(domain_e)?

Yes.

This is why I am calling a.value and float(kappa), as this does not include the fine mesh, just the values that are stored in the constant (either as a numpy array or as a single floating value.

Discoverer · March 6, 2023, 3:29pm

Thank you for reply. The version of dolfinx I use is 0.5.2, so maybe I should update it.

dokken · March 6, 2023, 3:43pm

Yes, you should update to 0.6.0-r1: Release v0.6.0 · FEniCS/dolfinx · GitHub

Discoverer · March 6, 2023, 3:59pm

Thats great.Sorry for another issue.Can I just usesudo apt update && apt install fenicsx to get the latest version or I need to download the source code on github?I am not very familiar with it.Thanks in advance for help.

dokken · March 6, 2023, 4:02pm

V0.6.0 is currently not available on apt-get. There has been some issues with keeping both DOLFIN and DOLFINx, see: ufl 2023.1 breaks old dolfin · Issue #151 · FEniCS/ufl · GitHub

I would suggest using conda or docker, see: GitHub - FEniCS/dolfinx: Next generation FEniCS problem solving environment

Discoverer · March 6, 2023, 4:05pm

Thank you very much.I am going to manage that.

Topic		Replies	Views
SEGV fault when interpolating function onto different mesh Errors mesh , dolfinx	5	173	February 12, 2024
Can't create mesh in parallel General	10	450	October 17, 2023
Nested Newton solver for mpc dolfinx	2	49	December 5, 2024
Segfault on when using MPIRUN with LU solver General	1	252	November 3, 2023
Interpolating a solution into a finer solution mesh	1	48	September 4, 2024

Cannot interpolate a function on a unit interval with SEGV error

Related topics