search

ComputationalPhysiology / oasisx

Next generation Open Source Navier Stokes solver using FEniCSx
https://computationalphysiology.github.io/oasisx
MIT License
6 stars 1 forks source link

Add Pressure boundary condition #7

Closed jorgensd closed 1 year ago

jorgensd commented 1 year ago

Add consistent treatment of pressure from external forces on outlets.

jorgensd commented 1 year ago

Issue with tentative velocity term, i.e. u_ab is written in the wrong way. Reference code for simpler CN AB solver:


from IPython import embed
import ufl
import dolfinx
import numpy as np
from mpi4py import MPI
from petsc4py import PETSc

class NaiveScheme():

    def __init__(self, mesh, u_deg, p_deg, dt, nu):
        self._mesh = mesh
        self.V = dolfinx.fem.VectorFunctionSpace(mesh, ("Lagrange", u_deg))
        self.Q = dolfinx.fem.FunctionSpace(mesh, ("Lagrange", p_deg))
        self.u = dolfinx.fem.Function(self.V, name="u")
        self.ut = dolfinx.fem.Function(self.V, name="ut")
        self.u1 = dolfinx.fem.Function(self.V, name="u1")
        self.u2 = dolfinx.fem.Function(self.V, name="u2")
        self.p = dolfinx.fem.Function(self.Q, name="p")
        self.phi = dolfinx.fem.Function(self.Q, name="phi")
        self.create_forms(dt, nu)
        self.create_solvers()

    def create_forms(self, k, nu):

        u = ufl.TrialFunction(self.V)
        v = ufl.TestFunction(self.V)
        dt = dolfinx.fem.Constant(self._mesh, k)
        dx = ufl.Measure("dx", domain=self._mesh)
        u_ab = 1.5 * self.u1 - 0.5 * self.u2

        self.a1 = 1./dt * ufl.inner(u, v) * dx
        self.a1 += ufl.inner(0.5*ufl.grad(u)*u_ab, v)*dx
        self.a1 += nu/2*ufl.inner(ufl.grad(u), ufl.grad(v))*dx
        self.L1 = 1./dt * ufl.inner(self.u1, v)*ufl.dx
        self.L1 -= ufl.inner((0.5*ufl.grad(self.u1)*u_ab), v)*dx
        self.L1 -= nu/2*ufl.inner(ufl.grad(self.u1), ufl.grad(v))*dx
        self.L1 += self.p * ufl.div(v)*dx

        # u_avg = 0.5 * (u + self.u1)
        # F = 1./dt * ufl.inner((u - self.u1), v) * dx
        # F += ufl.inner(ufl.dot(u_ab, ufl.grad(u_avg)), v) * dx
        # F += nu * ufl.inner(ufl.grad(u_avg), ufl.grad(v)) * dx
        # F -= self.p*ufl.div(v)*dx
        # self.a1, self.L1 = ufl.system(F)

        # w_time = dolfinx.fem.Constant(self._mesh, PETSc.ScalarType(3. / (2. * dt)))
        # w_diffusion = dolfinx.fem.Constant(self._mesh, PETSc.ScalarType(nu))
        # self.a1 = (w_time * ufl.inner(u, v) + w_diffusion
        #            * ufl.inner(ufl.grad(u), ufl.grad(v))) * dx
        # self.L1 = ufl.inner(self.p, ufl.div(v)) * dx
        # self.L1 += dolfinx.fem.Constant(mesh, PETSc.ScalarType(1. / (2. * dt))) *\
        #     ufl.inner(4 * self.u1 - self.u2, v) * dx

        # BDF2 with implicit Adams-Bashforth
        # bs = 2 * self.u1-self.u2
        # self.a1 += ufl.inner(ufl.grad(u) * bs, v) * dx
        # self.a1 += 0.5 * ufl.div(bs) * ufl.inner(u, v) * dx

        p = ufl.TrialFunction(self.Q)
        q = ufl.TestFunction(self.Q)
        self.a2 = ufl.inner(ufl.grad(p), ufl.grad(q))*ufl.dx
        self.L2 = -1./dt*ufl.div(self.ut)*q*ufl.dx

        # self.a2 = ufl.inner(ufl.grad(p), ufl.grad(q)) * dx
        # self.L2 = - w_time * ufl.inner(ufl.div(self.ut), q) * dx
        nullspace = PETSc.NullSpace().create(constant=True)

        F3 = 1/dt*ufl.inner(u-self.ut, v)*dx + ufl.inner(ufl.grad(self.phi), v)*dx
        self.a3, self.L3 = ufl.system(F3)

    def create_solvers(self):
        boundary_facets = dolfinx.mesh.exterior_facet_indices(self._mesh.topology)
        dofs = dolfinx.fem.locate_dofs_topological(
            self.V, self._mesh.topology.dim-1, boundary_facets)
        self.u_bc = dolfinx.fem.Function(self.V)
        bc = dolfinx.fem.dirichletbc(self.u_bc, dofs)
        self.problem1 = dolfinx.fem.petsc.LinearProblem(self.a1, self.L1, [bc], self.ut, petsc_options={"ksp_type": "preonly",
                                                                                                        "pc_type": "lu"})

        self.A2 = dolfinx.fem.petsc.assemble_matrix(dolfinx.fem.form(self.a2))
        nullspace = PETSc.NullSpace().create(constant=True, comm=self._mesh.comm)
        self.A2.setNearNullSpace(nullspace)
        self.A2.setNullSpace(nullspace)
        self.A2.assemble()
        self.problem2 = PETSc.KSP().create(self._mesh.comm)
        self.problem2.setOperators(self.A2)
        # opts = PETSc.Options()
        # opts["ksp_type"] = "preonly"
        # opts["pc_type"] = "lu"
        # opts["pc_factor_mat_solver_type"] = "mumps"
        # opts["mat_mumps_icntl_24"] = 1
        # opts["mat_mumps_icntl_25"] = 0
        # opts["ksp_error_if_not_converged"] = True
        self.problem2.setFromOptions()

        self.problem3 = dolfinx.fem.petsc.LinearProblem(self.a3, self.L3, [], self.u,  petsc_options={"ksp_type": "preonly",
                                                                                                      "pc_type": "lu"})

    def update_bc(self, bc_expr):
        self.u_bc.interpolate(bc_expr)

    def solve(self):
        uh = self.problem1.solve()

        b2 = dolfinx.fem.petsc.assemble_vector(dolfinx.fem.form(self.L2))
        b2.ghostUpdate(addv=PETSc.InsertMode.ADD, mode=PETSc.ScatterMode.REVERSE)

        self.A2.getNullSpace().remove(b2)
        self.problem2.solve(b2, self.phi.vector)
        self.phi.x.scatter_forward()
        vol = mesh.comm.allreduce(
            dolfinx.fem.assemble_scalar(dolfinx.fem.form(1*ufl.dx(domain=self._mesh))), op=MPI.SUM)
        phi_avg = mesh.comm.allreduce(
            dolfinx.fem.assemble_scalar(dolfinx.fem.form(self.phi*ufl.dx)),
            op=MPI.SUM)/vol
        self.phi.x.array[:] -= phi_avg
        u3 = self.problem3.solve()
        return (uh, self.phi, u3)

