STL Geometry: Blood Flow in Intracranial Aneurysm

Introduction

In this tutorial, you will import an STL file for a complicated geometry and use Modulus’ SDF library to sample points on the surface and the interior of the STL and train the PINNs to predict flow in this complex geometry. In this tutorial you will learn the following:

  1. How to import an STL file in Modulus and sample points in the interior and on the surface of the geometry.

Aneurysm STL file

Fig. 108 Aneurysm STL file

Note

This tutorial assumes that you have completed tutorial Introductory Example and have familiarized yourself with the basics of the Modulus APIs. Additionally, to use the modules described in this tutorial, make sure your system satisfies the requirements for SDF library (System Requirements).

For the interior sampling to work, ensure that the STL geometry is watertight. This requirement is not necessary for sampling points on the surface.

All the python scripts for this problem can be found at examples/aneurysm/.

Problem Description

This simulation, uses a no-slip boundary condition on the walls of the aneurysm \(u,v,w=0\). For the inlet, a parabolic flow where the flow goes in the normal direction of the inlet and has peak velocity 1.5, is used. The outlet has a zero pressure condition, \(p=0\). The kinematic viscosity of the fluid is \(0.025\) and the density is a constant \(1.0\).

Case Setup

In this tutorial, you will use Modulus’ Tessellation module to sample points using a STL geometry. The module works similar to Modulus’ geometry module. Which means you can use PointwiseInteriorConstraint and PointwiseBoundaryConstraint to sample points in the interior and the boundary of the geometry and define appropriate constraints. Separate STL files for each boundary of the geometry and another watertight geometry for sampling points in the interior of the geometry are required.

Importing the required packages

The list of required packages can be found below. Import Modulus’ Tessellation module to the sample points on the STL geometry.

import torch
import numpy as np
from sympy import Symbol, sqrt, Max

import modulus
from modulus.hydra import to_absolute_path, instantiate_arch, ModulusConfig
from modulus.solver import Solver
from modulus.domain import Domain
from modulus.domain.constraint import (
    PointwiseBoundaryConstraint,
    PointwiseInteriorConstraint,
    IntegralBoundaryConstraint,
)
from modulus.domain.validator import PointwiseValidator
from modulus.domain.monitor import PointwiseMonitor
from modulus.key import Key
from modulus.eq.pdes.navier_stokes import NavierStokes
from modulus.eq.pdes.basic import NormalDotVec
from modulus.utils.io import csv_to_dict
from modulus.geometry.tessellation import Tessellation

Using STL files to generate point clouds

Import the STL geometries using the Tessellation.from_stl() function. This function takes in the path of the STL geometry as input. You will need to specify the value of attribute airtight as False for the open surfaces (eg. boundary STL files).

Then these mesh objects can be used to create boundary or interior constraints similar to tutorial Introductory Example using the PointwiseBoundaryConstraint or PointwiseInteriorConstraint.

Note

For this tutorial, you can normalize the geometry by scaling it and centering it about the origin (0, 0, 0). This will help in speeding up the training process.

The code to sample using STL geometry, define all these functions, boundary conditions is shown below.

    # read stl files to make geometry
    point_path = to_absolute_path("./stl_files")
    inlet_mesh = Tessellation.from_stl(
        point_path + "/aneurysm_inlet.stl", airtight=False
    )
    outlet_mesh = Tessellation.from_stl(
        point_path + "/aneurysm_outlet.stl", airtight=False
    )
    noslip_mesh = Tessellation.from_stl(
        point_path + "/aneurysm_noslip.stl", airtight=False
    )
    integral_mesh = Tessellation.from_stl(
        point_path + "/aneurysm_integral.stl", airtight=False
    )
    interior_mesh = Tessellation.from_stl(
        point_path + "/aneurysm_closed.stl", airtight=True
    )

    # params
    nu = 0.025
    inlet_vel = 1.5

    # inlet velocity profile
    def circular_parabola(x, y, z, center, normal, radius, max_vel):
        centered_x = x - center[0]
        centered_y = y - center[1]
        centered_z = z - center[2]
        distance = sqrt(centered_x ** 2 + centered_y ** 2 + centered_z ** 2)
        parabola = max_vel * Max((1 - (distance / radius) ** 2), 0)
        return normal[0] * parabola, normal[1] * parabola, normal[2] * parabola

    # normalize meshes
    def normalize_mesh(mesh, center, scale):
        mesh = mesh.translate([-c for c in center])
        mesh = mesh.scale(scale)
        return mesh

    # normalize invars
    def normalize_invar(invar, center, scale, dims=2):
        invar["x"] -= center[0]
        invar["y"] -= center[1]
        invar["z"] -= center[2]
        invar["x"] *= scale
        invar["y"] *= scale
        invar["z"] *= scale
        if "area" in invar.keys():
            invar["area"] *= scale ** dims
        return invar

    # scale and normalize mesh and openfoam data
    center = (-18.40381048596882, -50.285383353981196, 12.848136936899031)
    scale = 0.4
    inlet_mesh = normalize_mesh(inlet_mesh, center, scale)
    outlet_mesh = normalize_mesh(outlet_mesh, center, scale)
    noslip_mesh = normalize_mesh(noslip_mesh, center, scale)
    integral_mesh = normalize_mesh(integral_mesh, center, scale)
    interior_mesh = normalize_mesh(interior_mesh, center, scale)

    # geom params
    inlet_normal = (0.8526, -0.428, 0.299)
    inlet_area = 21.1284 * (scale ** 2)
    inlet_center = (-4.24298030045776, 4.082857101816247, -4.637790193399717)
    inlet_radius = np.sqrt(inlet_area / np.pi)
    outlet_normal = (0.33179, 0.43424, 0.83747)
    outlet_area = 12.0773 * (scale ** 2)

