3D navier stokes on 600k cells mesh

General
1

Hi all,

I’m trying to run a CFD simulation using the Chorin’s projection method (from the demo) on a 600k-cell mesh.

I’m running into issues that look like memory issues.

I tried running in parrallel (8 cores) but it also kills the process:

Succesfully load mesh with 76465 cells
Succesfully load mesh with 76443 cells
Succesfully load mesh with 75092 cells
Succesfully load mesh with 78020 cells
Succesfully load mesh with 76799 cells
Succesfully load mesh with 75769 cells
Succesfully load mesh with 75959 cells
Succesfully load mesh with 77101 cells
Solving Navier-Stokes:   0%|          | 0/10 [00:00<?, ?it/s]
===================================================================================
=   BAD TERMINATION OF ONE OF YOUR APPLICATION PROCESSES
=   PID 526610 RUNNING AT lap-mathurin-01
=   EXIT CODE: 9
=   CLEANING UP REMAINING PROCESSES
=   YOU CAN IGNORE THE BELOW CLEANUP MESSAGES
===================================================================================
YOUR APPLICATION TERMINATED WITH THE EXIT STRING: Killed (signal 9)
This typically refers to a problem with your application.
Please see the FAQ page for debugging suggestions

Any pointers regarding this? Is this too big of a mesh? I’ve read application cases running on millions of cells so I’m a bit surprised.

Thanks for the help!

Rem

2

Without actually seeing the code, it is quite impossible to give any guidance.

It is also not clear from your post what version of dolfin you use.

Personally Ive ran Navier Stokes simulations with dolfinx with more than 4 millon cells, and general elasticity problems with 50 million cells with no issues.

3

Sure. This is the script I have. I didn’t want to share the large script if it was obvious that it wasn’t possible.

I’m using legacy fenics but would be keen on using dolfinx if it works!

This is mesh that I use.

from fenics import *
import fenics as fe
import numpy as np

# writes velocity field to an xdmf so that doesn't have to solve for it each time in the festim simulation

fluid_id = 6
inlet_id = 7
outlet_id = 8
vacuum_id = 9
walls_id = 10


def rho_lipb(T):  # units in kg/(m**3)
    return 10520.35 - 1.19051 * T


def visc_lipb(T):  # units (Pa s)
    return (
        0.01555147189 - 4.827051855e-05 * T + 5.641475215e-08 * T**2 - 2.2887e-11 * T**3
    )


