NVIDIA / modulus-sym

Framework providing pythonic APIs, algorithms and utilities to be used with Modulus core to physics inform model training as well as higher level abstraction for domain experts
https://developer.nvidia.com/modulus
Apache License 2.0
165 stars 68 forks source link

šŸ›[BUG]: potential Sympy compilation error when using parameterized geometry in FPGA example #58

Open hardik1311 opened 1 year ago

hardik1311 commented 1 year ago

Version

modulus.sym.version '1.0.0'

On which installation method(s) does this occur?

No response

Describe the issue

I am trying to run fpga example with parameterized heatsink geometry.

sympy.version '1.5.1'

no_slip = PointwiseBoundaryConstraint(
        nodes=flow_nodes,
        geometry=geo.geo,
        outvar={"u": 0, "v": 0, "w": 0},
        batch_size=cfg.batch_size.no_slip,
        quasirandom=cfg.custom.quasirandom,
        parameterization=geo.pr,
    )
    flow_domain.add_constraint(no_slip, "no_slip")

leads to a runtime error

Minimum reproducible example

Made some changes to the fpga examples/ 

Code inside fpga_geometry.py

# Copyright (c) 2023, NVIDIA CORPORATION & AFFILIATES. All rights reserved.
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
#     http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.

from sympy import Symbol, Eq, tanh
import numpy as np
import yaml

from modulus.sym.geometry.primitives_3d import Box, Channel, Plane
from modulus.sym.geometry import Parameterization, Parameter

# # check if the geometry is parameterized
# with open("./conf/config.yaml", "r") as file:
#     yaml_data = yaml.safe_load(file)
#     parameterized = yaml_data["custom"]["parameterized"]

# if parameterized:
#     # geometry parameterization
#     HS_height = Symbol("HS_height")
#     HS_length = Symbol("HS_length")
#     Source_dim = Symbol("Source_dim")
# else:
#     HS_length = 0.65
#     HS_height = 0.8625
#     Source_dim = 0.25

# define sympy varaibles to parametize domain curves
x, y, z = Symbol("x"), Symbol("y"), Symbol("z")
x_pos = Symbol("x_pos")
HS_height = Symbol("HS_height")
HS_length = Symbol("HS_length")
# Source_dim = Symbol("Source_dim")
# define range for each of the parameters
HS_height_range = (0.40625, 0.8625)
HS_length_range = (0.35, 0.65)
# Source_dim_range = (0.2, 0.3)
# parameter range dict
param_ranges = {
    HS_height: HS_height_range,
    HS_length: HS_length_range,
    # Source_dim: Source_dim_range,
}
fixed_param_ranges = {
    HS_height: 0.65,
    HS_length: 0.8625,
    # Source_dim: 0.25,
}
# geometry params for domain
channel_origin = (-2.5, -0.5, -0.5625)
channel_dim = (5.0, 1.0, 1.125)
heat_sink_base_origin = (-0.75, -0.5, -0.4375)
heat_sink_base_dim = (HS_length, 0.05, 0.875)
fin_origin = heat_sink_base_origin
fin_dim = (HS_length, HS_height, 0.0075)
total_fins = 17
flow_box_origin = (-0.85, -0.5, -0.5625)
flow_box_dim = (0.85, 1.0, 1.125)
source_origin = (-0.55, -0.5, -0.125)
source_dim = (0.25, 0.0, 0.25)