class U():
    def __init__(self, t, nu):
        self.t = t
        self.nu = nu

    def eval_x(self, x):
        return - np.cos(np.pi * x[0]) * np.sin(np.pi * x[1]) * np.exp(-2.0 * self.nu * np.pi**2 * self.t)

    def eval_y(self, x):
        return np.cos(np.pi * x[1]) * np.sin(np.pi * x[0]) * np.exp(-2.0 * self.nu * np.pi**2 * self.t)

    def eval(self, x):

        return (self.eval_x(x), self.eval_y(x))

N = 25
mesh = dolfinx.mesh.create_rectangle(
    MPI.COMM_WORLD, [[-1, -1], [1, 1]], [N, N], cell_type=dolfinx.mesh.CellType.triangle)
dim = mesh.topology.dim - 1
mesh.topology.create_connectivity(dim, dim+1)
facets = dolfinx.mesh.exterior_facet_indices(mesh.topology)
value = 3
values = np.full_like(facets, value, dtype=np.int32)
sort = np.argsort(facets)
facet_tags = dolfinx.mesh.meshtags(mesh, dim, facets[sort], values[sort])

dt = 1e-2
nu = 0.01
T = 1
solver = NaiveScheme(mesh, 2, 1, dt, nu)

# Create expression for error calculation and initial conditions
x = ufl.SpatialCoordinate(mesh)
p_time = dolfinx.fem.Constant(mesh, -dt/2)
man_p = -0.25 * (ufl.cos(2*ufl.pi*x[0])+ufl.cos(2*ufl.pi*x[1]))*ufl.exp(-4*ufl.pi**2*nu*p_time)
p_expr = dolfinx.fem.Expression(man_p, solver.Q.element.interpolation_points())

