STL Geometry: Blood Flow in Intracranial Aneurysm

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

aneurysm.png

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 Sym 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/.

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\).

In this tutorial, you will use Modulus Sym’ Tessellation module to sample points using a STL geometry. The module works similar to Modulus 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 Modulus Sym’ Tessellation module to the sample points on the STL geometry.

Copy
Copied!
            

import os import warnings import torch import numpy as np from sympy import Symbol, sqrt, Max import modulus.sym from modulus.sym.hydra import to_absolute_path, instantiate_arch, ModulusConfig from modulus.sym.solver import Solver from modulus.sym.domain import Domain from modulus.sym.domain.constraint import ( PointwiseBoundaryConstraint, PointwiseInteriorConstraint, IntegralBoundaryConstraint, ) from modulus.sym.domain.validator import PointwiseValidator from modulus.sym.domain.monitor import PointwiseMonitor from modulus.sym.key import Key from modulus.sym.eq.pdes.navier_stokes import NavierStokes from modulus.sym.eq.pdes.basic import NormalDotVec from modulus.sym.utils.io import csv_to_dict from modulus.sym.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.

Copy
Copied!
            

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)

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.

Copy
Copied!
            

# 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.

Copy
Copied!
            

# 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 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.

Copy
Copied!
            

# 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/modulus/resources/modulus_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, ) domain.add_monitor(pressure_monitor)

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

Copy
Copied!
            

python aneurysm.py

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.

val_errors.png

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 Sym successfully captures the flow field inside.

aneurysm_v_mag_labelled.png

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

aneurysm_p_labelled.png

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

aneurysm_streamlines3_crop.png

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

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.

aneurysm_transfer_learning.png

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

To use transfer learning in Modulus 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.

Copy
Copied!
            

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

Previous Interface Problem by Variational Method
Next Moving Time Window: Taylor Green Vortex Decay
© Copyright 2023, NVIDIA Modulus Team. Last updated on Jan 25, 2024.