def fluid_dynamics_sim_chorin(
    mesh, volume_markers, surface_markers, id_inlet, id_outlet, id_walls
):
    """
    Solves the Navier-Stokes equations using the Chorin's projection method
    See https://fenicsproject.org/pub/tutorial/html/._ftut1009.html#ftut1:NS for more details

    Args:
        mesh (fenics.Mesh): the mesh
        volume_markers (fenics.MeshFunction): the volume markers (subdomains)
        surface_markers (fenics.MeshFunction): the surface markers (boundaries)
        id_inlet (int): the id of the inlet boundary
        id_outlet (int): the id of the outlet boundary
        id_walls (int or list): the id(s) of the walls

    Returns:
        fenics.Function: the velocity field
        fenics.Function: the pressure field
    """
    V = fe.VectorFunctionSpace(mesh, "CG", 2)
    Q = fe.FunctionSpace(mesh, "CG", 1)

    u = fe.TrialFunction(V)
    p = fe.TrialFunction(Q)
    v = fe.TestFunction(V)
    q = fe.TestFunction(Q)

    u_n = fe.Function(V)
    u_ = fe.Function(V)
    p_n = fe.Function(Q)
    p_ = fe.Function(Q)

    # ##### Boundary conditions ##### #

    inlet_velocity = 2e-03  # units: m s-1
    outlet_pressure = 0  # units: Pa

    inflow = fe.DirichletBC(
        V, fe.Constant((inlet_velocity, 0.0, 0.0)), surface_markers, id_inlet
    )

    # make sure id_walls is a list
    if isinstance(id_walls, int):
        id_walls = [id_walls]

    # iterate through the walls
    walls = []
    for id_wall in id_walls:
        walls.append(
            fe.DirichletBC(V, fe.Constant((0.0, 0.0, 0.0)), surface_markers, id_wall)
        )

    pressure_outlet = fe.DirichletBC(
        Q, fe.Constant(outlet_pressure), surface_markers, id_outlet
    )
    bcu = [inflow] + walls
    bcp = [pressure_outlet]

    # ##### Solver ##### #
    dx = fe.Measure("dx", subdomain_data=volume_markers)
    ds = fe.Measure("ds", subdomain_data=surface_markers)

    t = 0
    total_time = 20
    dt = 0.1  # Time step size
    num_steps = int(total_time / dt)
    num_steps = 10

    k = fe.Constant(dt)
    n = fe.FacetNormal(mesh)
    U = 0.5 * (u_n + u)

    # LiPb
    T = 700  # TODO make this more generic
    mu = visc_lipb(T)
    rho = rho_lipb(T)

    def epsilon(u):
        return fe.sym(fe.nabla_grad(u))

    def sigma(u, p):
        return 2 * mu * epsilon(u) - p * fe.Identity(len(u))

    # ##### Solver ##### #

    # Tentative velocity step
    F1 = rho * fe.dot((u - u_n) / k, v) * dx
    F1 += rho * fe.dot(fe.dot(u_n, fe.nabla_grad(u_n)), v) * dx
    F1 += fe.inner(sigma(U, p_n), epsilon(v)) * dx
    F1 += fe.dot(p_n * n, v) * ds
    F1 -= fe.dot(mu * fe.nabla_grad(U) * n, v) * ds

    a1 = fe.lhs(F1)
    L1 = fe.rhs(F1)

    solver1 = fe.KrylovSolver("bicgstab", "hypre_amg")
    solver1.parameters["absolute_tolerance"] = 1e-08
    solver1.parameters["relative_tolerance"] = 1e-08
    solver1.parameters["maximum_iterations"] = 1000
    solver1.parameters["report"] = False
    solver1.parameters["monitor_convergence"] = False

    # Pressure update
    a2 = fe.dot(fe.nabla_grad(p), fe.nabla_grad(q)) * dx
    L2 = fe.dot(fe.nabla_grad(p_n), fe.nabla_grad(q)) * dx
    L2 -= (1 / k) * fe.div(u_) * q * dx

    solver2 = fe.KrylovSolver("bicgstab", "hypre_amg")
    solver2.parameters["absolute_tolerance"] = 1e-08
    solver2.parameters["relative_tolerance"] = 1e-08
    solver2.parameters["maximum_iterations"] = 1000
    solver2.parameters["report"] = False
    solver2.parameters["monitor_convergence"] = False

    # Velocity update
    a3 = fe.dot(u, v) * dx
    L3 = fe.dot(u_, v) * dx
    L3 -= k * fe.dot(fe.nabla_grad(p_ - p_n), v) * dx

    solver3 = fe.KrylovSolver("cg", "sor")
    solver3.parameters["absolute_tolerance"] = 1e-08
    solver3.parameters["relative_tolerance"] = 1e-08
    solver3.parameters["maximum_iterations"] = 1000
    solver3.parameters["report"] = False
    solver3.parameters["monitor_convergence"] = False

    # Assemble matrices
    A1 = fe.assemble(a1)
    A2 = fe.assemble(a2)
    A3 = fe.assemble(a3)

    [bc.apply(A1) for bc in bcu]
    [bc.apply(A2) for bc in bcp]

    velocity_file = fe.XDMFFile("velocity_field.xdmf")
    pressure_file = fe.XDMFFile("pressure_field.xdmf")

    max_u = []

    # Time-stepping
    append = False
    for i in range(num_steps):
        # Update current time step
        t += dt

        # Compute tentative velocity step
        b1 = fe.assemble(L1)
        [bc.apply(A1, b1) for bc in bcu]
        solver1.solve(A1, u_.vector(), b1)

        # Pressure correction
        b2 = fe.assemble(L2)
        [bc.apply(A2, b2) for bc in bcp]
        solver2.solve(A2, p_.vector(), b2)

        # Velocity correction
        b3 = fe.assemble(L3)
        [bc.apply(A3, b3) for bc in bcu]
        solver3.solve(A3, u_.vector(), b3)

        # Move to next time step
        u_n.assign(u_)
        p_n.assign(p_)

        max_u.append(u_.vector().max())
        np.savetxt("3D_case_max_u.txt", np.array(max_u))

        velocity_file.write_checkpoint(
            u_,
            "velocity",
            t,
            XDMFFile.Encoding.HDF5,
            append=append,
        )
        pressure_file.write_checkpoint(
            p_,
            "pressure",
            t,
            XDMFFile.Encoding.HDF5,
            append=append,
        )
        append = True
    return u_, p_