u_time = dolfinx.fem.Constant(mesh, 0.)
man_u = ufl.as_vector((
    - ufl.sin(ufl.pi*x[1])*ufl.cos(ufl.pi*x[0]),
    ufl.sin(ufl.pi*x[0])*ufl.cos(ufl.pi*x[1])))*ufl.exp(-2*ufl.pi**2*nu*u_time)

u_expr = dolfinx.fem.Expression(man_u, solver.V.element.interpolation_points())

# Create initial conditions at time 0 and -dt
u_time.value = -dt
solver.u2.interpolate(u_expr)
u_time.value = 0
solver.u1.interpolate(u_expr)
p_time.value = -dt/2.
solver.p.interpolate(p_expr)
bc_expr = U(0., nu)

# Create error computation
diff_u = solver.u-man_u
L2_u = dolfinx.fem.form(ufl.inner(diff_u, diff_u) * ufl.dx)
diff_p = solver.p - man_p
L2_p = dolfinx.fem.form(ufl.inner(diff_p, diff_p) * ufl.dx)

u_vtx = dolfinx.io.VTXWriter(mesh.comm, "u.bp", [solver.u])
ut_vtx = dolfinx.io.VTXWriter(mesh.comm, "ut.bp", [solver.ut])
p_vtx = dolfinx.io.VTXWriter(mesh.comm, "p.bp", [solver.p])
phi_vtx = dolfinx.io.VTXWriter(mesh.comm, "phi.bp", [solver.phi])
p_ex = dolfinx.fem.Function(solver.Q)
pex_vtx = dolfinx.io.VTXWriter(mesh.comm, "p_ex.bp", [p_ex])

i = 0

# p_vtx.write(float(0))
p_ex.interpolate(p_expr)
# pex_vtx.write(float(0))
while bc_expr.t <= T+1e-14:
    bc_expr.t += dt
    u_time.value += dt
    p_time.value += dt
    solver.update_bc(bc_expr.eval)

    u_t, phi, uh = solver.solve()

    solver.p.x.array[:] += phi.x.array
    solver.u2.x.array[:] = solver.u1.x.array[:]
    solver.u1.x.array[:] = solver.u.x.array[:]

    u_vtx.write(float(u_time.value))
    p_vtx.write(float(p_time.value))
    ut_vtx.write(float(u_time.value))
    phi_vtx.write(float(p_time.value))
    p_ex.interpolate(p_expr)
    pex_vtx.write(float(p_time.value))
    print(f"u_err {float(u_time.value):.2e}, {np.sqrt(mesh.comm.allreduce(dolfinx.fem.assemble_scalar(L2_u), op=MPI.SUM))}")
    print(f"p_err {float(p_time.value):.2e}, {np.sqrt(mesh.comm.allreduce(dolfinx.fem.assemble_scalar(L2_p), op=MPI.SUM))}")
    print("*"*10)