Defining the Equations, Networks and Nodes

This process is similar to other tutorials. In this problem you are only solving for laminar flow, so you can use only NavierStokes and NormalDotVec equations and define a network similar to tutorial Introductory Example. The code to generate the Network and required nodes is shown below.

    # make list of nodes to unroll graph on
    ns = NavierStokes(nu=nu * scale, rho=1.0, dim=3, time=False)
    normal_dot_vel = NormalDotVec(["u", "v", "w"])
    flow_net = instantiate_arch(
        input_keys=[Key("x"), Key("y"), Key("z")],
        output_keys=[Key("u"), Key("v"), Key("w"), Key("p")],
        cfg=cfg.arch.fully_connected,
    )
    nodes = (
        ns.make_nodes()
        + normal_dot_vel.make_nodes()
        + [flow_net.make_node(name="flow_network")]
    )

Setting up Domain and adding Constraints

Now that you have all the nodes and geometry elements defined, you can use the tesselated/mesh objects to create boundary or interior constraints similar to tutorial Introductory Example using the PointwiseBoundaryConstraint or PointwiseInteriorConstraint.

    # make aneurysm domain
    domain = Domain()
    # add constraints to solver
    # inlet
    u, v, w = circular_parabola(
        Symbol("x"),
        Symbol("y"),
        Symbol("z"),
        center=inlet_center,
        normal=inlet_normal,
        radius=inlet_radius,
        max_vel=inlet_vel,
    )
    inlet = PointwiseBoundaryConstraint(
        nodes=nodes,
        geometry=inlet_mesh,
        outvar={"u": u, "v": v, "w": w},
        batch_size=cfg.batch_size.inlet,
    )
    domain.add_constraint(inlet, "inlet")

    # outlet
    outlet = PointwiseBoundaryConstraint(
        nodes=nodes,
        geometry=outlet_mesh,
        outvar={"p": 0},
        batch_size=cfg.batch_size.outlet,
    )
    domain.add_constraint(outlet, "outlet")

    # no slip
    no_slip = PointwiseBoundaryConstraint(
        nodes=nodes,
        geometry=noslip_mesh,
        outvar={"u": 0, "v": 0, "w": 0},
        batch_size=cfg.batch_size.no_slip,
    )
    domain.add_constraint(no_slip, "no_slip")

    # interior
    interior = PointwiseInteriorConstraint(
        nodes=nodes,
        geometry=interior_mesh,
        outvar={"continuity": 0, "momentum_x": 0, "momentum_y": 0, "momentum_z": 0},
        batch_size=cfg.batch_size.interior,
    )
    domain.add_constraint(interior, "interior")

    # Integral Continuity 1
    integral_continuity = IntegralBoundaryConstraint(
        nodes=nodes,
        geometry=outlet_mesh,
        outvar={"normal_dot_vel": 2.540},
        batch_size=1,
        integral_batch_size=cfg.batch_size.integral_continuity,
        lambda_weighting={"normal_dot_vel": 0.1},
    )
    domain.add_constraint(integral_continuity, "integral_continuity_1")

    # Integral Continuity 2
    integral_continuity = IntegralBoundaryConstraint(
        nodes=nodes,
        geometry=integral_mesh,
        outvar={"normal_dot_vel": -2.540},
        batch_size=1,
        integral_batch_size=cfg.batch_size.integral_continuity,
        lambda_weighting={"normal_dot_vel": 0.1},
    )
    domain.add_constraint(integral_continuity, "integral_continuity_1")

Adding Validators and Monitors

