STL Geometry: Blood Flow in Intracranial Aneurysm#

Introduction#

In this tutorial, you will import an STL file for a complicated geometry and use PhysicsNeMo Sym’ 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 PhysicsNeMo Sym and sample points in the interior and on the surface of the geometry.

Aneurysm STL file

Fig. 141 Aneurysm STL file#

Note

This tutorial assumes that you have completed tutorial Introductory Example and have familiarized yourself with the basics of the PhysicsNeMo Sym APIs.

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 PhysicsNeMo Sym’ Tessellation module to sample points using a STL geometry. The module works similar to PhysicsNeMo Sym’ 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 PhysicsNeMo Sym’ Tessellation module to the sample points on the STL geometry.

# limitations under the License.

import os
import warnings

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

import physicsnemo.sym
from physicsnemo.sym.hydra import to_absolute_path, instantiate_arch, PhysicsNeMoConfig
from physicsnemo.sym.solver import Solver
from physicsnemo.sym.domain import Domain
from physicsnemo.sym.domain.constraint import (
    PointwiseBoundaryConstraint,
    PointwiseInteriorConstraint,
    IntegralBoundaryConstraint,
)
from physicsnemo.sym.domain.validator import PointwiseValidator
from physicsnemo.sym.domain.monitor import PointwiseMonitor
from physicsnemo.sym.key import Key
from physicsnemo.sym.eq.pdes.navier_stokes import NavierStokes
from physicsnemo.sym.eq.pdes.basic import NormalDotVec

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.

def run(cfg: PhysicsNeMoConfig) -> None:
    # 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)

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

    # 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},
    )

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 PhysicsNeMo Sym 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
    file_path = "./openfoam/aneurysm_parabolicInlet_sol0.csv"
    if os.path.exists(to_absolute_path(file_path)):
        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(file_path), 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)
    else:
        warnings.warn(
            f"Directory {file_path} does not exist. Will skip adding validators. Please download the additional files from NGC https://catalog.ngc.nvidia.com/orgs/nvidia/teams/physicsnemo/resources/physicsnemo_sym_examples_supplemental_materials"
        )

    # 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,
    )

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. 142 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. 142 Convergence plots for different point density#

Fig. 144 shows the pressure developed inside the aneurysm and the vein. A cross-sectional view in Fig. 143 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. 145 shows that PhysicsNeMo Sym successfully captures the flow field inside.

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

Fig. 143 Cross-sectional view aneurysm showing velocity magnitude. Left: PhysicsNeMo Sym. Center: OpenFOAM. Right: Difference#

Pressure across aneurysm. Left: PhysicsNeMo Sym. Center: OpenFOAM. Right: Difference

Fig. 144 Pressure across aneurysm. Left: PhysicsNeMo Sym. Center: OpenFOAM. Right: Difference#

Flow streamlines inside the aneurysm generated from PhysicsNeMo Sym simulation.

Fig. 145 Flow streamlines inside the aneurysm generated from PhysicsNeMo Sym 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. 146 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. 146 Transfer learning accelerates intracranial aneurysm simulations. Results are for two intracranial aneurysms with different sac shapes.#

To use transfer learning in PhysicsNeMo Sym, 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 :
  - physicsnemo_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