pex_vtx.close()
u_vtx.close()
p_vtx.close()
ut_vtx.close()
jorgensd commented 1 year ago 
from IPython import embed
import ufl
import dolfinx
import numpy as np
from mpi4py import MPI
from petsc4py import PETSc
import argparse
import logging

desc = "Taylor-Green convergence demo"
parser = argparse.ArgumentParser(description=desc,
                                 formatter_class=argparse.ArgumentDefaultsHelpFormatter)
parser.add_argument("-N", "--refinement", type=int, dest="Ns", action="append",
                    help="The number of elements in x and y direction", required=True)
parser.add_argument("-T0", "--T-start", dest="T_start", type=float,
                    help="Start time of simulation", default=0)
parser.add_argument("-T1", "--T-end", dest="T_end", type=float,
                    help="End time of simulation", default=1)
parser.add_argument("-dt", dest="dt", type=float, help="Time step", default=0.1)
parser.add_argument("-nu", dest="nu", type=float, help="Kinematic viscosity", default=0.01)
parser.add_argument("-u", dest="u_deg", type=int, help="Degree of velocity space", default=2)
parser.add_argument("-p", dest="p_deg", type=int, help="Degree of pressure space", default=1)
parser.add_argument("-lm", "--low-memory", dest="lm", action="store_true",
                    default=False, help="Use low memory version of Oasisx")
inputs = parser.parse_args()

logger = logging.getLogger("Oasisx")