class FPGAChannel(object):
    def __init__(self, parameterized: bool = False):
        if parameterized:
            pr = Parameterization(param_ranges)
            self.pr = param_ranges
        else:
            pr = Parameterization(fixed_param_ranges)
            self.pr = fixed_param_ranges        
        # define geometry
        # channel
        self.channel = Channel(
            channel_origin,
            (
                channel_origin[0] + channel_dim[0],
                channel_origin[1] + channel_dim[1],
                channel_origin[2] + channel_dim[2],
            ),
            parameterization=pr,
        )
        # fpga heat sink
        heat_sink_base = Box(
            heat_sink_base_origin,
            (
                heat_sink_base_origin[0] + heat_sink_base_dim[0],  # base of heat sink
                heat_sink_base_origin[1] + heat_sink_base_dim[1],
                heat_sink_base_origin[2] + heat_sink_base_dim[2],
            ),
            parameterization=pr,
        )
        fin_center = (
            fin_origin[0] + fin_dim[0] / 2,
            fin_origin[1] + fin_dim[1] / 2,
            fin_origin[2] + fin_dim[2] / 2,
        )
        fin = Box(
            fin_origin,
            (
                fin_origin[0] + fin_dim[0],
                fin_origin[1] + fin_dim[1],
                fin_origin[2] + fin_dim[2],
            ),
            parameterization=pr,
        )
        # gap between fins
        gap = (heat_sink_base_dim[2] - fin_dim[2]) / (total_fins - 1)
        fin = fin.repeat(
            gap,
            repeat_lower=(0, 0, 0),
            repeat_higher=(0, 0, total_fins - 1),
            center=fin_center,
        )
        self.fpga = heat_sink_base + fin
        # entire geometry
        self.geo = self.channel - self.fpga
        # inlet and outlet
        self.inlet = Plane(
            channel_origin,
            (
                channel_origin[0],
                channel_origin[1] + channel_dim[1],
                channel_origin[2] + channel_dim[2],
            ),
            -1,
            parameterization=pr,
        )
        self.outlet = Plane(
            (channel_origin[0] + channel_dim[0], channel_origin[1], channel_origin[2]),
            (
                channel_origin[0] + channel_dim[0],
                channel_origin[1] + channel_dim[1],
                channel_origin[2] + channel_dim[2],
            ),
            1,
            parameterization=pr,
        )
        # planes for integral continuity
        x_pos_range = {x_pos: (-0.75, 0.0)}
        self.integral_plane = Plane(
            (x_pos, channel_origin[1], channel_origin[2]),
            (x_pos, channel_origin[1] + channel_dim[1], channel_origin[2] + channel_dim[2]),
            1,
            parameterization=Parameterization(x_pos_range),
        )

Code inside fpga_flow.py

# Copyright (c) 2023, NVIDIA CORPORATION & AFFILIATES. All rights reserved.
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
#     http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.

import os
import warnings

import sys
import torch
from sympy import Symbol, Eq, Abs, tanh, And, Or
import numpy as np

import modulus.sym
from modulus.sym.hydra import to_absolute_path, instantiate_arch, ModulusConfig
from modulus.sym.utils.io import csv_to_dict
from modulus.sym.solver import Solver
from modulus.sym.domain import Domain
from modulus.sym.geometry.primitives_3d import Box, Channel, Plane
from modulus.sym.domain.constraint import (
    PointwiseBoundaryConstraint,
    PointwiseInteriorConstraint,
    IntegralBoundaryConstraint,
)
from modulus.sym.domain.validator import PointwiseValidator
from modulus.sym.domain.inferencer import PointwiseInferencer
from modulus.sym.domain.monitor import PointwiseMonitor
from modulus.sym.key import Key
from modulus.sym.node import Node
from modulus.sym.eq.pdes.navier_stokes import NavierStokes, Curl
from modulus.sym.eq.pdes.basic import NormalDotVec, GradNormal
from modulus.sym.models.fully_connected import FullyConnectedArch
from modulus.sym.models.fourier_net import FourierNetArch
from modulus.sym.models.siren import SirenArch
from modulus.sym.models.modified_fourier_net import ModifiedFourierNetArch
from modulus.sym.models.dgm import DGMArch

from fpga_geometry import *

