Solver returns inf on first step

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