class NaiveScheme():

    def __init__(self, mesh, u_deg, p_deg, dt, nu):
        self._mesh = mesh
        self.V = dolfinx.fem.VectorFunctionSpace(mesh, ("Lagrange", u_deg))
        self.Q = dolfinx.fem.FunctionSpace(mesh, ("Lagrange", p_deg))
        self.u = dolfinx.fem.Function(self.V, name="u")
        self.ut = dolfinx.fem.Function(self.V, name="ut")
        self.u1 = dolfinx.fem.Function(self.V, name="u1")
        self.u2 = dolfinx.fem.Function(self.V, name="u2")
        self.p = dolfinx.fem.Function(self.Q, name="p")
        self.phi = dolfinx.fem.Function(self.Q, name="phi")
        self.create_forms(dt, nu)
        self.create_solvers()

    def create_forms(self, k, nu):

        u = ufl.TrialFunction(self.V)
        v = ufl.TestFunction(self.V)
        dt = dolfinx.fem.Constant(self._mesh, k)
        dx = ufl.Measure("dx", domain=self._mesh)
        u_ab = 1.5 * self.u1 - 0.5 * self.u2

        self.a1 = 1./dt * ufl.inner(u, v) * dx
        self.a1 += ufl.inner(0.5*ufl.grad(u)*u_ab, v)*dx
        self.a1 += nu/2*ufl.inner(ufl.grad(u), ufl.grad(v))*dx
        self.L1 = 1./dt * ufl.inner(self.u1, v)*ufl.dx
        self.L1 -= ufl.inner((0.5*ufl.grad(self.u1)*u_ab), v)*dx
        self.L1 -= nu/2*ufl.inner(ufl.grad(self.u1), ufl.grad(v))*dx
        self.L1 += self.p * ufl.div(v)*dx

        # u_avg = 0.5 * (u + self.u1)
        # F = 1./dt * ufl.inner((u - self.u1), v) * dx
        # F += ufl.inner(ufl.dot(u_ab, ufl.grad(u_avg)), v) * dx
        # F += nu * ufl.inner(ufl.grad(u_avg), ufl.grad(v)) * dx
        # F -= self.p*ufl.div(v)*dx
        # self.a1, self.L1 = ufl.system(F)

        # w_time = dolfinx.fem.Constant(self._mesh, PETSc.ScalarType(3. / (2. * dt)))
        # w_diffusion = dolfinx.fem.Constant(self._mesh, PETSc.ScalarType(nu))
        # self.a1 = (w_time * ufl.inner(u, v) + w_diffusion
        #            * ufl.inner(ufl.grad(u), ufl.grad(v))) * dx
        # self.L1 = ufl.inner(self.p, ufl.div(v)) * dx
        # self.L1 += dolfinx.fem.Constant(mesh, PETSc.ScalarType(1. / (2. * dt))) *\
        #     ufl.inner(4 * self.u1 - self.u2, v) * dx

        # BDF2 with implicit Adams-Bashforth
        # bs = 2 * self.u1-self.u2
        # self.a1 += ufl.inner(ufl.grad(u) * bs, v) * dx
        # self.a1 += 0.5 * ufl.div(bs) * ufl.inner(u, v) * dx

        p = ufl.TrialFunction(self.Q)
        q = ufl.TestFunction(self.Q)
        self.a2 = ufl.inner(ufl.grad(p), ufl.grad(q))*ufl.dx
        self.L2 = -1./dt*ufl.div(self.ut)*q*ufl.dx

        # self.a2 = ufl.inner(ufl.grad(p), ufl.grad(q)) * dx
        # self.L2 = - w_time * ufl.inner(ufl.div(self.ut), q) * dx
        nullspace = PETSc.NullSpace().create(constant=True)

        F3 = 1/dt*ufl.inner(u-self.ut, v)*dx + ufl.inner(ufl.grad(self.phi), v)*dx
        self.a3, self.L3 = ufl.system(F3)

    def create_solvers(self):
        boundary_facets = dolfinx.mesh.exterior_facet_indices(self._mesh.topology)
        dofs = dolfinx.fem.locate_dofs_topological(
            self.V, self._mesh.topology.dim-1, boundary_facets)
        self.u_bc = dolfinx.fem.Function(self.V)
        bc = dolfinx.fem.dirichletbc(self.u_bc, dofs)
        self.problem1 = dolfinx.fem.petsc.LinearProblem(self.a1, self.L1, [bc], self.ut, petsc_options={"ksp_type": "preonly",
                                                                                                        "pc_type": "lu"})

        self.A2 = dolfinx.fem.petsc.assemble_matrix(dolfinx.fem.form(self.a2))
        nullspace = PETSc.NullSpace().create(constant=True, comm=self._mesh.comm)
        self.A2.setNearNullSpace(nullspace)
        self.A2.setNullSpace(nullspace)
        self.A2.assemble()
        self.problem2 = PETSc.KSP().create(self._mesh.comm)
        self.problem2.setOperators(self.A2)
        opts = PETSc.Options()
        opts["ksp_type"] = "preonly"
        opts["pc_type"] = "lu"
        opts["pc_factor_mat_solver_type"] = "mumps"
        opts["mat_mumps_icntl_24"] = 1
        opts["mat_mumps_icntl_25"] = 0
        opts["ksp_error_if_not_converged"] = True
        self.problem2.setFromOptions()

        self.problem3 = dolfinx.fem.petsc.LinearProblem(self.a3, self.L3, [], self.u,  petsc_options={"ksp_type": "preonly",
                                                                                                      "pc_type": "lu"})

    def update_bc(self, bc_expr):
        self.u_bc.interpolate(bc_expr)

    def solve(self):
        uh = self.problem1.solve()

        b2 = dolfinx.fem.petsc.assemble_vector(dolfinx.fem.form(self.L2))
        b2.ghostUpdate(addv=PETSc.InsertMode.ADD, mode=PETSc.ScatterMode.REVERSE)

        self.A2.getNullSpace().remove(b2)
        self.problem2.solve(b2, self.phi.vector)
        self.phi.x.scatter_forward()
        vol = mesh.comm.allreduce(
            dolfinx.fem.assemble_scalar(dolfinx.fem.form(1*ufl.dx(domain=self._mesh))), op=MPI.SUM)
        phi_avg = mesh.comm.allreduce(
            dolfinx.fem.assemble_scalar(dolfinx.fem.form(self.phi*ufl.dx)),
            op=MPI.SUM)/vol
        self.phi.x.array[:] -= phi_avg
        u3 = self.problem3.solve()
        return (uh, self.phi, u3)