@modulus.sym.main(config_path="conf", config_name="config")
def run(cfg: ModulusConfig) -> None:
    # params for simulation
    # fluid params
    nu = 0.02
    rho = 1
    inlet_vel = 1.0
    volumetric_flow = 1.125

    # make list of nodes to unroll graph on
    ns = NavierStokes(nu=nu, rho=rho, dim=3, time=False)
    normal_dot_vel = NormalDotVec()
    equation_nodes = ns.make_nodes() + normal_dot_vel.make_nodes()

    # determine inputs outputs of the network
    # determine inputs outputs of the network
    input_keys = [Key("x"), Key("y"), Key("z")]
    if cfg.custom.parameterized:
        input_keys += [Key("HS_height"), Key("HS_length")]
        monitor_fixed_param_ranges = {
            HS_height: lambda batch_size: np.full(
                (batch_size, 1), np.random.uniform(*HS_height_range)
            ),
            HS_length: lambda batch_size: np.full(
                (batch_size, 1), np.random.uniform(*HS_length_range)
            ),
        }
    else:
        monitor_fixed_param_ranges = {}
    # determine inputs outputs of the network
    output_keys = [Key("u"), Key("v"), Key("w"), Key("p")]
    # select the network and the specific configs
    if cfg.custom.arch == "FullyConnectedArch":
        flow_net = FullyConnectedArch(
            input_keys=input_keys,
            output_keys=output_keys,
            adaptive_activations=cfg.custom.adaptive_activations,
        )
    elif cfg.custom.arch == "FourierNetArch":
        flow_net = FourierNetArch(
            input_keys=input_keys,
            output_keys=output_keys,
            adaptive_activations=cfg.custom.adaptive_activations,
        )
    elif cfg.custom.arch == "SirenArch":
        flow_net = SirenArch(
            input_keys=input_keys,
            output_keys=output_keys,
            normalization={"x": (-2.5, 2.5), "y": (-2.5, 2.5), "z": (-2.5, 2.5)},
        )
    elif cfg.custom.arch == "ModifiedFourierNetArch":
        flow_net = ModifiedFourierNetArch(
            input_keys=input_keys,
            output_keys=output_keys,
            adaptive_activations=cfg.custom.adaptive_activations,
        )
    elif cfg.custom.arch == "DGMArch":
        flow_net = DGMArch(
            input_keys=input_keys,
            output_keys=output_keys,
            layer_size=128,
            adaptive_activations=cfg.custom.adaptive_activations,
        )
    else:
        sys.exit(
            "Network not configured for this script. Please include the network in the script"
        )

    flow_nodes = equation_nodes + [flow_net.make_node(name="flow_network")]

    # make flow domain
    flow_domain = Domain()
    geo = FPGAChannel(parameterized=cfg.custom.parameterized)

    # inlet
    def channel_sdf(x, y, z):
        sdf = geo.channel.sdf({"x": x, "y": y, "z": z}, {})
        return sdf["sdf"]

    constraint_inlet = PointwiseBoundaryConstraint(
        nodes=flow_nodes,
        geometry=geo.inlet,
        outvar={"u": inlet_vel, "v": 0, "w": 0},
        batch_size=cfg.batch_size.inlet,
        criteria=Eq(x, channel_origin[0]),
        lambda_weighting={"u": channel_sdf, "v": 1.0, "w": 1.0},  # weight zero on edges
        quasirandom=cfg.custom.quasirandom,
        parameterization=geo.pr,
    )
    flow_domain.add_constraint(constraint_inlet, "inlet")

    # outlet
    constraint_outlet = PointwiseBoundaryConstraint(
        nodes=flow_nodes,
        geometry=geo.outlet,
        outvar={"p": 0},
        batch_size=cfg.batch_size.outlet,
        criteria=Eq(x, channel_origin[0] + channel_dim[0]),
        quasirandom=cfg.custom.quasirandom,
        parameterization=geo.pr,
    )
    flow_domain.add_constraint(constraint_outlet, "outlet")

    # no slip for channel walls
    print(type(geo.pr))
    print(geo.pr)
    no_slip = PointwiseBoundaryConstraint(
        nodes=flow_nodes,
        geometry=geo.geo,
        outvar={"u": 0, "v": 0, "w": 0},
        batch_size=cfg.batch_size.no_slip,
        quasirandom=cfg.custom.quasirandom,
        parameterization=geo.pr,
    )
    flow_domain.add_constraint(no_slip, "no_slip")

    # flow interior low res away from fpga
    lr_interior = PointwiseInteriorConstraint(
        nodes=flow_nodes,
        geometry=geo.geo,
        outvar={"continuity": 0, "momentum_x": 0, "momentum_y": 0, "momentum_z": 0},
        batch_size=cfg.batch_size.lr_interior,
        criteria=Or(x < flow_box_origin[0], x > (flow_box_origin[0] + flow_box_dim[0])),
        lambda_weighting={
            "continuity": Symbol("sdf"),
            "momentum_x": Symbol("sdf"),
            "momentum_y": Symbol("sdf"),
            "momentum_z": Symbol("sdf"),
        },
        quasirandom=cfg.custom.quasirandom,
        parameterization=geo.pr,
    )
    flow_domain.add_constraint(lr_interior, "lr_interior")

    # flow interiror high res near fpga
    hr_interior = PointwiseInteriorConstraint(
        nodes=flow_nodes,
        geometry=geo.geo,
        outvar={"continuity": 0, "momentum_x": 0, "momentum_z": 0, "momentum_y": 0},
        batch_size=cfg.batch_size.hr_interior,
        criteria=And(
            x > flow_box_origin[0], x < (flow_box_origin[0] + flow_box_dim[0])
        ),
        lambda_weighting={
            "continuity": Symbol("sdf"),
            "momentum_x": Symbol("sdf"),
            "momentum_y": Symbol("sdf"),
            "momentum_z": Symbol("sdf"),
        },
        quasirandom=cfg.custom.quasirandom,
        parameterization=geo.pr,
    )
    flow_domain.add_constraint(hr_interior, "hr_interior")

    # integral continuity
    def integral_criteria(invar, params):
        sdf = geo.sdf(invar, params)
        return np.greater(sdf["sdf"], 0)

    integral_continuity = IntegralBoundaryConstraint(
        nodes=flow_nodes,
        geometry=geo.integral_plane,
        outvar={"normal_dot_vel": volumetric_flow},
        batch_size=cfg.batch_size.num_integral_continuity,
        integral_batch_size=cfg.batch_size.integral_continuity,
        criteria=integral_criteria,
        lambda_weighting={"normal_dot_vel": 1.0},        
        quasirandom=cfg.custom.quasirandom,
        parameterization={**geo.pr, **{x_pos: (-0.75, 0.0)}},
    )
    flow_domain.add_constraint(integral_continuity, "integral_continuity")

    # flow data
    file_path = "../openfoam/fpga_heat_fluid0.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_rgh": "p",
        }
        openfoam_var = csv_to_dict(to_absolute_path(file_path), mapping)
        openfoam_var["x"] = openfoam_var["x"] + channel_origin[0]
        openfoam_var["y"] = openfoam_var["y"] + channel_origin[1]
        openfoam_var["z"] = openfoam_var["z"] + channel_origin[2]
        if cfg.custom.parametrized:            
            openfoam_var["HS_height"] = (
                np.ones_like(openfoam_var["x"])
                * fixed_param_ranges[Symbol("HS_height")]
            )
            openfoam_var["HS_length"] = (
                np.ones_like(openfoam_var["x"])
                * fixed_param_ranges[Symbol("HS_length")]
            )
            openfoam_invar_numpy = {
                key: value
                for key, value in openfoam_var.items()
                if key in ["x", "y", "z", "HS_height", "HS_length"]
            }
        else:
            openfoam_invar_numpy = {
                key: value
                for key, value in openfoam_var.items()
                if key in ["x", "y", "z"]
            }
        openfoam_outvar_numpy = {
            key: value
            for key, value in openfoam_var.items()
            if key in ["u", "v", "w", "p"]
        }
        openfoam_validator = PointwiseValidator(
            nodes=flow_nodes,
            invar=openfoam_invar_numpy,
            true_outvar=openfoam_outvar_numpy,
            requires_grad=cfg.custom.exact_continuity,
        )
        flow_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
    invar_front_pressure = geo.integral_plane.sample_boundary(
        1024,
        parameterization={**fixed_param_ranges,**{x_pos: heat_sink_base_origin[0] - heat_sink_base_dim[0]}},
    )
    pressure_monitor = PointwiseMonitor(
        invar_front_pressure,
        output_names=["p"],
        metrics={"front_pressure": lambda var: torch.mean(var["p"])},
        nodes=flow_nodes,
    )
    flow_domain.add_monitor(pressure_monitor)
    invar_back_pressure = geo.integral_plane.sample_boundary(
        1024,
        parameterization={**fixed_param_ranges,**{x_pos: heat_sink_base_origin[0] + 2 * heat_sink_base_dim[0]}},
    )
    pressure_monitor = PointwiseMonitor(
        invar_back_pressure,
        output_names=["p"],
        metrics={"back_pressure": lambda var: torch.mean(var["p"])},
        nodes=flow_nodes,
    )
    flow_domain.add_monitor(pressure_monitor)

    # make solver
    flow_slv = Solver(cfg, flow_domain)

    # start flow solver
    flow_slv.solve()

