Hi @npstrike and @ngeneva
I wanted to revisit this interesting problem and share my approach, which might help others tackling similar challenges or serve as a base tutorial for CFD predictions using NVIDIA Modulus. Using the geometry provided by @npstrike, I scaled it down to a box of 1m and attempted to reproduce the flow behavior. To simplify the setup, I used the following parameters:
scale = 1
nu = 0.025
inlet_vel = 0.5 m/s
After just 100k iterations, the results look promising in terms of capturing the expected flow and velocity behavior inside the T-junction. I’m in the process of validating the Modulus output against OpenFOAM simulations. I’ll update this post with the validation results once complete.
import math
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.utils.io.vtk import var_to_polyvtk
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, Curl
from modulus.sym.eq.pdes.basic import NormalDotVec
from modulus.sym.utils.io import csv_to_dict
from modulus.sym.geometry.tessellation import Tessellation
from modulus.sym.eq.pdes.turbulence_zero_eq import ZeroEquation
from modulus.sym.domain.inferencer import VoxelInferencer
@modulus.sym.main(config_path="conf", config_name="config_stubbed")
def run(cfg: ModulusConfig) -> None:
# read stl files to make geometry
point_path = to_absolute_path("./stl_new")
# Base imports (dimensions are in mm)
wall_pipe_mesh = Tessellation.from_stl(point_path + "/WALL_PIPE.stl", airtight=False)
inlet_main_mesh = Tessellation.from_stl(point_path + "/INLET_MAIN.stl", airtight=False)
inlet_branch_mesh = Tessellation.from_stl(point_path + "/INLET_BRANCH.stl", airtight=False)
outlet_mesh = Tessellation.from_stl(point_path + "/OUTLET.stl", airtight=False)
internal_volume_mesh = Tessellation.from_stl(point_path + "/INTERNAL_VOLUME.stl", airtight=True)
integral_main_mesh = Tessellation.from_stl(point_path + "/INTEGRAL_MAIN.stl", airtight=False)
integral_branch_mesh = Tessellation.from_stl(point_path + "/INTEGRAL_BRANCH.stl", airtight=False)
integral_outlet_mesh = Tessellation.from_stl(point_path + "/INTEGRAL_OUTLET.stl", airtight=False)
# normalize meshes
def normalize_mesh(mesh, center, scale):
mesh = mesh.translate([-c for c in center])
mesh = mesh.scale(scale)
return mesh
# scale and normalize mesh and openfoam data
center = (0.05213, 0.081092, 0.0)
scale = 1
nu = 0.025
inlet_vel = 0.5
inlet_branch_vel = inlet_vel
inlet_main_mesh = normalize_mesh(inlet_main_mesh, center, scale)
inlet_branch_mesh = normalize_mesh(inlet_branch_mesh, center, scale)
outlet_mesh = normalize_mesh(outlet_mesh, center, scale)
wall_pipe_mesh = normalize_mesh(wall_pipe_mesh, center, scale)
integral_main_mesh = normalize_mesh(integral_main_mesh, center, scale)
integral_branch_mesh = normalize_mesh(integral_branch_mesh, center, scale)
integral_outlet_mesh = normalize_mesh(integral_outlet_mesh, center, scale)
internal_volume_mesh = normalize_mesh(internal_volume_mesh, center, scale)
# geom params
inlet_main_normal = (1, 0, 0)
inlet_main_area = 0.0786428 * (scale**2)
inlet_radius = 0.316435 / 2
inlet_main_center = ((-0.429621 - center[0]) * scale, (0.009655 - center[1]) * scale, (0.00018 - center[2]) * scale)
inlet_branch_normal = (0, 0, -1)
inlet_branch_area = 0.0888933 * (scale**2)
inlet_branch_radius = 0.336426 / 2
inlet_branch_center = ((0.0 - center[0]) * scale, (0.488395 - center[1]) * scale, (0.0 - center[2]) * scale)
outlet_normal = (1, 0, 0)
outlet_area = 0.0786428 * (scale**2)
outlet_radius = 0.316435 / 2
outlet_center = ((0.570938 - center[0]) * scale, (0.009579 - center[1]) * scale, (0.001158 - center[2]) * scale)
# make domain
domain = Domain()
# 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")]
)
weight_integrals = 0.1
# add constraints to solver
# inlet_main
inlet_main = PointwiseBoundaryConstraint(
nodes=nodes,
geometry=inlet_main_mesh,
outvar={"u": inlet_vel, "v": 0, "w": 0},
batch_size=cfg.batch_size.inlet,
)
domain.add_constraint(inlet_main, "inlet_main")
# inlet_branch
inlet_branch = PointwiseBoundaryConstraint(
nodes=nodes,
geometry=inlet_branch_mesh,
outvar={"u": 0, "v": -inlet_vel, "w": 0},
batch_size=cfg.batch_size.inlet,
)
domain.add_constraint(inlet_branch, "inlet_branch")
# outlet
outlet = PointwiseBoundaryConstraint(
nodes=nodes,
geometry=outlet_mesh,
outvar={"p": 0},
batch_size=cfg.batch_size.outlet,
)
domain.add_constraint(outlet, "outlet")
# no slip
wall_pipe = PointwiseBoundaryConstraint(
nodes=nodes,
geometry=wall_pipe_mesh,
outvar={"u": 0, "v": 0, "w": 0},
batch_size=cfg.batch_size.no_slip,
)
domain.add_constraint(wall_pipe, "no_slip")
# interior
interior = PointwiseInteriorConstraint(
nodes=nodes,
geometry=internal_volume_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 -inlet main
integral_continuity = IntegralBoundaryConstraint(
nodes=nodes,
geometry=integral_main_mesh,
outvar={"normal_dot_vel": 0.0393214},
batch_size=1,
integral_batch_size=cfg.batch_size.integral_continuity,
lambda_weighting={"normal_dot_vel": weight_integrals},
)
domain.add_constraint(integral_continuity, "integral_main")
# Integral Continuity -integral branch
integral_continuity = IntegralBoundaryConstraint(
nodes=nodes,
geometry=integral_branch_mesh,
outvar={"normal_dot_vel": -0.04444665},
batch_size=1,
integral_batch_size=cfg.batch_size.integral_continuity,
lambda_weighting={"normal_dot_vel": weight_integrals},
)
domain.add_constraint(integral_continuity, "integral_branch") ##check for the name --integral_continuity_1
# Integral Continuity -integral outlet
integral_continuity = IntegralBoundaryConstraint(
nodes=nodes,
geometry=integral_outlet_mesh,
outvar={"normal_dot_vel": 0.08376805},
batch_size=1,
integral_batch_size=cfg.batch_size.integral_continuity,
lambda_weighting={"normal_dot_vel": weight_integrals},
)
domain.add_constraint(integral_continuity, "integral_outlet")
# add validation data
file_path = "./openfoam/openfoamOutput.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_main_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()
Feedback and suggestions are always welcome!
Thank you.