class U():
    def __init__(self, t, nu):
        self.t = t
        self.nu = nu

    def eval_x(self, x):
        return - np.cos(np.pi * x[0]) * np.sin(np.pi * x[1]) * np.exp(-2.0 * self.nu * np.pi**2 * float(self.t))

    def eval_y(self, x):
        return np.cos(np.pi * x[1]) * np.sin(np.pi * x[0]) * np.exp(-2.0 * self.nu * np.pi**2 * float(self.t))

    def eval(self, x):

        return (self.eval_x(x), self.eval_y(x))

dt = inputs.dt
nu = inputs.nu
assert (inputs.T_start < inputs.T_end)
T_end = inputs.T_end
T_start = inputs.T_start
num_steps = int((T_end-T_start)//dt)

assert inputs.u_deg > inputs.p_deg

solver_options = {"tentative": {"ksp_type": "preonly", "pc_type": "lu"},
                  "pressure": {"ksp_type": "preonly", "pc_type": "lu"},
                  "scalar": {"ksp_type": "preonly", "pc_type": "lu"}}

space_errors = np.empty((2, len(inputs.Ns)), dtype=np.float64)
hs = np.empty(len(inputs.Ns), dtype=np.float64)
for n, N in enumerate(inputs.Ns):

    mesh = dolfinx.mesh.create_rectangle(
        MPI.COMM_WORLD, [[-1, -1], [1, 1]], [N, N], cell_type=dolfinx.mesh.CellType.triangle)
    dim = mesh.topology.dim - 1
    mesh.topology.create_connectivity(dim, dim+1)
    facets = dolfinx.mesh.exterior_facet_indices(mesh.topology)
    value = 3
    values = np.full_like(facets, value, dtype=np.int32)
    sort = np.argsort(facets)
    facet_tags = dolfinx.mesh.meshtags(mesh, dim, facets[sort], values[sort])

    solver = NaiveScheme(mesh, inputs.u_deg, inputs.p_deg, dt, nu)

    # Create expression for error calculation and initial conditions
    x = ufl.SpatialCoordinate(mesh)
    p_time = dolfinx.fem.Constant(mesh, -dt/2)
    man_p = -0.25 * (ufl.cos(2*ufl.pi*x[0])+ufl.cos(2*ufl.pi*x[1]))*ufl.exp(-4*ufl.pi**2*nu*p_time)
    p_expr = dolfinx.fem.Expression(man_p, solver.Q.element.interpolation_points())

    u_time = dolfinx.fem.Constant(mesh, 0.)
    man_u = ufl.as_vector((
        - ufl.sin(ufl.pi*x[1])*ufl.cos(ufl.pi*x[0]),
        ufl.sin(ufl.pi*x[0])*ufl.cos(ufl.pi*x[1])))*ufl.exp(-2*ufl.pi**2*nu*u_time)

    u_expr = dolfinx.fem.Expression(man_u, solver.V.element.interpolation_points())

    # Create initial conditions at time 0 and -dt
    u_time.value = -dt
    solver.u2.interpolate(u_expr)
    u_time.value = 0
    solver.u1.interpolate(u_expr)
    p_time.value = -dt/2.
    solver.p.interpolate(p_expr)
    bc_expr = U(u_time, nu)

    # Create error computation
    diff_u = solver.u-man_u
    L2_u = dolfinx.fem.form(ufl.inner(diff_u, diff_u) * ufl.dx)
    diff_p = solver.p - man_p
    L2_p = dolfinx.fem.form(ufl.inner(diff_p, diff_p) * ufl.dx)

    u_vtx = dolfinx.io.VTXWriter(mesh.comm, "u.bp", [solver.u])
    ut_vtx = dolfinx.io.VTXWriter(mesh.comm, "ut.bp", [solver.ut])
    p_vtx = dolfinx.io.VTXWriter(mesh.comm, "p.bp", [solver.p])
    phi_vtx = dolfinx.io.VTXWriter(mesh.comm, "phi.bp", [solver.phi])
    p_ex = dolfinx.fem.Function(solver.Q)
    pex_vtx = dolfinx.io.VTXWriter(mesh.comm, "p_ex.bp", [p_ex])

    i = 0

    error_space_time = np.empty((2, num_steps), dtype=np.float64)
    for i in range(num_steps):
        u_time.value += dt
        p_time.value += dt
        solver.update_bc(bc_expr.eval)

        u_t, phi, uh = solver.solve()

        solver.p.x.array[:] += phi.x.array
        solver.u2.x.array[:] = solver.u1.x.array[:]
        solver.u1.x.array[:] = solver.u.x.array[:]

        L2_u_loc = dolfinx.fem.assemble_scalar(L2_u)
        error_u = mesh.comm.allreduce(L2_u_loc, op=MPI.SUM)
        L2_p_loc = dolfinx.fem.assemble_scalar(L2_p)
        error_p = mesh.comm.allreduce(L2_p_loc, op=MPI.SUM)
        logger.debug(f"{i=} {float(u_time.value)}, {error_u=}")
        logger.debug(f"{i=} {float(p_time.value)}, {error_p=}")

        error_space_time[:, i] = [error_u, error_p]
        logger.debug("*"*10)

    hmax_loc = np.max(mesh.h(mesh.topology.dim, np.arange(
        mesh.topology.index_map(mesh.topology.dim).size_local, dtype=np.int32)))
    hmax = mesh.comm.allreduce(hmax_loc, op=MPI.MAX)
    space_time_u_L2 = np.sqrt(dt*np.sum(error_space_time[0, :]))
    space_time_p_L2 = np.sqrt(dt*np.sum(error_space_time[1, :]))

    logger.setLevel(logging.INFO)
    logger.info(f"{hmax=} {space_time_u_L2=} {space_time_p_L2=}")
    hs[n] = hmax
    space_errors[:, n] = [space_time_u_L2, space_time_p_L2]

order = np.argsort(hs)[::-1]

logger.setLevel(logging.INFO)
logger.info(f"{hmax=} {space_time_u_L2=} {space_time_p_L2=}")
hs[n] = hmax
space_errors[:, n] = [space_time_u_L2, space_time_p_L2]

order = np.argsort(hs)[::-1]
hs = hs[order]

space_errors[0, :] = space_errors[0, order]
space_errors[1, :] = space_errors[1, order]
logger.info(
    f"Convergence rates u: {np.log(space_errors[0, 1:] / space_errors[0, :-1]) / np.log(hs[1:]/hs[:-1])}")
logger.info(
    f"Convergence rates p: {np.log(space_errors[1, 1:] / space_errors[1, :-1]) / np.log(hs[1:]/hs[:-1])}")

gives

root@dokken-XPS-9320:~/shared/oasisx/demo# python3 mwe2.py -N 15 -N 30 -N 60 -dt 0.01
INFO:Oasisx:hmax=0.18856180831641278 space_time_u_L2=0.009238574053141763 space_time_p_L2=0.011946135014718771
INFO:Oasisx:hmax=0.09428090415820647 space_time_u_L2=0.0005878265124763631 space_time_p_L2=0.002767463093457577
INFO:Oasisx:hmax=0.047140452079103314 space_time_u_L2=6.124957864698556e-05 space_time_p_L2=0.000681973256128087
INFO:Oasisx:hmax=0.047140452079103314 space_time_u_L2=6.124957864698556e-05 space_time_p_L2=0.000681973256128087
INFO:Oasisx:Convergence rates u: [3.97420786 3.26261861]
INFO:Oasisx:Convergence rates p: [2.10990795 2.02077701]