Fenicsx is excellent, but I have been tormented by the boundary conditions seting.
I have two equations
- \phi '' -4 A' \phi ' -8 \phi^3 +6 \phi^5 =0
- A '' - 2/3 \phi'^2 =0
the coordinate x \in [0,10] , with the boundary condations
A(0) =0, A'(0) = 0, \phi(0) =0, \phi'(10) =0
The bcs had been set by the following codes but I’m not sure if they are right.
Where \phi \to u1 , \ A \to u2 .
r_0 = 0.0
r_b = 10.0
# the coordinate of right boundary
msh = mesh.create_interval(comm=MPI.COMM_WORLD, points=(r_0, r_b), nx=100)
fdim = msh.topology.dim - 1
P1 = ufl.FiniteElement("Lagrange", msh.ufl_cell(), 1)
V = FunctionSpace(msh, ufl.MixedElement(P1, P1))
# the left boundary
facets_l = mesh.locate_entities_boundary(msh, dim=fdim, marker=lambda x: np.isclose(x[0], r_0))
dofs_u1 = fem.locate_dofs_topological(V=V.sub(0), entity_dim=fdim, entities=facets_l)
dofs_u2 = fem.locate_dofs_topological(V=V.sub(1), entity_dim=fdim, entities=facets_l)
# the Dirichleb boundary conditions
bc_u1 = fem.dirichletbc(value=ScalarType(0), dofs=dofs_u1, V=V.sub(0)) # u1 = phi
bc_u2 = fem.dirichletbc(value=ScalarType(0), dofs=dofs_u2, V=V.sub(1)) # u2 = A
bcs = [bc_u1, bc_u2]
The values of are dofs_u1 = [0] , dofs_u2 = [2] that confued me.
And the whole code are
import numpy as np
import ufl
from dolfinx import mesh, fem
from dolfinx.fem import Function, FunctionSpace
from dolfinx.fem.petsc import NonlinearProblem
from dolfinx.nls.petsc import NewtonSolver
from ufl import dx, split, TestFunctions
from mpi4py import MPI
from petsc4py import PETSc
from petsc4py.PETSc import ScalarType
# parmeters
r_0 = 0.0
r_b = 10.0 # the coordinate of right boundary
msh = mesh.create_interval(comm=MPI.COMM_WORLD, points=(r_0, r_b), nx=100)
fdim = msh.topology.dim - 1
P1 = ufl.FiniteElement("Lagrange", msh.ufl_cell(), 1)
V = FunctionSpace(msh, ufl.MixedElement(P1, P1))
# the left boundary
facets_l = mesh.locate_entities_boundary(
msh, dim=fdim, marker=lambda x: np.isclose(x[0], r_0))
dofs_u1 = fem.locate_dofs_topological(V=V.sub(0), entity_dim=fdim, entities=facets_l)
dofs_u2 = fem.locate_dofs_topological(V=V.sub(1), entity_dim=fdim, entities=facets_l)
# the Dirichleb boundary conditions
bc_u1 = fem.dirichletbc(value=ScalarType(0), dofs=dofs_u1, V=V.sub(0)) # u1 = phi
bc_u2 = fem.dirichletbc(value=ScalarType(0), dofs=dofs_u2, V=V.sub(1)) # u2 = A
bcs = [bc_u1, bc_u2]
# Define functions
u = Function(V)
# Split system functions to access components
#phi, m, sigma = u1, u2, u3
u1, u2 = split(u)
# Define test functions
v_1, v_2 = TestFunctions(V)
# Define variational problem
F1 = v_1.dx(0)*u1.dx(0) + (-4.0 * u2.dx(0)*u1.dx(0) - 8.0*u1**3+6.0*u1**5)*v_1
F2 = v_2.dx(0)*u2.dx(0) - 2.0/3.0*u1.dx(0)**2*v_2
F = (F1+F2)*dx
# Solve the system.
#
# t1 = data[:,1]
# t2 = data[:,2]
# t3 = data[:,3]
u.x.array[:] = np.random.random(202)
# Create nonlinear problem and Newton solver
problem = NonlinearProblem(F, u, bcs)
solver = NewtonSolver(MPI.COMM_WORLD, problem)
solver.convergence_criterion = "incremental"
solver.rtol = 1e-5
# We can customize the linear solver used inside the NewtonSolver by
# modifying the PETSc options
ksp = solver.krylov_solver
opts = PETSc.Options()
option_prefix = ksp.getOptionsPrefix()
opts[f"{option_prefix}ksp_type"] = "preonly"
opts[f"{option_prefix}pc_type"] = "lu"
opts[f"{option_prefix}pc_factor_mat_solver_type"] = "mumps"
ksp.setFromOptions()
# -
r = solver.solve(u)
The above codes can’t reach the a convegence result.
