Some questions about solving the nonlinear vector-valued problem

[Smith_Jack]

Hello.
I met a PETSc Segmentation fault error when I was trying to solve the following equations:

The backward Euler formula is used.

import ufl 
import dolfinx
from dolfinx.fem.petsc import NonlinearProblem
from dolfinx.nls.petsc import NewtonSolver
from petsc4py import PETSc
from mpi4py import MPI
import numpy as np 

t_ini = 0
t_fin = 2.0
n_step = 40
dt = (t_fin - t_ini) / n_step

nx = 32
L = 2.0
rank = MPI.COMM_WORLD.Get_rank()

domain = dolfinx.mesh.create_interval(MPI.COMM_WORLD, nx, points=(0.0, L))
P1 = ufl.FiniteElement("Lagrange", domain.ufl_cell(), 1)

#The VectorElement is used when the space is n-dim and the problem is also n-dim
#Here, the space is 1 dim, so we use MixedElement
V = dolfinx.fem.FunctionSpace(domain, ufl.MixedElement(P1, P1))

fdim = domain.topology.dim - 1

#right bounday condition is Dirchlet condition, so we extract the right boundary only.
#x[0] is the x-coordinate, x[1] is the y-coordinate (In this case, there is no y-coordinate)
boundary_facets = dolfinx.mesh.locate_entities_boundary(domain, fdim, lambda x: np.isclose(L, x[0]))

bc1 = dolfinx.fem.dirichletbc(PETSc.ScalarType(1.0), 
                              dolfinx.fem.locate_dofs_topological(V.sub(0), fdim, boundary_facets), 
                              V.sub(0))
bc2 = dolfinx.fem.dirichletbc(PETSc.ScalarType(0.0), 
                              dolfinx.fem.locate_dofs_topological(V.sub(1), fdim, boundary_facets), 
                              V.sub(1))
# u^{n}
u_n = dolfinx.fem.Function(V)
u_n.name = "u_n"
u_n1, u_n2 = u_n.split()
u_n1.interpolate(lambda x: np.ones(x[0].shape))
u_n2.interpolate(lambda x: np.zeros(x[0].shape))
# u^{n+1}
u = dolfinx.fem.Function(V)
u1, u2 = u.split()
v1, v2 = ufl.TestFunction(V)

G = -ufl.exp(u1-u2)+ufl.exp(-(u1-u2)) 
F1 = u1*v1*ufl.dx + dt*ufl.dot(ufl.grad(u1), ufl.grad(v1))*ufl.dx + dt*ufl.sin(u1)*v1*ufl.ds - dt*G*v1*ufl.dx - u_n1*v1*ufl.dx
F2 = u2*v2*ufl.dx + dt*ufl.dot(ufl.grad(u2), ufl.grad(v2))*ufl.dx + dt*(-0.1*u1*u2)*v2*ufl.ds + dt*G*v2*ufl.dx - u_n2*v2*ufl.dx
F = F1 + F2

problem = NonlinearProblem(F, u, bcs=[bc1, bc2])
solver = NewtonSolver(MPI.COMM_WORLD, problem)
solver.convergence_criterion = "incremental"
solver.rtol = 1e-6

ksp = solver.krylov_solver
opts = PETSc.Options()
option_prefix = ksp.getOptionsPrefix()
opts[f"{option_prefix}ksp_type"] = "cg"
opts[f"{option_prefix}pc_type"] = "gamg"
opts[f"{option_prefix}pc_factor_mat_solver_type"] = "mumps"
ksp.setFromOptions()

n, converged = solver.solve(u)

xdmf = dolfinx.io.XDMFFile(domain.comm, "diffusion.xdmf", "w")
xdmf.write_mesh(domain)

t = t_ini 
while t < t_fin:
    n, converged = solver.solve(u)
    if (rank == 0):
        print(f"Number of interations: {n:d}")
    u_n.x.array[:] = u.x.array
    t += dt
    xdmf.write_function(u, t)
xdmf.close()

Except this, I also have two other questions:
(1) What I really want to solve is a set of coupled diffusion-convection equations. The left-boundary condition are the fluxes of chemical species which are calculated by an ODE solver. How should I write the variational form when the boundary condition is numerical ？
(2) I tried a one component equation, the solution is correct but the XDMFFile writer can’t write the solution correctly when the program is runned in parallel. It can’t put the solution points in correct order. The following is my code, the solution in “diffusion.xdmf” file is correct when I run the script in serial, but it is in wrong order when I run the script in parallel.

import ufl 
import dolfinx
from dolfinx.fem.petsc import NonlinearProblem
from dolfinx.nls.petsc import NewtonSolver
from petsc4py import PETSc
from mpi4py import MPI
import numpy as np 

t_ini = 0
t_fin = 2.0
n_step = 20
dt = (t_fin - t_ini) / n_step
nx = 32
L = 2.0
rank = MPI.COMM_WORLD.Get_rank()

