Laminar Flow in stenotic U-Bend

This notebook demonstrates the functionality of the trajectory module cmrsim.trajectory.FlowTrajectoryModule, alongside the Re-seeding dataset to simulate particle trajectories inside a stenotic U-Bend. The geometry and the pre-simulated velocity field is given as example resource as a vtk-Mesh.

Imports

# Load TF and check for GPUs
import os
os.environ['TF_CPP_MIN_LOG_LEVEL'] = '2'
import tensorflow as tf
gpus = tf.config.list_physical_devices('GPU')
for gpu in gpus:
    tf.config.experimental.set_memory_growth(gpu, True)
logical_gpus = tf.config.list_logical_devices('GPU')
print(len(gpus), "Physical GPUs,", len(logical_gpus), "Logical GPUs")

# 3rd Party dependencies
from tqdm.notebook import tqdm
from pint import Quantity
import pyvista
import numpy as np
import matplotlib.pyplot as plt
from copy import deepcopy
from IPython.display import display, HTML
import base64
%matplotlib inline

# Project library cmrsim
import sys
sys.path.append("../../")
sys.path.append("../")
import local_functions



sys.path.append("../../../cmrseq")
import cmrseq
import cmrsim
import cmrsim.utils.particle_properties as part_factory

1 Physical GPUs, 1 Logical GPUs

Define default settings for rendering the mesh

custom_camera_pos = [(0.211, 0.267, 0.081),  (0.0262, 0.0261, 0.0172), (0.802, -0.570, -0.178)]
custom_theme = pyvista.themes.DefaultTheme()
custom_theme.background = 'white'
custom_theme.font.family = 'times'
custom_theme.font.color = 'black'
custom_theme.outline_color = "white"
custom_theme.edge_color = 'white'
custom_theme.transparent_background = True

Load Mesh file

The mesh is loaded and geometrically transformed/rotated such that the stanchions of the U-Bend points in z-direction. The modifications are performed using the pyvista package.

# Load UBend Example.
mesh = pyvista.read("../example_resources/sUbend_219_9.vtu")

# The velocity [m/s] is stored in the mesh as U1. In cmrsim all definitions are done in [ms] therefore
# We add the velocity field as scaled data to the mesh
mesh["velocity"] = Quantity(mesh["U"], "m/s").m_as("m/ms")

# To match the usual geometry (Head <-> Foot along z-axis) flip the axes of the mesh
mesh.points = mesh.points[:,[2,1,0]]
mesh["velocity"] = mesh["velocity"][:,[2,1,0]]
global_mesh_offset = [0, 0.04, -0.05]
mesh.translate(global_mesh_offset, inplace=True)
mesh = mesh.cell_data_to_point_data()

# clip to restricted length to reduce mesh-size
mesh.clip(normal='-z', origin=[0., 0., -0.1], inplace=True)

Header	Data Arrays
UnstructuredGrid Information N Cells 1034880 N Points 1061341 X Bounds -1.603e-02, 1.602e-02 Y Bounds -6.312e-02, 5.603e-02 Z Bounds -1.000e-01, 9.140e-02 N Arrays 4	Name Field Type N Comp Min Max U Points float64 3 -2.534e+00 5.380e+00 p Points float64 1 -3.984e+00 9.466e+00 RST Points float64 6 -9.636e-01 2.484e+00 velocity Points float64 3 -2.534e-03 5.380e-03

Render the mesh and display image.

pyvista.close_all()
plotter = pyvista.Plotter(off_screen=True, window_size=(500, 500), theme=custom_theme)
local_functions.add_custom_axes(plotter)
plotter.camera_position = custom_camera_pos
plotter.add_mesh(mesh, scalars="velocity", cmap="twilight", opacity=0.5,
                 scalar_bar_args={"title":"Velocity [m/ms]"}, clim=[0., 0.006])
plotter.add_mesh(mesh.slice_orthogonal(), show_scalar_bar=False, cmap="twilight", clim=[0., 0.007])
img = plotter.screenshot("mesh.png")
b64 = base64.b64encode(open("mesh.png",'rb').read()).decode('ascii')
display(HTML(f'<img src="data:image/gif;base64,{b64}" />'))

Instantiate the Trajectory module

The trajectory module needs a velocity field definition on a regular grid, as the look up and trilinear interpolation is implemented using this assumption. The most convenient way to achieve this, is to use the RefillingFlowDataset class. For more detailed instructions checkout the corresponding notebook and the API reference.

# Define Reseeding slice parameters
seeding_slice_position = Quantity([0, 0.048, -0.045], "m")
seeding_slice_normal = np.array([1., 0., 0.])
seeding_slice_thickness = Quantity(50, "mm")
seeding_pixel_spacing = np.array((0.25e-3, 0.25e-3, 0.25e-3))  # resolution of seeding mesh, mm
lookup_mesh_spacing = np.array((0.5e-3, 0.5e-3, 0.5e-3)) # resolution of lookup mesh, mm
seeding_slice_bb = Quantity([5, 5.1, 2], "cm")
seeding_particle_density = 10 # per 1/mm*3

