Hello everyone,
I have two equations (for testing simple examples):
u1'=0,
u2'=0.
The boundary condition is
u1(0)=1, u2(0)=3
Therefore, the solution is “u1 = 1, u2 = 3”, however, we use the “NewtonSolver”, although the program runs, the results are wrong. (When only one equation or “u1(0) = u2(0)”, the result is correct). Could someone please let me know what I am doing wrong with my code. I think maybe there’s something wrong with the boundary conditions, but I don’t know how to change them. Here is the code
# bc
def left_boundary(x):
return np.isclose(x[0], 0.0)
def right_boundary(x):
return np.isclose(x[0], 1.0)
tdim = msh.topology.dim
msh.topology.create_connectivity(tdim-1, tdim)
left_facets = dolfinx.mesh.locate_entities_boundary(msh, tdim-1, left_boundary)
right_facets = dolfinx.mesh.locate_entities_boundary(msh, tdim-1, right_boundary)
ME0 = ME.sub(0)
ME1 = ME.sub(1)
V0, V0_to_ME0 = ME0.collapse()
V1, V1_to_ME1 = ME1.collapse()
u0 = dolfinx.fem.Function(V0)
u0.x.array[:] = 1
dofs0 = dolfinx.fem.locate_dofs_topological((ME0, V0), msh.topology.dim - 1, left_facets)
bc0 = dolfinx.fem.dirichletbc(u0, dofs0, ME0)
u1 = dolfinx.fem.Function(V0)
u1.x.array[:] = 3
dofs1 = dolfinx.fem.locate_dofs_topological((ME1, V1), msh.topology.dim - 1, left_facets)
bc1= dolfinx.fem.dirichletbc(u1, dofs1, ME1)
and the whole code is
import os
from PIL.ImageChops import constant
from fontTools.misc.transform import Scale
try:
from petsc4py import PETSc
import dolfinx
if not dolfinx.has_petsc:
print("This demo requires DOLFINx to be compiled with PETSc enabled.")
exit(0)
except ModuleNotFoundError:
print("This demo requires petsc4py.")
exit(0)
import ufl
from dolfinx.nls.petsc import NewtonSolver
from dolfinx import default_real_type
from dolfinx.io import XDMFFile
from dolfinx.fem import (functionspace,
dirichletbc,
Constant,
Function,
locate_dofs_topological,
locate_dofs_geometrical,
FunctionSpace, DirichletBC, )
from dolfinx.mesh import (create_interval,
locate_entities_boundary,)
from dolfinx.fem.petsc import NonlinearProblem, LinearProblem
from mpi4py import MPI
import numpy as np
from ufl import (SpatialCoordinate,
TrialFunctions,
TestFunctions,
inner,
grad,
dx,
sin)
from basix.ufl import element,mixed_element
import pandas as pd
import matplotlib.pyplot as plt
#Creat mesh
msh = create_interval(
MPI.COMM_WORLD,
nx=100,
points=[0, 1],
)
nodes = msh.geometry.x[:, 0]
P1 = element("Lagrange", msh.basix_cell(), 1, dtype=default_real_type)
P2 = element("Lagrange", msh.basix_cell(), 1, dtype=default_real_type)
ME = functionspace(msh, mixed_element([P1, P2]))
v1, v2 = ufl.TestFunctions(ME)
u = Function(ME) # current solution
u0 = Function(ME) # solution from previous converged step
# functions
u1, u2 = ufl.split(u)
u10, u20 = ufl.split(u0)
# Zero u
u.x.array[:] = 0.0
# Interpolate initial condition
u.sub(0).interpolate(lambda x: np.full_like(x[0], 1.0))
u.sub(1).interpolate(lambda x: np.full_like(x[0], 3.0))
u.x.scatter_forward()
x = SpatialCoordinate(msh)
#noslip = np.zeros(msh.geometry.dim, dtype=PETSc.ScalarType) # type: ignore
#facets = locate_entities_boundary(msh, msh.topology.dim - 1, noslip_boundary)
#bc0 = dirichletbc(noslip, locate_dofs_topological(V, msh.topology.dim - 1, facets), V)
# Weak statement of the equations
F1 = u1.dx(0) * v1 * dx
F2 = u2.dx(0) * v2 * dx
F = F1 + F2
# Collect Dirichlet boundary conditions
# bc
def left_boundary(x):
return np.isclose(x[0], 0.0)
def right_boundary(x):
return np.isclose(x[0], 1.0)
tdim = msh.topology.dim
msh.topology.create_connectivity(tdim-1, tdim)
left_facets = dolfinx.mesh.locate_entities_boundary(msh, tdim-1, left_boundary)
right_facets = dolfinx.mesh.locate_entities_boundary(msh, tdim-1, right_boundary)
ME0 = ME.sub(0)
ME1 = ME.sub(1)
V0, V0_to_ME0 = ME0.collapse()
V1, V1_to_ME1 = ME1.collapse()
u0 = dolfinx.fem.Function(V0)
u0.x.array[:] = 1
dofs0 = dolfinx.fem.locate_dofs_topological((ME0, V0), msh.topology.dim - 1, left_facets)
bc0 = dolfinx.fem.dirichletbc(u0, dofs0, ME0)
u1 = dolfinx.fem.Function(V0)
u1.x.array[:] = 3
dofs1 = dolfinx.fem.locate_dofs_topological((ME1, V1), msh.topology.dim - 1, left_facets)
bc1= dolfinx.fem.dirichletbc(u1, dofs1, ME1)
# Create nonlinear problem and Newton solver
problem = NonlinearProblem(F, u, bcs=[bc0, bc1])
#
solver = NewtonSolver(MPI.COMM_WORLD, problem)
solver.convergence_criterion = "incremental"
solver.rtol = 1e-8
solver.atol = 1e-10
solver.max_it = 25
# We can customize the linear solver used inside the NewtonSolver by
# modifying the PETSc options
ksp = solver.krylov_solver
opts = PETSc.Options() # type: ignore
option_prefix = ksp.getOptionsPrefix()
opts[f"{option_prefix}ksp_type"] = "preonly"
opts[f"{option_prefix}pc_type"] = "lu"
sys = PETSc.Sys() # type: ignore
# For factorisation prefer superlu_dist, then MUMPS, then default
if sys.hasExternalPackage("superlu_dist"):
opts[f"{option_prefix}pc_factor_mat_solver_type"] = "superlu_dist"
elif sys.hasExternalPackage("mumps"):
opts[f"{option_prefix}pc_factor_mat_solver_type"] = "mumps"
ksp.setFromOptions()
# solver
dolfinx.log.set_log_level(dolfinx.log.LogLevel.INFO)
num_iterations, converged = solver.solve(u)
assert converged, "Solver did not converge"
#
msh.topology.create_connectivity(1, 0)
nodes = msh.geometry.x[:, 0]
sorted_nodes = np.argsort(nodes)
#
u1_values = u.sub(0).x.array[:]
u2_values = u.sub(1).x.array[:]
#
rank = MPI.COMM_WORLD.rank
if rank == 0:
x_vals = nodes[sorted_nodes]
u1_sorted = u1_values[sorted_nodes]
u2_sorted = u2_values[sorted_nodes]
#
df = pd.DataFrame({'x': x_vals, 'u1': u1_sorted, 'u2': u2_sorted})
df.to_csv("solution.csv", index=False)
print("已保存为 solution.csv")
#
plt.figure(figsize=(8, 5))
plt.plot(x_vals, u1_sorted, label="u1", linewidth=2)
plt.plot(x_vals, u2_sorted, label="u2", linewidth=2)
plt.xlabel("x")
plt.ylabel("Solution")
plt.title("Solutions u1 and u2 vs x")
plt.legend()
plt.grid(True)
plt.tight_layout()
plt.savefig("solution_plot.png", dpi=300)
plt.show()
import sys
sys.exit(0)