if __name__ == "__main__":

    volume_file = "mesh_cells.xdmf"
    boundary_file = "mesh_boundaries.xdmf"
    mesh = fe.Mesh()
    fe.XDMFFile(volume_file).read(mesh)

    # Read tags for volume elements
    volume_markers = fe.MeshFunction("size_t", mesh, mesh.topology().dim())
    fe.XDMFFile(volume_file).read(volume_markers)

    # Read tags for surface elements
    # (can also be used for applying DirichletBC)
    surface_markers = fe.MeshValueCollection("size_t", mesh, mesh.topology().dim() - 1)
    fe.XDMFFile(boundary_file).read(surface_markers, "f")
    surface_markers = fe.MeshFunction("size_t", mesh, surface_markers)

    print("Succesfully load mesh with " + str(len(volume_markers)) + " cells")

    fluid_dynamics_sim_chorin(
        mesh=mesh,
        volume_markers=volume_markers,
        surface_markers=surface_markers,
        id_inlet=inlet_id,
        id_outlet=outlet_id,
        id_walls=[walls_id, vacuum_id],
    )

I

4

Have you tried to locate what line of thai script segfaults on 8 procs? It seems like it is happening prior to calling solve, thus alot of code could be removed from the snippet above.

5

Yep, I’ve reduced the code to this and the solve seems to be the culprit

import fenics as fe
import numpy as np

# writes velocity field to an xdmf so that doesn't have to solve for it each time in the festim simulation

fluid_id = 6
inlet_id = 7
outlet_id = 8
vacuum_id = 9
walls_id = 10


def fluid_dynamics_sim_chorin(
    mesh, volume_markers, surface_markers, id_inlet, id_outlet, id_walls
):
    """
    Solves the Navier-Stokes equations using the Chorin's projection method
    See https://fenicsproject.org/pub/tutorial/html/._ftut1009.html#ftut1:NS for more details

    Args:
        mesh (fenics.Mesh): the mesh
        volume_markers (fenics.MeshFunction): the volume markers (subdomains)
        surface_markers (fenics.MeshFunction): the surface markers (boundaries)
        id_inlet (int): the id of the inlet boundary
        id_outlet (int): the id of the outlet boundary
        id_walls (int or list): the id(s) of the walls

    Returns:
        fenics.Function: the velocity field
        fenics.Function: the pressure field
    """
    V = fe.VectorFunctionSpace(mesh, "CG", 2)
    Q = fe.FunctionSpace(mesh, "CG", 1)

    u = fe.TrialFunction(V)
    p = fe.TrialFunction(Q)
    v = fe.TestFunction(V)
    q = fe.TestFunction(Q)

    u_n = fe.Function(V)
    u_ = fe.Function(V)
    p_n = fe.Function(Q)
    p_ = fe.Function(Q)

    # ##### Boundary conditions ##### #

    inlet_velocity = 2e-03  # units: m s-1
    outlet_pressure = 0  # units: Pa

    inflow = fe.DirichletBC(
        V, fe.Constant((inlet_velocity, 0.0, 0.0)), surface_markers, id_inlet
    )

    # make sure id_walls is a list
    if isinstance(id_walls, int):
        id_walls = [id_walls]

    # iterate through the walls
    walls = []
    for id_wall in id_walls:
        walls.append(
            fe.DirichletBC(V, fe.Constant((0.0, 0.0, 0.0)), surface_markers, id_wall)
        )

    pressure_outlet = fe.DirichletBC(
        Q, fe.Constant(outlet_pressure), surface_markers, id_outlet
    )
    bcu = [inflow] + walls
    bcp = [pressure_outlet]

    # ##### Solver ##### #
    dx = fe.Measure("dx", subdomain_data=volume_markers)
    ds = fe.Measure("ds", subdomain_data=surface_markers)

    dt = 0.1  # Time step size

    k = fe.Constant(dt)
    n = fe.FacetNormal(mesh)
    U = 0.5 * (u_n + u)

    # LiPb
    mu = 1
    rho = 1

    def epsilon(u):
        return fe.sym(fe.nabla_grad(u))

    def sigma(u, p):
        return 2 * mu * epsilon(u) - p * fe.Identity(len(u))

    # ##### Solver ##### #

    # Tentative velocity step
    F1 = rho * fe.dot((u - u_n) / k, v) * dx
    F1 += rho * fe.dot(fe.dot(u_n, fe.nabla_grad(u_n)), v) * dx
    F1 += fe.inner(sigma(U, p_n), epsilon(v)) * dx
    F1 += fe.dot(p_n * n, v) * ds
    F1 -= fe.dot(mu * fe.nabla_grad(U) * n, v) * ds

    a1 = fe.lhs(F1)
    L1 = fe.rhs(F1)

    solver1 = fe.KrylovSolver("bicgstab", "hypre_amg")
    solver1.parameters["absolute_tolerance"] = 1e-08
    solver1.parameters["relative_tolerance"] = 1e-08
    solver1.parameters["maximum_iterations"] = 1000
    solver1.parameters["report"] = False
    solver1.parameters["monitor_convergence"] = False

    # Pressure update
    a2 = fe.dot(fe.nabla_grad(p), fe.nabla_grad(q)) * dx
    L2 = fe.dot(fe.nabla_grad(p_n), fe.nabla_grad(q)) * dx
    L2 -= (1 / k) * fe.div(u_) * q * dx

    solver2 = fe.KrylovSolver("bicgstab", "hypre_amg")
    solver2.parameters["absolute_tolerance"] = 1e-08
    solver2.parameters["relative_tolerance"] = 1e-08
    solver2.parameters["maximum_iterations"] = 1000
    solver2.parameters["report"] = False
    solver2.parameters["monitor_convergence"] = False

    # Velocity update
    a3 = fe.dot(u, v) * dx
    L3 = fe.dot(u_, v) * dx
    L3 -= k * fe.dot(fe.nabla_grad(p_ - p_n), v) * dx

    solver3 = fe.KrylovSolver("cg", "sor")
    solver3.parameters["absolute_tolerance"] = 1e-08
    solver3.parameters["relative_tolerance"] = 1e-08
    solver3.parameters["maximum_iterations"] = 1000
    solver3.parameters["report"] = False
    solver3.parameters["monitor_convergence"] = False

    # Assemble matrices
    A1 = fe.assemble(a1)
    A2 = fe.assemble(a2)
    A3 = fe.assemble(a3)

    [bc.apply(A1) for bc in bcu]
    [bc.apply(A2) for bc in bcp]

    # Compute tentative velocity step
    b1 = fe.assemble(L1)
    [bc.apply(A1, b1) for bc in bcu]
    solver1.solve(A1, u_.vector(), b1)  # <-- culprit!

    # # Pressure correction
    # b2 = fe.assemble(L2)
    # [bc.apply(A2, b2) for bc in bcp]
    # solver2.solve(A2, p_.vector(), b2)

    # # Velocity correction
    # b3 = fe.assemble(L3)
    # [bc.apply(A3, b3) for bc in bcu]
    # solver3.solve(A3, u_.vector(), b3)

    # # Move to next time step
    # u_n.assign(u_)
    # p_n.assign(p_)

    return u_, p_