domain = dolfinx.mesh.create_interval(MPI.COMM_WORLD, nx, points=(0.0, L))
V = dolfinx.fem.FunctionSpace(domain, ("Lagrange", 1))
fdim = domain.topology.dim - 1
boundary_facets = dolfinx.mesh.locate_entities_boundary(domain, fdim, lambda x: np.isclose(L, x[0]))
bc = dolfinx.fem.dirichletbc(PETSc.ScalarType(1.0), dolfinx.fem.locate_dofs_topological(V, fdim, boundary_facets), V)

u_n = dolfinx.fem.Function(V)
u_n.name = "u_n"
u_n.interpolate(lambda x: np.ones(x[0].shape))

u = dolfinx.fem.Function(V)
v = ufl.TestFunction(V)
a = u*v*ufl.dx + dt*ufl.dot(ufl.grad(u), ufl.grad(v))*ufl.dx + dt*ufl.sin(u)*v*ufl.ds
L = u_n*v*ufl.dx
F = a - L 

problem = NonlinearProblem(F, u, bcs=[bc])
solver = NewtonSolver(MPI.COMM_WORLD, problem)
solver.convergence_criterion = "incremental"
solver.rtol = 1e-6

ksp = solver.krylov_solver
opts = PETSc.Options()
option_prefix = ksp.getOptionsPrefix()
opts[f"{option_prefix}ksp_type"] = "cg"
opts[f"{option_prefix}pc_type"] = "gamg"
opts[f"{option_prefix}pc_factor_mat_solver_type"] = "mumps"
ksp.setFromOptions()

xdmf = dolfinx.io.XDMFFile(domain.comm, "diffusion.xdmf", "w")
xdmf.write_mesh(domain)

t = t_ini 
while t < t_fin:
    n, converged = solver.solve(u)
    if (rank == 0):
        print(f"Number of interations: {n:d}")
    u_n.x.array[:] = u.x.array
    t += dt
    xdmf.write_function(u, t)
xdmf.close()

[dokken]

There are a few things wrong here:
1.

@Smith_Jack

u1, u2 = u.split()

should be u1, u2 = ufl.split(u).
Secondly, to use the XDMFFile, you should use a VectorElement. There is no real reason not to use a vector element, as you can supply the dimension to it:

@Smith_Jack

V = dolfinx.fem.FunctionSpace(domain, ufl.MixedElement(P1, P1))

