Hello,
I am trying to solve a 2D system of equations with equal order approximation, starting with constant flow. First I tried with Dirichlet boundary conditions(https://fenicsproject.discourse.group/t/custom-newton-solver-how-to-implement-mixed-element-dirichlet-bcs/10332 ) and was unsucessful. I moved to all flux boundary conditions in an attempt to simplify the problem setup. However, the first newton step doesn’t make sense at all.
I am looking for guidance as to why the solver would return infinity on the first step. I simplified the code as much as possible:
import numpy as np
from dolfinx import mesh, fem, nls, plot, cpp,log
import ufl
from mpi4py import MPI
from petsc4py import PETSc
import pyvista
# problem parameters
rhoRef = 1
gamHt = 1.4
Pr = 0.72
Mach=1.1
R=288
cv = R/(gamHt - 1)
Pin = (rhoRef/gamHt)*(1/Mach**2)
dx = 0.2
dt = dx/2 # for CFL to hold
nx, ny = int(1/dx), int(1/dx)
#--------------------------------------------------------------------------------
# Domain
#--------------------------------------------------------------------------------
domain = mesh.create_unit_square(MPI.COMM_WORLD, nx, ny, mesh.CellType.triangle)
# Extract Boundaries
#--------------------------------------------------------------------------------
tdim = domain.topology.dim
fdim = tdim - 1
domain.topology.create_connectivity(fdim, tdim)
n = ufl.FacetNormal(domain)
boundaries = [(1, lambda x: np.isclose(x[0], 0)), # inlet
(2, lambda x: np.isclose(x[0], 1)), # outlet
(3, lambda x: np.isclose(x[1], 0)), # bottom
(4, lambda x: np.isclose(x[1], 1))] # top
facet_indices, facet_markers = [], []
for (marker, locator) in boundaries:
facets = mesh.locate_entities(domain, 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 = mesh.meshtags(domain, fdim, facet_indices[sorted_facets], facet_markers[sorted_facets])
# create boundary measures
ds = ufl.Measure("ds", domain=domain, subdomain_data=facet_tag)
#--------------------------------------------------------------------------------
# weak form setup
#--------------------------------------------------------------------------------
# Vector function space
W = fem.VectorFunctionSpace(domain, ("CG",1),dim=4)
# variables
U_n = fem.Function(W)
U_np1 = fem.Function(W)
w0, w1, w2, w3 = ufl.TestFunctions(W)
rho_n, ru1_n, ru2_n, E_n = ufl.split(U_n)
rho_np1, ru1_np1, ru2_np1 , E_np1 = ufl.split(U_np1)
# variables at midpoint
rho = (rho_n + rho_np1)/2
ru1 = (ru1_n + ru1_np1)/2
ru2 = (ru2_n + ru2_np1)/2
E = (E_n + E_np1)/2
# Pressure variable at midpoint
p = (gamHt-1)*(E - (rho/2)*( (ru1/rho)**2 + (ru2/rho)**2 ) )
# fluxes
F1 = ufl.as_vector([ru1,\
-((gamHt - 3)/2) * (ru1**2)/rho + (gamHt - 1)*(E - (ru2**2)/(2*rho)),\
ru1*ru2/rho,\
((gamHt - 1)/(2*rho**2)) * ru1*(ru1**2 + ru2**2) + gamHt*ru1*E/rho])
F2 = ufl.as_vector([ru2,\
ru1*ru2/rho,\
-((gamHt - 3)/2) * (ru2**2)/rho + (gamHt - 1)*(E - (ru1**2)/(2*rho)),\
((gamHt - 1)/(2*rho**2)) * ru2*(ru1**2 + ru2**2) + gamHt*ru2*E/rho])
# Residual
#--------------------------------------------------------------------------------
WtVec = ufl.as_vector([w0,w1,w2,w3])
# Galerkin
Res = ufl.inner(WtVec, (U_np1 - U_n)/dt)*ufl.dx
Res = -ufl.inner(WtVec.dx(0),F1)*ufl.dx - ufl.inner(WtVec.dx(1),F2)*ufl.dx
# boundaray terms
u1_in = fem.Constant(domain, PETSc.ScalarType(1))
u2_in = fem.Constant(domain, PETSc.ScalarType(0))
rho_in = fem.Constant(domain, PETSc.ScalarType(1))
p_in = fem.Constant(domain, PETSc.ScalarType(Pin))
E_in = fem.Constant(domain, PETSc.ScalarType(p_in/(gamHt - 1) + (rho_in/2)*(u1_in**2 + u2_in**2)))
n = ufl.FacetNormal(domain)
Fw2 = ufl.as_vector([0, 0, p, 0]) # p is defined in terms of solution vars above
Res += ufl.inner(WtVec, F1*n[0])*ds(2) # Outlet
Res += ufl.inner(WtVec,Fw2*n[1])*ds(3) # top
Res += ufl.inner(WtVec,Fw2*n[1])*ds(4) # bottom
F1_in = ufl.as_vector([rho_in*u1_in,\
-((gamHt - 3)/2) * ((rho_in*u1_in)**2)/rho_in + (gamHt - 1)*(E_in - ((rho_in*u2_in)**2)/(2*rho_in)),\
rho_in*u1_in*rho_in*u2_in/rho_in,\
((gamHt - 1)/(2*rho_in**2)) * rho_in*u1_in*((rho_in*u1_in)**2 + (rho_in*u2_in)**2) + gamHt*rho_in*E_in/rho_in])
Res += ufl.inner(WtVec,F1_in*n[0])*ds(1) # inlet
# initial conditions
U_n.sub(0).interpolate(fem.Expression(fem.Constant(domain, PETSc.ScalarType(1)), W.sub(0).element.interpolation_points()))
U_n.sub(1).interpolate(fem.Expression(fem.Constant(domain, PETSc.ScalarType(1)), W.sub(1).element.interpolation_points()))
U_n.sub(2).interpolate(fem.Expression(fem.Constant(domain, PETSc.ScalarType(0)), W.sub(2).element.interpolation_points()))
U_n.sub(3).interpolate(fem.Expression(E_in, W.sub(3).element.interpolation_points()))
#--------------------------------------------------------------------------------
# Time Loop
#--------------------------------------------------------------------------------
J = ufl.derivative(Res, U_np1)
residual = fem.form(Res)
jacobian = fem.form(J)
# Iteration dependent structures
A = fem.petsc.create_matrix(jacobian)
L = fem.petsc.create_vector(residual)
solver = PETSc.KSP().create(domain.comm)
solver.setOperators(A)
DeltaA = fem.Function(W)
#--------------------------------------------------------------------------------
# Time Loop
#--------------------------------------------------------------------------------
for niter in range(1):
# Solve NON-linear problem
#--------------------------------------------------------------------------------
A = fem.petsc.create_matrix(jacobian)
L = fem.petsc.create_vector(residual)
solver = PETSc.KSP().create(domain.comm)
solver.setOperators(A)
# begin newton loop
U_np1.x.array[:] = U_n.x.array
i = 0
while i < 10:
# Assemble Jacobian and residual
with L.localForm() as loc_L:
loc_L.set(0)
A.zeroEntries()
fem.petsc.assemble_matrix(A, jacobian)
A.assemble()
fem.petsc.assemble_vector(L, residual)
# Scale residual by -1
L.scale(-1)
# Solve
solver.solve(L, DeltaA.vector)
# Update
U_np1.x.array[:] += DeltaA.x.array
i += 1
# Compute norm of update
correction_norm = DeltaA.vector.norm(0)
print(f"Iteration {i}: Correction norm {correction_norm}")
if correction_norm < 1e-10:
break
# update for next time step
U_n.x.array[:] = U_np1.x.array
and the output of the correction_norm iterations is:
Iteration 1: Correction norm inf
Iteration 2: Correction norm inf
Iteration 3: Correction norm inf
Iteration 4: Correction norm inf
Iteration 5: Correction norm inf
Iteration 6: Correction norm inf
Iteration 7: Correction norm inf
Iteration 8: Correction norm inf
Iteration 9: Correction norm inf
Iteration 10: Correction norm inf