Hello, I’m attempting to use modulus to emulate a 2 gpm water flow in a ~1/2in pipe with two 90 degree bends. The water flows into the inlet
, and out of the inlet_exit
. Unfortunately, I am having significant difficulties getting it to converge onto the correct solution. If I render it’s output in ParaView while training I can see it occasionally produces a result close to what I want but it will not stay there. I would appreciate any guidance that can be offered, as I am attempting to evaluate whether this tool will be useful for my company. I have an embarrassing 2 weeks invested into this, but there’s a good chance it’s a relatively basic error that I’m making but cannot see.
This project is based off of the anuerism sample.
Thanks!
Project files
my_flow.zip (668.3 KB)
Run Environment
Host OS | Ubuntu 22.04 |
Modulus Version | 22.09 |
Modulus Type | Docker Image |
GPU | RTX 3090 |
CPU | AMD Ryzen 3950x |
RAM | 64GB |
Run file (from my_flow.zip
):
from sympy import Symbol, sqrt, Max
[my_flow.zip|attachment](upload://eQhFUZNLmYul9FqWq7gBQ3nPCH6.zip) (668.3 KB)
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.key import Key
from modulus.eq.pdes.navier_stokes import NavierStokes
from modulus.eq.pdes.basic import NormalDotVec
from modulus.geometry.tessellation import Tessellation
@modulus.main(config_path="conf", config_name="config")
def run(cfg: ModulusConfig) -> None:
# read stl files to make geometry
point_path = to_absolute_path("./STL")
mm_to_in = 0.0393701
# stl's are provided in inches (except INLET_EXIT.stl which is in mm)
wall_inlet_tube_mesh = Tessellation.from_stl(point_path + "/WALL_INLET_TUBE.stl", airtight=False)
interior = Tessellation.from_stl(point_path + "/WALL_INLET_TUBE.stl", airtight=True)
inlet_mesh = Tessellation.from_stl(point_path + "/INLET.stl", airtight=False)
inlet_exit_mesh = Tessellation.from_stl(point_path + "/INLET_EXIT.stl", airtight=False)
inlet_exit_mesh = inlet_exit_mesh.scale(mm_to_in) # Convert model from mm to inches
# Use inlet_exit mesh to create planes throughout the tube at various heights
helper_1_offset = 50
helper_2_offset = 25
inlet_helper_1_mesh = inlet_exit_mesh.translate([c for c in (0,helper_1_offset,0)])
inlet_helper_2_mesh = inlet_exit_mesh.translate([c for c in (0,helper_2_offset,0)])
# 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
# scale and normalize mesh and openfoam data
center = (5.433142, 30.313574, 0.000000)
scale = (1 / 54) # stl is 60in tall, 1524mm = 60in
inlet_mesh = normalize_mesh(inlet_mesh, center, scale)
wall_inlet_tube_mesh = normalize_mesh(wall_inlet_tube_mesh, center, scale)
interior = normalize_mesh(interior, center, scale)
inlet_exit_mesh = normalize_mesh(inlet_exit_mesh, center, scale)
inlet_helper_1_mesh = normalize_mesh(inlet_helper_1_mesh, center, scale)
inlet_helper_2_mesh = normalize_mesh(inlet_helper_2_mesh, center, scale)
# params
inlet_normal = (0.0, 1.0, 0.0)
inlet_center = ((7.0 - center[0]) * scale, (29.0 - center[1]) * scale, (0.0 - center[2]) * scale)
inlet_radius = 0.585 / 2 * scale
inlet_vel = 2.39 * 12 * scale # feet per sec -> in per sec
outlet_normal = (0.0, -1.0, 0.0)
outlet_center = (((101.599998 * mm_to_in) - center[0]) * scale, ((-553.211975 * mm_to_in) - center[1]) * scale, ((0.0 * mm_to_in) - center[2]) * scale)
outlet_radius = 0.585 / 2 * scale
outlet_vel = inlet_vel # feet per sec -> in per sec
domain = Domain()
# make list of nodes to unroll graph on
nu = 0.564 # Kinematic viscosity calculated here: https://www.omnicalculator.com/physics/water-viscosity
rho = 0.9885e-6 # Density Calculated here https://www.omnicalculator.com/physics/water-viscosity
ns = NavierStokes(nu=nu * scale, rho=rho, 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")]
)
# add constraints to solver
u_inlet, v_inlet, w_inlet = circular_parabola(
Symbol("x"),
Symbol("y"),
Symbol("z"),
center=inlet_center,
normal=inlet_normal,
radius=inlet_radius,
max_vel=inlet_vel,
)
u_outlet, v_outlet, w_outlet = circular_parabola(
Symbol("x"),
Symbol("y"),
Symbol("z"),
center=outlet_center,
normal=outlet_normal,
radius=outlet_radius,
max_vel=outlet_vel,
)
inlet = PointwiseBoundaryConstraint(
nodes=nodes,
geometry=inlet_mesh,
outvar={"u": u_inlet, "v": v_inlet, "w": w_inlet},
batch_size=cfg.batch_size.inlet
)
domain.add_constraint(inlet, "inlet")
inlet_exit = PointwiseBoundaryConstraint(
nodes=nodes,
geometry=inlet_exit_mesh,
# Testing various outvar's
# outvar={"u": u_outlet, "v": v_outlet, "w": w_outlet},
# outvar={"p": 0},
outvar={"u": u_outlet, "v": v_outlet, "w": w_outlet, "p": 0},
batch_size=cfg.batch_size.inlet_exit,
)
domain.add_constraint(inlet_exit, "inlet_exit")
inlet_tube = PointwiseBoundaryConstraint(
nodes=nodes,
geometry=wall_inlet_tube_mesh,
outvar={"u": 0, "v": 0, "w": 0},
batch_size=cfg.batch_size.inlet_tube,
)
domain.add_constraint(inlet_tube, "inlet_tube")
internal_volume = PointwiseInteriorConstraint(
nodes=nodes,
geometry=interior,
outvar={"continuity": 0, "momentum_x": 0, "momentum_y": 0, "momentum_z": 0},
batch_size=cfg.batch_size.internal_volume,
)
domain.add_constraint(internal_volume, "internal_volume")
# Uncommenting this will give my best attempt at integral continuity, but I'm not convinced it was right
# volumetric flow calculated using 2 gpm: https://www.xconvert.com/unit-converter/gallons-per-minute-to-cubic-inches-per-second
# in3_per_sec = 7.7
#
# integral_lambda_weight = 1e-6 # Picked experimentally to match the scale of aggregate loss
# integral_continuity = IntegralBoundaryConstraint(
# nodes=nodes,
# geometry=inlet_helper_1_mesh,
# outvar={"normal_dot_vel": -in3_per_sec},
# batch_size=1,
# integral_batch_size=cfg.batch_size.integral_continuity,
# lambda_weighting={"normal_dot_vel": integral_lambda_weight},
# )
# domain.add_constraint(integral_continuity, "integral_continuity")
#
# integral_continuity = IntegralBoundaryConstraint(
# nodes=nodes,
# geometry=inlet_helper_2_mesh,
# outvar={"normal_dot_vel": -in3_per_sec},
# batch_size=1,
# integral_batch_size=cfg.batch_size.integral_continuity,
# lambda_weighting={"normal_dot_vel": integral_lambda_weight},
# )
# domain.add_constraint(integral_continuity, "integral_continuity")
# make solver
slv = Solver(cfg, domain)
# start solver
slv.solve()
if __name__ == "__main__":
run()
Config (from my_flow.zip
)
defaults :
- modulus_default
- arch:
- fully_connected
- scheduler: tf_exponential_lr
- optimizer: adam
- loss: sum
- _self_
run_mode: 'train'
scheduler:
decay_rate: 0.95
decay_steps: 150000
training:
rec_results_freq : 1000
rec_constraint_freq: 1000
max_steps : 150000
batch_size:
inlet: 100
inlet_helper: 100
inlet_exit: 100
inlet_tube: 5000
internal_volume: 600
integral_continuity: 50