Add Pressure boundary condition #7

jorgensd · 2022-11-14T09:47:56Z

Add consistent treatment of pressure from external forces on outlets.

jorgensd · 2023-03-23T19:05:38Z

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 · 2023-03-24T11:57:15Z

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]

jorgensd added 9 commits November 9, 2022 16:23

Restructure solver to use IBP on pressure term of tentative velocity eq

cbc7cef

Add some potential boundary conditions

49ef0c8

Various updates

7e5fc6e

Add fix for non-zero pressure

2cb8bd5

Start working on taylor-green. still debugging

7d90a44

Various updates to make project easier to use

df7580a

Various updates

6ac9ef1

Fix ci

98d3bf2

Update test

413c539

jorgensd added 4 commits March 24, 2023 09:00

Debug variant

260d4c2

Cleanup

e1953b6

More adjustments

29f902b

Rewrite demo for convergence study

369394e

jorgensd added 6 commits March 24, 2023 12:35

Fix pressure convergence

49d19cc

Flake8 and mypy

5a5c0b5

Try fixing docs

b5322fb

Fix indentation

5f9b92f

Typing of demos

9ff63b0

Fixes to notebooks

82c6b0c

jorgensd merged commit cd7bb72 into main Mar 24, 2023

jorgensd deleted the dokken/reorganize branch July 4, 2023 11:32

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Add Pressure boundary condition #7

Add Pressure boundary condition #7

jorgensd commented Nov 14, 2022

jorgensd commented Mar 23, 2023

jorgensd commented Mar 24, 2023

Add Pressure boundary condition #7

Add Pressure boundary condition #7

Conversation

jorgensd commented Nov 14, 2022

jorgensd commented Mar 23, 2023

jorgensd commented Mar 24, 2023