if __name__ == "__main__":
    run()

and conf

# Copyright (c) 2023, NVIDIA CORPORATION & AFFILIATES. All rights reserved.
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
#     http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.

# This script allows you to try various combinations of Modulus
# features by changing only this config.yaml file. For best
# performance, we recommend the below defaults for each architecture.
# You can modify them by editing the  correponding entries in this 
# file
#
# Arch                          Start Lr    Max Steps   Decay Steps
# FullyConnectedArch            1.00E-03    1500000         15000      
# FourierNetArch                1.00E-03    400000          7500       
# ModifiedFourierNetArch    1.00E-03    400000          7500       
# SirenArch                     2.00E-05    500000          5000       
# DGMArch                       1.00E-03        1500000         15000           

# WARNING: Setting "exact_continuity" to true or setting the arch
# as "ModifiedFourierNetArch" increases the memory requirements of the 
# problem. Batchsizes may need to be reduced for such cases.  

defaults:
  - modulus_default
  - scheduler: tf_exponential_lr
  - optimizer: adam
  - loss: sum
  - _self_

scheduler: 
  decay_rate: 0.95
  decay_steps: 7500            # Change this based on arch chosen

optimizer: 
  lr: 1e-3                      # Change this based on arch chosen

training:
  rec_validation_freq: 1000
  rec_inference_freq: 2000
  rec_monitor_freq: 1000
  rec_constraint_freq: 2000
  max_steps: 400000            # Change this based on arch chosen

