Hello,
i advanced towards my goal, but i’m facing a difference in the solution found betwenn g(UFL) and g(Numpy).
I changed my FEniCSx version, I’m using FEniCSx 0.8.0,
My system is a square in 2d planar, where top border has uh = 1 (Dirichlet BC) and bottom border is Neumann BC with g = -uh². Using g in ufl format works fine, but i need to use a g from a numpy array.
I used the links provided by Dokken, thanks again they were usefull, to have a custom Newton solver and to use the external operator package.
With this 2 combined i was finally able to have convergence when i use g from numpy!
Unfortunately i have a non negligeable difference in the results after convergence, also the convergence rate is not the same see fig below :
i think there is something that i didn’t implement, but i really can’t see what! i would like advice on that, g (UFL) is working fine but how is it different from what i did with g (numpy)?
the MWE with g numpy and commmented g ufl:
Thanks for reading
"""
Created on Tue Jun 18 13:02:24 2024
In 2d planar square system, solve laplacian(uh) = 0 in Omega, uh = 1 in d_omega_D (top border), g = -uh² in d_omega_N (Neumann bottom border)
g use external operator package
"""
import dolfinx
from dolfinx import default_scalar_type, log, fem
from dolfinx.fem import (Constant, Function, functionspace,petsc,
assemble_scalar, dirichletbc, form, locate_dofs_geometrical,locate_dofs_topological)
from dolfinx.mesh import create_unit_square, locate_entities_boundary, meshtags, locate_entities
from dolfinx.plot import vtk_mesh
from mpi4py import MPI
import ufl
from ufl import * #SpatialCoordinate, TestFunction, TrialFunction, dot, ds, dx, grad, inner, Measure
from petsc4py.PETSc import ScalarType
from petsc4py import PETSc
import numpy as np
import pyvista
import matplotlib.pyplot as plt
from dolfinx_external_operator import (
FEMExternalOperator,
evaluate_external_operators,
evaluate_operands,
replace_external_operators)
import basix
def g_impl(uuh):
output = -1*pow(uuh,2) # NUMPY
return output.reshape(-1)# The output must be returned flattened to one dimension
def dgdu_impl(uuh):
aa = -2*uuh # NUMPY
return aa.reshape(-1) # The output must be returned flattened to one dimension
def g_external(derivatives):
if derivatives == (0,): # no derivation, the function itself
return g_impl
elif derivatives == (1,): # the derivative with respect to the operand `uh`
return dgdu_impl
else:
return NotImplementedError
##################### Code #######################################
mesh = create_unit_square(MPI.COMM_WORLD, 10, 10)
V = functionspace(mesh, ("CG", 1))
uh = Function(V)
v = TestFunction(V)
a = inner(grad(uh), grad(v)) * dx
x = SpatialCoordinate(mesh)
##################### ds fabrication for Neumann BC, bottom of square ####
boundaries = [(1, lambda x: x[1]<=0.001)]
facet_indices, facet_markers = [], []
fdim = mesh.topology.dim - 1
for (marker, locator) in boundaries:
facets = locate_entities(mesh, fdim, locator)
facet_indices.append(facets)
facet_markers.append(np.full_like(facets, marker))
facet_indices = np.hstack(facet_indices).astype(np.int32)
facet_markers = np.hstack(facet_markers).astype(np.int32)
sorted_facets = np.argsort(facet_indices)
facet_tag = meshtags(mesh, fdim, facet_indices[sorted_facets], facet_markers[sorted_facets])
ds = Measure("ds", domain=mesh, subdomain_data=facet_tag)
##################### Dirichlet BC on top +1V #####
dofs_up = locate_dofs_geometrical(V, lambda x: x[1]>=0.999)
u_up = Function(V)
u_up.interpolate(lambda x: x[0]*0+ 1)
bc_up = dirichletbc(u_up, dofs_up)
bcs = [bc_up]
##### dofs for NEUMANN BC
fdim = mesh.topology.dim - 1
facets = locate_entities_boundary(mesh, fdim, lambda m: m[1] <=0.001)
dofs = locate_dofs_topological(V, fdim, facets)
####### g(uh ufl) :
# g = -1*pow(uh,2) ### Full UFL formulation for comparaison purpose
##########################################################
####### g with external operator package :
quadrature_degree = 2
Qe = basix.ufl.quadrature_element(mesh.topology.cell_name(), degree=quadrature_degree, value_shape=())
Q = fem.functionspace(mesh, Qe)
dx = Measure("dx", metadata={"quadrature_scheme": "default", "quadrature_degree": quadrature_degree})
g = FEMExternalOperator(uh, function_space=Q) # g is now Symbolic not numpy involved
g.external_function = g_external
##########################################################
du = dolfinx.fem.Function(V)
maxiter = 10
i=0
correct_list=[]
L1 = inner(g,v) * ds # g symbolic used, this should be Neumann BC
F=a-L1
while i < maxiter:
######### Below work with g external operator
######### Assemble Jacobian and residual
F_replaced, F_external_operators = replace_external_operators(F)
evaluated_operands = evaluate_operands(F_external_operators)
_ = evaluate_external_operators(F_external_operators, evaluated_operands)
F_compiled = fem.form(F_replaced) # Residual
#################### Jacobian ###############################
J = derivative(F, uh)# define the derivative
J_expanded = ufl.algorithms.expand_derivatives(J)# actual value calculations of the derivative
J_replaced, J_external_operators = replace_external_operators(J_expanded)
_ = evaluate_external_operators(J_external_operators, evaluated_operands)
J_compiled = fem.form(J_replaced) # Jacobian var of the prgrm that works, jacobian = dolfinx.fem.form(J)
b_vector = fem.assemble_vector(F_compiled)
A_matrix = fem.assemble_matrix(J_compiled)
A = dolfinx.fem.petsc.create_matrix(J_compiled)
L = dolfinx.fem.petsc.create_vector(F_compiled)
solver = PETSc.KSP().create(mesh.comm)
solver.setOperators(A)
# # # dolfinx.fem.petsc.assemble_matrix(A, jacobian) ## replaced by :
dolfinx.fem.petsc.assemble_matrix(A, J_compiled, bcs=bcs)
A.assemble() ### important
# # # dolfinx.fem.petsc.assemble_vector(L, residual) ## replaced by :
dolfinx.fem.petsc.assemble_vector(L, F_compiled)
L.ghostUpdate(addv=PETSc.InsertMode.ADD, mode=PETSc.ScatterMode.REVERSE)
L.scale(-1)
# Compute b - J(u_D-u_(i-1))
dolfinx.fem.petsc.apply_lifting(L, [J_compiled], [bcs], x0=[uh.vector], scale=1)
# Set du|_bc = u_{i-1}-u_D
dolfinx.fem.petsc.set_bc(L, bcs, uh.vector, 1.0)
L.ghostUpdate(addv=PETSc.InsertMode.INSERT_VALUES, mode=PETSc.ScatterMode.FORWARD)
# Solve linear problem
solver.solve(L, du.vector) ## error if A = create matrix but not if A = assemble
du.x.scatter_forward()
# Update u_{i+1} = u_i + delta u_i
uh.x.array[:] += du.x.array
# Compute norm of update
correction_norm = du.vector.norm(0)
correct_list.append(correction_norm)
######## Below work with g UFL , to use comment above until while loop, uncomment g UFL and fully comment g external operator
# L1 = inner(g,v) * ds
# F=a-L1
# residual = dolfinx.fem.form(F)
# J = ufl.derivative(F, uh)
# jacobian = dolfinx.fem.form(J)
# A = dolfinx.fem.petsc.create_matrix(jacobian)
# L = dolfinx.fem.petsc.create_vector(residual)
# solver = PETSc.KSP().create(mesh.comm)
# solver.setOperators(A)
# dolfinx.fem.petsc.assemble_matrix(A, jacobian, bcs=bcs)
# A.assemble()
# dolfinx.fem.petsc.assemble_vector(L, residual)
# L.ghostUpdate(addv=PETSc.InsertMode.ADD, mode=PETSc.ScatterMode.REVERSE)
# L.scale(-1)
# # Compute b - J(u_D-u_(i-1))
# dolfinx.fem.petsc.apply_lifting(L, [jacobian], [bcs], x0=[uh.vector], scale=1)
# # Set du|_bc = u_{i-1}-u_D
# dolfinx.fem.petsc.set_bc(L, bcs, uh.vector, 1.0)
# L.ghostUpdate(addv=PETSc.InsertMode.INSERT_VALUES, mode=PETSc.ScatterMode.FORWARD)
# # Solve linear problem
# solver.solve(L, du.vector)
# du.x.scatter_forward()
# # Update u_{i+1} = u_i + delta u_i
# uh_old=uh.copy()
# uh.x.array[:] += du.x.array
# # Compute norm of update
# correction_norm = du.vector.norm(0)
# correct_list.append(correction_norm)
##########################################################################
print(f"Iteration {i}: Correction norm {correction_norm}")
i += 1
if correction_norm < 1e-8:
break
plt.figure(1)
plt.semilogy(range(i),correct_list,'-+')
plt.title('correction norm (iteration)')