should be
V = dolfinx.fem.FunctionSpace(domain, ufl.VectorElement("Lagrange", mesh.ufl_cell(), 1, dim=2)

[kumar]

Hello Mr.Dokken, I have a question regarding the above statement, is there any problem in using

P1 = ufl.FiniteElement("Lagrange", domain.ufl_cell(), 1)

Instead of the below for the scalar values considered to be dof,

ufl.VectorElement("Lagrange", mesh.ufl_cell(), 1, dim=2

Are they both the same for representing a scalar valued dof ? Or they make any differences?

Thanks for your continuous support.

[dokken]

The main difference is that if you create a mixed element of a set of similar elements, DOLFINx will not optimize the dof layout for them being identical in each component.
By using VectorElement, you are instructing DOLFINx that each component of the element uses the same basis function. This optimizes the dof-maps, index-maps etc over block size.

[kumar]

Thanks for your response. Sorry, but I have few questions regarding the above statements as it is new for me.

If that is the case, then in your provided code, dim=2 represents the two scalar values of the dof that we wanted to have. Also now here the vector element works as mixed space with the code shown below?

V = dolfinx.fem.FunctionSpace(domain, ufl.VectorElement("Lagrange", mesh.ufl_cell(), 1, dim=2)

Additionally Mr.Dokken, if someone have a dof of two scalar values and one vector values, lets say chemical potential(mu), concentration(c) were scalar and displacement(u) as vector. Then is it right to represent it like

V = ufl.VectorElement("Lagrange", mesh.ufl_cell(),1) #For deformation map
P1 = ufl.FiniteElement("Lagrange", mesh.ufl_cell(),1) 
ME = fem.FunctionSpace(domain, ufl.MixedElement([V, P1, P1]))

As I have no other idea than to use it like above code. Is it also needed to be corrected? Please help in the same regard.

Thanks for your support.

[dokken]

@kumar

If that is the case, then in your provided code, dim=2 represents the two scalar values of the dof that we wanted to have. Also now here the vector element works as mixed space with the code shown below?

Yes

@kumar

Additionally Mr.Dokken, if someone have a dof of two scalar values and one vector values, lets say chemical potential(mu), concentration(c) were scalar and displacement(u) as vector. Then is it right to represent it like

V = ufl.VectorElement("Lagrange", mesh.ufl_cell(),1) #For deformation map
P1 = ufl.FiniteElement("Lagrange", mesh.ufl_cell(),1) 
ME = fem.FunctionSpace(domain, ufl.MixedElement([V, P1, P1]))

As I have no other idea than to use it like above code. Is it also needed to be corrected? Please help in the same regard.

You can use the MixedElement as you suggest, it is not introducing any bugs, it just means that the layout will be different, compared to using ufl.VectorElement("Lagrange", mesh.ufl_cell(), 1, dim=4), which you could in turn split into your sub-vector and two scalar spaces.

[Smith_Jack]

Thanks for your detailed and kind reply.
For the XMDFFile writer, I said “it writes the solution in wrong order”, in fact it’s not the case. I just have to sort the mesh point when I read the .h5 file.

[Smith_Jack]

I paste my code here for someone who want a simple case of time-dependent problem.

# -*- coding:utf-8 -*-
import ufl 
import dolfinx
from dolfinx.fem.petsc import NonlinearProblem
from dolfinx.nls.petsc import NewtonSolver
from petsc4py import PETSc
from mpi4py import MPI
import numpy as np 

# Here, we solve a two-component vector-valued heat equation
#\begin{equation}
#    \begin{cases}
#    \frac{\partial u_1}{\partial t} &= \frac{\partial^2 u_1}{\partial x^2} - \exp(u_1 - u_2) + \exp(-(u_1 - u_2)) \\
#    \frac{\partial u_2}{\partial t} &= \frac{\partial^2 u_2}{\partial x^2} + \exp(u_1 - u_2) - \exp(-(u_1 - u_2)) \\
#    \end{cases}
#\end{equation}
#With the boundary condation:
#\begin{equation}
#    \begin{cases}
#        \frac{\partial u_1}{\partial x}(0, t) &= \sin(u_1(x = 0)) \\
#        \frac{\partial u_2}{\partial x}(0, t) &= -0.1u_1(x=0)u_2(x=0)
#    \end{cases}
#    \quad
#    \begin{cases}
#        u_1(1, t) = 1 \\
#        u_2(1, t) = 0
#    \end{cases}
#\end{equation}
#We use backward Euler formula

t_ini = 0
t_fin = 1.1
n_step = 11
dt = (t_fin - t_ini) / n_step

nx = 64
L = 1.0
rank = MPI.COMM_WORLD.Get_rank()

domain = dolfinx.mesh.create_interval(MPI.COMM_WORLD, nx, points=(0.0, L))
V = dolfinx.fem.FunctionSpace(domain, ufl.VectorElement("Lagrange", domain.ufl_cell(), 1, dim=2))

fdim = domain.topology.dim - 1

#right bounday condition is Dirchlet condition, so we extract the right boundary only.
#x[0] is the x-coordinate, x[1] is the y-coordinate (In this case, there is no y-coordinate)
boundary_facets = dolfinx.mesh.locate_entities_boundary(domain, fdim, lambda x: np.isclose(L, x[0]))

bc1 = dolfinx.fem.dirichletbc(PETSc.ScalarType(1.0), 
                              dolfinx.fem.locate_dofs_topological(V.sub(0), fdim, boundary_facets), 
                              V.sub(0))
bc2 = dolfinx.fem.dirichletbc(PETSc.ScalarType(0.0), 
                              dolfinx.fem.locate_dofs_topological(V.sub(1), fdim, boundary_facets), 
                              V.sub(1))


u_n = dolfinx.fem.Function(V)
u = dolfinx.fem.Function(V)
u_n1, u_n2 = ufl.split(u_n)
u1, u2 = ufl.split(u)

u_n.sub(0).interpolate(lambda x: np.ones(x[0].shape))
u_n.sub(1).interpolate(lambda x: np.zeros(x[0].shape))
u_n.x.scatter_forward()

v1, v2 = ufl.TestFunction(V)
G = -ufl.exp(u1-u2)+ufl.exp(-(u1-u2)) 
F1 = u1*v1*ufl.dx + dt*ufl.dot(ufl.grad(u1), ufl.grad(v1))*ufl.dx + dt*ufl.sin(u1)*v1*ufl.ds - dt*G*v1*ufl.dx - u_n1*v1*ufl.dx
F2 = u2*v2*ufl.dx + dt*ufl.dot(ufl.grad(u2), ufl.grad(v2))*ufl.dx + dt*(-0.1*u1*u2)*v2*ufl.ds + dt*G*v2*ufl.dx - u_n2*v2*ufl.dx
F = F1 + F2


problem = NonlinearProblem(F, u, bcs=[bc1, bc2])
solver = NewtonSolver(MPI.COMM_WORLD, problem)
solver.convergence_criterion = "residual"
solver.rtol = 1e-6

ksp = solver.krylov_solver
opts = PETSc.Options()
option_prefix = ksp.getOptionsPrefix()
opts[f"{option_prefix}ksp_type"] = "cg"
opts[f"{option_prefix}pc_type"] = "gamg"
opts[f"{option_prefix}pc_factor_mat_solver_type"] = "mumps"
ksp.setFromOptions()

file = dolfinx.io.XDMFFile(MPI.COMM_WORLD, "diffu_test3/u_tot.xdmf", "w")
file.write_mesh(domain)

t = t_ini 
while t < t_fin:
    n, converged = solver.solve(u)
    if (rank == 0):
        print(f"Number of interations: {n:d}")
    u_n.x.array[:] = u.x.array
    file.write_function(u, t)
    t += dt

file.close()