if __name__ == "__main__":

    volume_file = "mesh_cells.xdmf"
    boundary_file = "mesh_facets.xdmf"
    mesh = fe.Mesh()
    fe.XDMFFile(volume_file).read(mesh)

    # Read tags for volume elements
    volume_markers = fe.MeshFunction("size_t", mesh, mesh.topology().dim())
    fe.XDMFFile(volume_file).read(volume_markers)

    # Read tags for surface elements
    # (can also be used for applying DirichletBC)
    surface_markers = fe.MeshValueCollection("size_t", mesh, mesh.topology().dim() - 1)
    fe.XDMFFile(boundary_file).read(surface_markers, "f")
    surface_markers = fe.MeshFunction("size_t", mesh, surface_markers)

    print("Succesfully load mesh with " + str(len(volume_markers)) + " cells")

    fluid_dynamics_sim_chorin(
        mesh=mesh,
        volume_markers=volume_markers,
        surface_markers=surface_markers,
        id_inlet=inlet_id,
        id_outlet=outlet_id,
        id_walls=[walls_id, vacuum_id],
    )
6

Is there a reason for applying the bc twice?

Have you Also tried using gmres with a Jacobi smoother instead of AMG?

1 Like
7

No reason for the double BC application. Probably a copy paste error.

I change solver1 to

solver1 = fe.KrylovSolver("gmres", "ilu")

And it seems to work now!
Was it what you had in mind?

8

I was thinking solver1 = fe.KrylovSolver("gmres", "jacobi")

9

This worked well too.

Would you use it for all three solvers or just the first one?

10

My goto options for the three steps of ipcs are

solver_options = {"tentative": {"ksp_type": "gmres", "pc_type": "jacobi", "ksp_error_if_not_converged": True}, 
                  "pressure":
                  {"ksp_type": "cg", "pc_type": "hypre", "pc_hypre_type": "boomeramg", "ksp_atol": 1e-9,
                               "ksp_rtol": 1e-9, "ksp_error_if_not_converged": True},
                  "scalar":
                  {"ksp_type": "cg", "pc_type": "jacobi", "ksp_atol": 1e-9,
                      "ksp_rtol": 1e-9, "ksp_error_if_not_converged": True},

                  }
1 Like