custom:
  arch: "ModifiedFourierNetArch"
  quasirandom: false
  exact_continuity: false
  adaptive_activations: false
  parameterized: true

network_dir: "network_checkpoint_flow"
initialization_network_dir: ""
save_filetypes: "vtk"

batch_size:
  inlet: 560
  outlet: 560
  no_slip: 20000
  lr_interior: 2500
  hr_interior: 2500
  integral_continuity: 11250
  num_integral_continuity: 5

graph:
  func_arch: true

Relevant log output

Error executing job with overrides: []
AttributeError: 'Add' object has no attribute 'rint'
The above exception was the direct cause of the following exception:
Traceback (most recent call last):
  File "/home/hardik/.local/lib/python3.8/site-packages/numpy/core/fromnumeric.py", line 57, in _wrapfunc
    return bound(*args, **kwds)
TypeError: loop of ufunc does not support argument 0 of type Add which has no callable rint method
During handling of the above exception, another exception occurred:
AttributeError: 'Add' object has no attribute 'rint'
The above exception was the direct cause of the following exception:
Traceback (most recent call last):
  File "fpga_flow.py", line 304, in <module>
    run()
  File "/home/hardik/.local/lib/python3.8/site-packages/modulus/sym/hydra/utils.py", line 104, in func_decorated
    _run_hydra(
  File "/home/hardik/.local/lib/python3.8/site-packages/hydra/_internal/utils.py", line 394, in _run_hydra
    _run_app(
  File "/home/hardik/.local/lib/python3.8/site-packages/hydra/_internal/utils.py", line 457, in _run_app
    run_and_report(
  File "/home/hardik/.local/lib/python3.8/site-packages/hydra/_internal/utils.py", line 223, in run_and_report
    raise ex
  File "/home/hardik/.local/lib/python3.8/site-packages/hydra/_internal/utils.py", line 220, in run_and_report
    return func()
  File "/home/hardik/.local/lib/python3.8/site-packages/hydra/_internal/utils.py", line 458, in <lambda>
    lambda: hydra.run(
  File "/home/hardik/.local/lib/python3.8/site-packages/hydra/_internal/hydra.py", line 132, in run
    _ = ret.return_value
  File "/home/hardik/.local/lib/python3.8/site-packages/hydra/core/utils.py", line 260, in return_value
    raise self._return_value
  File "/home/hardik/.local/lib/python3.8/site-packages/hydra/core/utils.py", line 186, in run_job
    ret.return_value = task_function(task_cfg)
  File "fpga_flow.py", line 154, in run
    no_slip = PointwiseBoundaryConstraint(
  File "/home/hardik/.local/lib/python3.8/site-packages/modulus/sym/domain/constraint/continuous.py", line 291, in __init__
    invar = geometry.sample_boundary(
  File "/home/hardik/.local/lib/python3.8/site-packages/modulus/sym/geometry/geometry.py", line 488, in sample_boundary
    [
  File "/home/hardik/.local/lib/python3.8/site-packages/modulus/sym/geometry/geometry.py", line 489, in <listcomp>
    curve.approx_area(parameterization, criteria=closed_boundary_criteria)
  File "/home/hardik/.local/lib/python3.8/site-packages/modulus/sym/geometry/curve.py", line 148, in approx_area
    computed_criteria = criteria(s, p)
  File "/home/hardik/.local/lib/python3.8/site-packages/modulus/sym/geometry/geometry.py", line 480, in boundary_criteria
    return self.boundary_criteria(invar, criteria=criteria, params=params)
  File "/home/hardik/.local/lib/python3.8/site-packages/modulus/sym/geometry/geometry.py", line 401, in boundary_criteria
    sdf_normal_plus = self.sdf(
  File "/home/hardik/.local/lib/python3.8/site-packages/modulus/sym/geometry/geometry.py", line 701, in sub_sdf
    computed_sdf_2 = sdf_2(invar, params, compute_sdf_derivatives)
  File "/home/hardik/.local/lib/python3.8/site-packages/modulus/sym/geometry/geometry.py", line 668, in add_sdf
    computed_sdf_2 = sdf_2(invar, params, compute_sdf_derivatives)
  File "/home/hardik/.local/lib/python3.8/site-packages/modulus/sym/geometry/geometry.py", line 348, in repeat_sdf
    np.maximum(np.around(clamped_invar[d] / spacing), rl), rh
  File "<__array_function__ internals>", line 200, in around
  File "/home/hardik/.local/lib/python3.8/site-packages/numpy/core/fromnumeric.py", line 3337, in around
    return _wrapfunc(a, 'round', decimals=decimals, out=out)
  File "/home/hardik/.local/lib/python3.8/site-packages/numpy/core/fromnumeric.py", line 66, in _wrapfunc
    return _wrapit(obj, method, *args, **kwds)
  File "/home/hardik/.local/lib/python3.8/site-packages/numpy/core/fromnumeric.py", line 43, in _wrapit
    result = getattr(asarray(obj), method)(*args, **kwds)
TypeError: loop of ufunc does not support argument 0 of type Add which has no callable rint method

Environment details

No response

Other/Misc.

No response