Issues Emulating Simple Water Flow

[npstrike]

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

[ngeneva]

Hi @npstrike

I’m having a look at this now and will follow up early next week with any items I see. Thanks for the comprehensive details.

[npstrike]

Thanks @ngeneva, I appreciate it. I’ve also experimented with an additional integral plane along the tube, I’ll include the stl for that here in case it’s beneficial:

INTERNAL_INLET.stl (20.8 KB)

[ngeneva]

Hi @npstrike

So I did just some initial looking at the outputs and I am seeing some oddities with the points getting sampled from the STL file. I’m curious if you see the same. Namely with this section here:

I’m not sure whats going on exactly, but looking at the interior constraint I do see considerabley less points sampled in this area. Plotting the interior_volume vtp file in Paraview gives the following. I personally get a lot less points at the end of the pipe section.

I added a simple voxel inferencer and viewing this in Paraview with a threshold filter, seems the SDF isn’t right in the region either (the SDF is returning some incorrect value and these points are getting masked out as outside the pipe):

Inferencer code here (I edited the center to be center = (5.433142, 6.5, 0.000000) so its centered more around 0):

    from modulus.sym.domain.inferencer import VoxelInferencer

    mask_fn = (lambda x, y, z: interior.sdf({"x": x, "y": y, "z": z}, {})["sdf"] < 0)

    voxel_inferencer = VoxelInferencer(
        bounds = [[-0.04, 0.04], [-0.6, 0.6], [-0.02, 0.02]],
        npoints = [128, 512, 64],
        nodes=nodes,
        output_names=["u", "v", "w", "p"],
        export_map={"U": ["u", "v", "w"], "p": ["p"]},
        mask_fn = mask_fn,
        requires_grad=False,
        batch_size=1024,
    )
    domain.add_inferencer(voxel_inferencer, "vox_inf")

Edit: looking at that cut off point from the Voxel inferencer on the STL file I see some double layer. I’m guessing this is confusing pysdf any causing some sampling issues.

[npstrike]

@ngeneva yes, I’ve seen the same strange cutoff but I wasn’t sure what to make of it. I wasn’t aware of the VoxelInferencer feature, that’s a good tip! This odd behavior is likely because this is just a portion of a larger model, and I was trying to reduce the complexity for debugging by stripping it down to this tube. It would appear that I have had the opposite effect. I really appreciate you looking into this, I think I have some ideas on how I can move forwards now :-). If I have further questions can I send them your way?

[ngeneva]

Hi @npstrike

Great! PySdf does tend be to pretty picky with how the STL is formed. It can get confused on whats “inside” the mesh if things have some unenclosed parts. Feel free to send additional questions, I can continue to provide support (bandwidth permitting of course).

By the way, if this example (laminar flow through a curved pipe) starts working well and is able to be shared publicly, we would be very happy to have it as part of our example collection as well as potentially publicize it and in one our release blog posts. We can discuss more in the future if there’s interest.

[npstrike]

I have a question regarding creative use of the boundary constraint. One thing I’m playing with is using different types of helpers, similar to the use of the IntegralConstraint. However, what’s more intuitive for me and my use case (in-compressible flow) is to specify u,v, and w to a constraints outvar rather than the norm_dot_vel that’s proposed for IntegralConstraints .

As such I’ve been experimenting with using PointwiseBoundaryConstraints defined along the interior of the volume and using that to provide known velocity vectors. Even if this is odd, are there any reasons that this should fundamentally not work?

I’m still working getting my case to produce satisfactory output and I’d like to know if this is a red herring.

Thanks!

As for the example problem, I doubt my company would allow me to share the full sample of what I’m actually doing, however if I can get our sample working I can likely get a revised version of this simple pipe I provided above using corrected .stl’s if you guys are interested. I’m personally a big supporter of community development, but am limited by our IP policies. I’m happy to contribute what I can though if it helps the project.

[ngeneva]

As such I’ve been experimenting with using PointwiseBoundaryConstraint s defined along the interior of the volume and using that to provide known velocity vectors. Even if this is odd, are there any reasons that this should fundamentally not work?

If this is prior physical information that you would know ahead of time for the problem, then you should definitely try to incorporate it to speed up convergence. Similar to numerical solvers, there’s multiple approaches for boundary conditions that work better in some problems and worse in others. Experimenting with different constraint methods is exactly what we do when developing these examples.

[npstrike]

Ok, thanks! Yes, this is information I know prior, however I wasn’t sure if using a PointwiseBoundaryConstraint might confuse it somehow as I’ve only seen it used as a barrier between the scope of the run and the external world that’s not emulated.

Thanks for the feedback :)