The process of adding validation data and monitors is similar to previous tutorials. This example uses the simulation from OpenFOAM for validating the Modulus results. Also, you can create a monitor for pressure drop across the aneurysm to monitor the convergence and compare against OpenFOAM data. The code to generate the these domains is shown below.

    # add validation data
    mapping = {
        "Points:0": "x",
        "Points:1": "y",
        "Points:2": "z",
        "U:0": "u",
        "U:1": "v",
        "U:2": "w",
        "p": "p",
    }
    openfoam_var = csv_to_dict(
        to_absolute_path("openfoam/aneurysm_parabolicInlet_sol0.csv"), mapping
    )
    openfoam_invar = {
        key: value for key, value in openfoam_var.items() if key in ["x", "y", "z"]
    }
    openfoam_invar = normalize_invar(openfoam_invar, center, scale, dims=3)
    openfoam_outvar = {
        key: value for key, value in openfoam_var.items() if key in ["u", "v", "w", "p"]
    }
    openfoam_validator = PointwiseValidator(
        nodes=nodes, invar=openfoam_invar, true_outvar=openfoam_outvar, batch_size=4096
    )
    domain.add_validator(openfoam_validator)

    # add pressure monitor
    pressure_monitor = PointwiseMonitor(
        inlet_mesh.sample_boundary(16),
        output_names=["p"],
        metrics={"pressure_drop": lambda var: torch.mean(var["p"])},
        nodes=nodes,
    )
    domain.add_monitor(pressure_monitor)

    # make solver
    slv = Solver(cfg, domain)

    # start solver
    slv.solve()


if __name__ == "__main__":
    run()

Training the model

Once the python file is setup, the training can be simply started by executing the python script.

python aneurysm.py

Results and Post-processing

We use this tutorial to give an example of overfitting of training data in the PINNs. Fig. 109 shows the comparison of the validation error plots achieved for two different point densities. The case using 10 M points shows an initial convergence which later diverges even when the training error keeps reducing. This implies that the network is overfitting the sampled points while sacrificing the accuracy of flow in between them. Increasing the points to 20 M solves that problem and the flow field is generalized to a better resolution.

Convergence plots for different point density

Fig. 109 Convergence plots for different point density

Fig. 111 shows the pressure developed inside the aneurysm and the vein. A cross-sectional view in Fig. 110 shows the distribution of velocity magnitude inside the aneurysm. One of the key challenges of this problem is getting the flow to develop inside the aneurysm sac and the streamline plot in Fig. 112 shows that Modulus successfully captures the flow field inside.

Cross-sectional view aneurysm showing velocity magnitude. Left: Modulus. Center: OpenFOAM. Right: Difference

Fig. 110 Cross-sectional view aneurysm showing velocity magnitude. Left: Modulus. Center: OpenFOAM. Right: Difference

Pressure across aneurysm. Left: Modulus. Center: OpenFOAM. Right: Difference

Fig. 111 Pressure across aneurysm. Left: Modulus. Center: OpenFOAM. Right: Difference

Flow streamlines inside the aneurysm generated from Modulus simulation.

Fig. 112 Flow streamlines inside the aneurysm generated from Modulus simulation.

Accelerating the Training of Neural Network Solvers via Transfer Learning

Numerous applications in science and engineering require repetitive simulations, such as simulation of blood flow in different patient specific models. Traditional solvers simulate these models independently and from scratch. Even a minor change to the model geometry (such as an adjustment to the patient specific medical image segmentation) requires a new simulation. Interestingly, and unlike the traditional solvers, neural network solvers can transfer knowledge across different neural network models via transfer learning. In transfer learning, the knowledge acquired by a (source) trained neural network model for a physical system is transferred to another (target) neural network model that is to be trained for a similar physical system with slightly different characteristics (such as geometrical differences). The network parameters of the target model are initialized from the source model, and are retrained to cope with the new system characteristics without having the neural network model trained from scratch. This transfer of knowledge effectively reduces the time to convergence for neural network solvers. As an example, Fig. 113 shows the application of transfer learning in training of neural network solvers for two intracranial aneurysm models with different sac shapes.

Transfer learning accelerates intracranial aneurysm simulations. Results are for two intracranial aneurysms with different sac shapes.

Fig. 113 Transfer learning accelerates intracranial aneurysm simulations. Results are for two intracranial aneurysms with different sac shapes.

To use transfer learning in Modulus, set 'initialize_network_dir' in the configs to the source model network checkpoint. Also, since in transfer learning you fine-tune the source model instead of training from scratch, use a relatively smaller learning rate compared to a full run, with smaller number of iterations and faster decay, as shown below.

defaults :
  - modulus_default
  - arch:
      - fully_connected
  - scheduler: tf_exponential_lr
  - optimizer: adam
  - loss: sum
  - _self_

scheduler:
  decay_rate: 0.95
  #decay_steps: 15000 # full run
  decay_steps: 6000   # TL run

network_dir : "network_checkpoint_target"
initialization_network_dir : "../aneurysm/network_checkpoint_source/"

training:
  rec_results_freq : 10000
  rec_constraint_freq: 50000
  #max_steps: 1500000 # full run
  max_steps: 400000   # TL run

batch_size:
  inlet: 1100
  outlet: 650
  no_slip: 5200
  interior: 6000
  integral_continuity: 310