particle_creation_fncs = {"magnetization": part_factory.norm_magnetization(),
                          "T1": part_factory.uniform(default_value=300.),
                          "T2": part_factory.uniform(default_value=100.),
                          "M0": part_factory.uniform(default_value=1)}

dataset = cmrsim.datasets.RefillingFlowDataset(mesh,
                                               particle_creation_callables=particle_creation_fncs,
                                               slice_position=seeding_slice_position.m_as("m"),
                                               slice_normal=seeding_slice_normal,
                                               slice_thickness=seeding_slice_thickness.m_as("m"),
                                               slice_bounding_box=seeding_slice_bb.m_as("m"),
                                               lookup_map_spacing=lookup_mesh_spacing,
                                               seeding_vol_spacing=seeding_pixel_spacing,
                                               field_list=[("velocity", 3), ])

Updated Mesh
Updated Slice Position

Now actually instantiate the trajectory module by using the Dataset-field-definitions

velocity_field_3d, map_dimensions = dataset.get_lookup_table()

# Setup trajectory module
trajectory_module = cmrsim.datasets.FlowTrajectory(velocity_field_3d,
                                                   map_dimensions=map_dimensions.astype(np.float32),
                                                   additional_dimension_mapping=dataset.field_list[1:],
                                                   device="GPU:0")

Simulate Particle Trajectories

To simulate actual trajetories, we first need to seed particles into the velocity fields. A reasonable assumption is, that the original volume is filled with a uniform density of partilces. Furthermore, in-flow and out-flow is taken care of by calling the RefillingFlowDataset instance every 5 milliseconds, mimicking a sequence with a corresponding TR.

particle_density = Quantity(1, "1/mm^3")
positions, properties = dataset.initial_filling(particle_density.m)
print(f"Filling resulted in a total of: {positions.shape[0]} particles")

Filling resulted in a total of: 257957 particles

Now numerically solve the kinetic equation for particles in a laminar flow-field by calling the increment_particles function in a loop to advance the position of all seeded particles and interleave the reseeding as described above.

For illustration the intermediate positions and velocities are saved every 0.2 milliseconds.

n_TR = 2
steps = 25
substeps = 70
dt = tf.constant(Quantity(10, "us").m_as("ms"))
r = positions
props = properties

trajec, vel = [], []
for tr in tqdm(range(n_TR)):
    pos_per_step = np.empty([steps, *r.shape], dtype=np.float32)
    vel_per_step = np.empty([steps, *r.shape], dtype=np.float32)

    for step in tqdm(range(steps), leave=False):
        for _ in range(substeps):
            r, fields = trajectory_module.increment_particles(r, dt, return_velocities=True)
        pos_per_step[step] = r.numpy()
        vel_per_step[step] = fields["velocity"].numpy()

    trajec.append(pos_per_step), vel.append(vel_per_step)
    r, props, _ = dataset(particle_density.m, residual_particle_pos=r.numpy(),
                               particle_properties=props, reseed_threshold=0.3,
                               distance_tolerance=Quantity(30, "cm").m_as("m"))

Animate Trajectories

import importlib
importlib.reload(local_functions)

timing = np.arange(0, steps * n_TR) * (substeps * dt.numpy())

plotter = pyvista.Plotter(off_screen=True, window_size=(2400, 950), theme=custom_theme, shape=(1, 2))
plotter.subplot(0, 0)
local_functions.add_custom_axes(plotter)
plotter.camera_position = custom_camera_pos
plotter.add_mesh(mesh, scalars="velocity", cmap="twilight", opacity=0.5,
                 scalar_bar_args={"title":"Velocity [m/ms]"}, clim=[0., 0.006])
plotter.add_mesh(mesh.slice_orthogonal(), show_scalar_bar=False, cmap="twilight", clim=[0., 0.007])
local_functions.add_refilling_box(plotter, seeding_slice_bb, seeding_slice_thickness, seeding_slice_position)
plotter.add_text("C", position="upper_left")

plotter.subplot(0, 1)
local_functions.add_custom_axes(plotter)
plotter.camera_position = custom_camera_pos
plotter.add_text("D", position="upper_left")

local_functions.animate_trajectories(plotter, [_[0:] for t in trajec for _ in t], timing,
                                     filename="animation.gif",
                                     scalars=[_[0:] for v in vel for _ in np.linalg.norm(v, axis=-1)],
                                     mesh_kwargs=dict(cmap="twilight", opacity=0.5, clim=[0, 0.007],
                                                      scalar_bar_args={"title":"Velocity [m/ms]"}),
                                     text_kwargs=dict(color="black", position='upper_right'))

plotter.close()

b64 = base64.b64encode(open("animation.gif",'rb').read()).decode('ascii')
display(HTML(f'<img src="data:image/gif;base64,{b64}" />'))