Could not unroll graph!

Physics Simulation Modulus Physics-ML Model Framework Technical Support (Modulus Only)
1

I created my PDE and nodes like:

# Define the PDE
class wavePDE (PDE):
name = "AcousticWaveEquation "
def __init__(self, u="u", c="c", S="S", dim=3, time=True, mixed_form=False):
    # set params
    self.u = u
    self.dim = dim
    self.time = time
    self.mixed_form = mixed_form

    # coordinates
    x, y, z = Symbol("x"), Symbol("y"), Symbol("z")

    # time
    t = Symbol("t")

    # make input variables
    input_variables = {"x": x, "y": y, "z": z, "t": t}
    if self.dim == 1:
        input_variables.pop("y")
        input_variables.pop("z")
    elif self.dim == 2:
        input_variables.pop("z")
    if not self.time:
        input_variables.pop("t")

    # Scalar function (make u function)
    assert type(u) == str, "u needs to be string"
    u = Function(u)(*input_variables)

    # wave speed coefficient
    if type(c) is str:
        c = Function(c)(*input_variables)
    elif type(c) in [float, int]:
        c = Number(c)
        
    # Source Term
    if type(S) is str:
        S = Function(S)(*input_variables)
    elif type(S) in [float, int]:
        S = Number(S)

    # set equations
    self.equations = {}

    if not self.mixed_form: 
        self.equations["wave_PDE"] = (
            u.diff(t, 2)
            - c**2 * u.diff(x, 2)
            - c**2 * u.diff(y, 2)
            - c**2 * u.diff(z, 2)
            - S
        )
    elif self.mixed_form:
        u_x = Function("u_x")(*input_variables)
        u_y = Function("u_y")(*input_variables)
        if self.dim == 3:
            u_z = Function("u_z")(*input_variables)
        else:
            u_z = Number(0)
        if self.time:
            u_t = Function("u_t")(*input_variables)
        else:
            u_t = Number(0)

        self.equations["wave_PDE"] = (
            u_t.diff(t)
            - c**2 * u_x.diff(x)
            - c**2 * u_y.diff(y)
            - c**2 * u_z.diff(z)
        )
        self.equations["compatibility_u_x"] = u.diff(x) - u_x
        self.equations["compatibility_u_y"] = u.diff(y) - u_y
        self.equations["compatibility_u_z"] = u.diff(z) - u_z
        self.equations["compatibility_u_xy"] = u_x.diff(y) - u_y.diff(x)
        self.equations["compatibility_u_xz"] = u_x.diff(z) - u_z.diff(x)
        self.equations["compatibility_u_yz"] = u_y.diff(z) - u_z.diff(y)
        if self.dim == 2:
            self.equations.pop("compatibility_u_z")
            self.equations.pop("compatibility_u_xz")
            self.equations.pop("compatibility_u_yz")




wave_eq = wavePDE(u="u", c="c", S="S", dim=2, time=True)

# define networks
wave_net = instantiate_arch(
    input_keys=[Key("x"), Key("y"), Key("t")],
    output_keys=[Key("u")],
    cfg=cfg.arch.fully_connected,
)


#Finally, Create the nodes
nodes = wave_eq.make_nodes() + [wave_net.make_node("wave_network"),]



# add interior constraint
interior = PointwiseInteriorConstraint(
      nodes=nodes,
      geometry=geo,
      outvar={"wave_PDE": 0 },
      batch_size=4096,
      bounds={x: (0, L), y: (0, L)},
      #lambda_weighting={"wave_PDE": 0.0001},
      parameterization=time_range,
)
DefinedDomain.add_constraint(interior, "Interior")

# add open boundary constraint
bc = PointwiseBoundaryConstraint(
     nodes=nodes,
     geometry=geo,
     outvar={"u": 0},
     batch_size=1024,
     #lambda_weighting={"u": 0.01 * time_length},
     parameterization=time_range,
)
DefinedDomain.add_constraint(bc, "BC")

# initial condition
IC = PointwiseInteriorConstraint(
     nodes=nodes,
     geometry=geo,
     outvar={"u": 0},
     batch_size=1024,
     #lambda_weighting={"u": 1.0},
     parameterization={t: 0.0},
)
DefinedDomain.add_constraint(IC, "IC")

But adding my constraints fail with error:

####################################
could not unroll graph!
This is probably because you are asking to compute a value that is not an output of any node
####################################
invar: [x, y, sdf, area, t]
requested var: [wave_PDE]
computable var: [x, y, sdf, area, t, u]
####################################
Nodes in graph: 
node: Sympy Node: wave_PDE
evaluate: SympyToTorch
inputs: [S, c]
derivatives: [u__t__t, u__x__x, u__y__y]
outputs: [wave_PDE]
optimize: False
node: Arch Node: wave_network
evaluate: FullyConnectedArch
inputs: [x, y, t]
derivatives: []
outputs: [u]
optimize: True
####################################
---------------------------------------------------------------------------
RuntimeError                              Traceback (most recent call last)
Cell In [27], line 4
      1 # Specifying Constraints (bcs/ics/interior) sampling
      2 
      3 # add interior constraint
----> 4 interior = PointwiseInteriorConstraint(
      5     nodes=nodes,
      6     geometry=geo,
      7     outvar={"wave_PDE": 0 },
      8     batch_size=4096,
      9     bounds={x: (0, L), y: (0, L)},
     10     #lambda_weighting={"wave_PDE": 0.0001},
     11     parameterization=time_range,
     12 )
     13 DefinedDomain.add_constraint(interior, "Interior")
     15 # add open boundary constraint

File ~/.conda/envs/modulus_env/lib/python3.8/site-packages/modulus-22.9-py3.8.egg/modulus/domain/constraint/continuous.py:493, in PointwiseInteriorConstraint.__init__(self, nodes, geometry, outvar, batch_size, bounds, criteria, lambda_weighting, parameterization, fixed_dataset, compute_sdf_derivatives, importance_measure, batch_per_epoch, quasirandom, num_workers, loss, shuffle)
    486     dataset = ContinuousPointwiseIterableDataset(
    487         invar_fn=invar_fn,
    488         outvar_fn=outvar_fn,
    489         lambda_weighting_fn=lambda_weighting_fn,
    490     )
    492 # initialize constraint
--> 493 super().__init__(
    494     nodes=nodes,
    495     dataset=dataset,
    496     loss=loss,
    497     batch_size=batch_size,
    498     shuffle=shuffle,
    499     drop_last=True,
    500     num_workers=num_workers,
    501 )

File ~/.conda/envs/modulus_env/lib/python3.8/site-packages/modulus-22.9-py3.8.egg/modulus/domain/constraint/constraint.py:55, in Constraint.__init__(self, nodes, dataset, loss, batch_size, shuffle, drop_last, num_workers)
     44 self.dataloader = iter(
     45     Constraint.get_dataloader(
     46         dataset=self.dataset,
   (...)
     51     )
     52 )
     54 # construct model from nodes
---> 55 self.model = Graph(
     56     nodes,
     57     Key.convert_list(self.dataset.invar_keys),
     58     Key.convert_list(self.dataset.outvar_keys),
     59 )
     60 self.model.to(self.device)
     61 if self.manager.distributed:
     62     # https://pytorch.org/docs/master/notes/cuda.html#id5

File ~/.conda/envs/modulus_env/lib/python3.8/site-packages/modulus-22.9-py3.8.egg/modulus/graph.py:94, in Graph.__init__(self, nodes, invar, req_names, diff_nodes, func_arch, func_arch_allow_partial_hessian)
     92 if not set(req_names_no_diff).issubset(self.computable_names):
     93     _print_graph_unroll_error(nodes, invar, req_names)
---> 94     raise RuntimeError("Failed Unrolling Graph")
     96 # compute only necessary nodes for req_names
     97 # Walk backwards from the output nodes in the graph and keep adding required inputs
     98 # until all inputs are available in invar
     99 nodes = copy(nodes)

RuntimeError: Failed Unrolling Graph

Any suggestions on how to fix this?

2

Hi @cpe.sk

Thanks for your interest in Modulus.
Let me see if I can help you better understand this error:

####################################
could not unroll graph!
This is probably because you are asking to compute a value that is not an output of any node
####################################
invar: [x, y, sdf, area, t]
requested var: [wave_PDE]
computable var: [x, y, sdf, area, t, u]
####################################

The symbolic graph failed to unroll (for PointwiseInteriorConstraint as we see right under this print out); this tells you high level info here. What the input variables are, the ones that need to be computed, and the intermediate ones that can currently be computed…

So we know that we want wave_PDE but we only have the quantities x, y, sdf, area, t, u (u is from your NN which is revealed below:

Nodes in graph:
node: Sympy Node: wave_PDE
evaluate: SympyToTorch
inputs: [S, c]
derivatives: [u__t__t, u__x__x, u__y__y]
outputs: [wave_PDE]
optimize: False
node: Arch Node: wave_network
evaluate: FullyConnectedArch
inputs: [x, y, t]
derivatives:
outputs: [u]
optimize: True
####################################

Here we see 2 nodes: wave_PDE and wave_network. The wave PDE we want to execute to get output residual but can’t, so we take a look at the inputs / derivatives:
inputs: [S, c]
derivatives: [u__t__t, u__x__x, u__y__y]

Okay, so we can likely compute those derivatives because we have u from wave_network with inputs t,x,y. The problem here is that S,c don’t exist, so the computation isn’t possible.

Now that we know what the issue is, we can work on a solution. In your PointwiseInteriorConstraint the geometry is providing x and y. And the parameterization=time_range, is providing the t. So what you should do is add the missing variables S, c to your parameterization dictionary. E.g.
parameterization={Symbol("t"): (0.0,1.0), Symbol("S"): 1.0, Symbol("c"):1.0}

For more info, a good example of this is in the 3 fin example where complex parameterization are defined in the geometry file then used in the constraint.

1 Like
3

Thank you very much for the professional and detailed explanation. I sincerely appreciate it. It helped me understand a lot.
I have a follow-up question in this case:
“S” is the Source function (five points producing a signal from the specified locations) in my case which is supposed to be:

 # Add the source term constraint
 locx1, locx2, locx3, locx4, locx5 = 1.2, 0.7, 0.5, 0.7, 1.2
 locy1, locy2, locy3, locy4, locy5 = 4.3, 3.4, 2.5, 1.6, 0.7
 Amp = 2                                                                         #Amplitude of the source
 sigma = 0.2                                                                     #Width of the Gaussian impulse
 freq = 1                                                                        #Frequency of the source signal

 S_invar = {}
 S_invar["x"] = np.expand_dims(mesh_x.flatten(), axis=-1)
 S_invar["y"] = np.expand_dims(mesh_y.flatten(), axis=-1)
 S_invar["t"] = np.expand_dims(mesh_t.flatten(), axis=-1)

 GP = np.exp(-0.5*( ((S_invar["x"]-locx1)/sigma)**2 + ((S_invar["y"]-locy1)/sigma)**2 ))\
      +np.exp(-0.5*( ((S_invar["x"]-locx2)/sigma)**2 + ((S_invar["y"]-locy2)/sigma)**2 ))\
          +np.exp(-0.5*( ((S_invar["x"]-locx3)/sigma)**2 + ((S_invar["y"]-locy3)/sigma)**2 ))\
              +np.exp(-0.5*( ((S_invar["x"]-locx4)/sigma)**2 + ((S_invar["y"]-locy4)/sigma)**2 ))\
                  +np.exp(-0.5*( ((S_invar["x"]-locx5)/sigma)**2 + ((S_invar["y"]-locy5)/sigma)**2 ))        #Gaussian impulse formula
#S = GP*Amp*torch.sin(2*np.pi*freq*S_invar["t"])

S_outvar = {}
S_outvar["S"] = GP*Amp*np.sin(2*np.pi*freq*S_invar["t"])

and “c” is also a function like:

mesh_x, mesh_y, mesh_t = np.meshgrid(
    np.linspace(0, 5, 100), np.linspace(0, 5, 100), np.linspace(0, 5, 100), indexing="ij"
)
wave_speed_invar = {}
wave_speed_invar["x"] = np.expand_dims(mesh_x.flatten(), axis=-1)
wave_speed_invar["y"] = np.expand_dims(mesh_y.flatten(), axis=-1)
wave_speed_outvar = {}
wave_speed_outvar["c"] = np.tanh(80 * (wave_speed_invar["y"] - 0.5)) / 2 + 1.5

How can I correctly implement them? Please advise. Thanks again!

4

If your S and C are functions that depend on x,y,t then you need to provide Modulus a node to connect between x,y,t and s,c.

One option is to use the Node.from_sympy function. Theres a few example of this in the provided examples such as the turbulent channel problem.

1 Like
5

Thanks. What if I wanted to solve the 2D wave equation with a point source placed at coordinates (x,y), what is the best way to define the source point?

The case I’m trying to solve is a 2D wave equation with a source point function that looks something like this:
f(x=xs, y=ys) = sin (4 pi f t)

The difficulties I’m facing:

  1. Creating the source point(s) somewhere within the defined domain (for example at x= 0.5, y= 0.2), not on the boundary.
  2. Testing different frequencies f in the source function via parameterization while maintaining time dependency. Without the different frequency values, I use parameterization for t just like the wave_2d example in seismic waves, so is it still possible to add another parameterization for frequency f in addition to time t?

At this point I tried solving a 1D wave equation with a time dependent source point placed at one of the boundaries and it was successful. However, with the difficulties I mentioned, I do not know how to proceed to parameterizing the frequency f or expanding to 2D with source points within the domain.

validator_u
validator_u1481Ă—412 45.5 KB

I hope I explained my problem well. Please assist. Also, please ask me for further explanation if I missed something.

Update: I have looked at the example in Interface Problem by Variational Method — Modulus 22.09 documentation about using a Dirac function for creating a point source, but I don’t fully understand how to apply it to the wave equation class.

6

Hi @cpe.sk

  1. Creating the source point(s) somewhere within the defined domain (for example at x= 0.5, y= 0.2), not on the boundary.

If you want to add a source to the 2D wave equation / change the PDE, you’ll need to create a custom PDE class for your problem. This should be simple enough for the wave equations because you can copy the existing class and add your forcing equation into the definition of that equation. The DiracDelta function can work for a point source, although this may require having a constraint specifically to train on the source point.

The 1D wave equation example in our user guide, you’ve been working with creates a custom PDE.

so is it still possible to add another parameterization for frequency f in addition to time t ?

Yep, should be possible from my understanding here. Granted training will most likely take longer to converge and you may need a larger model. Just add it as another input to your network and define a range to sample f from in your parameterization

7

Thank you for your reply!
I haven’t tried parameterization for this example yet, but I’m trying to implement the source points and I got the the following error:

Traceback (most recent call last):
      File "wave_2d_source.py", line 171, in run
        wave_source_outvar["s_func"] = Amp*sin(2*np.pi*freq*wave_source_invar["t"])*(\
      File "C:\ProgramData\Anaconda3\envs\modulus_env\lib\site-packages\sympy\core\cache.py", line 96, in wrapper
        retval = func(*args, **kwargs)
      File "C:\ProgramData\Anaconda3\envs\modulus_env\lib\site-packages\sympy\core\function.py", line 473, in __new__
        result = super(Function, cls).__new__(cls, *args, **options)
      File "C:\ProgramData\Anaconda3\envs\modulus_env\lib\site-packages\sympy\core\cache.py", line 96, in wrapper
        retval = func(*args, **kwargs)
      File "C:\ProgramData\Anaconda3\envs\modulus_env\lib\site-packages\sympy\core\function.py", line 288, in __new__
        evaluated = cls.eval(*args)
      File "C:\ProgramData\Anaconda3\envs\modulus_env\lib\site-packages\sympy\functions\elementary\trigonometric.py", line 299, in eval
        if arg.could_extract_minus_sign():
    AttributeError: 'ImmutableDenseNDimArray' object has no attribute 'could_extract_minus_sign'

My original attempt code is :

import numpy as np
from sympy import Symbol, sin, exp, Eq, Number, Function

import modulus
from modulus.hydra import to_absolute_path, ModulusConfig, instantiate_arch
from modulus.solver import Solver
from modulus.domain import Domain
#from modulus.geometry.primitives_1d import Line1D
from modulus.geometry.primitives_2d import Rectangle
from modulus.domain.constraint import (
    PointwiseBoundaryConstraint,
    PointwiseInteriorConstraint,
    PointwiseConstraint,
)

from modulus.domain.validator import PointwiseValidator
from modulus.domain.inferencer import PointwiseInferencer
from modulus.key import Key
from modulus.node import Node
#from wave_equation import WaveEquation1D
#from modulus.eq.pdes.wave_equation import WaveEquation
from modulus.utils.io import (
    ValidatorPlotter,
    InferencerPlotter,
)
from modulus.eq.pde import PDE



class Wave2DEquation(PDE):
    """
    Wave equation

    Parameters
    ==========
    u : str
        The dependent variable.
    c : float, Sympy Symbol/Expr, str
        Wave speed coefficient. If `c` is a str then it is
        converted to Sympy Function of form 'c(x,y,z,t)'.
        If 'c' is a Sympy Symbol or Expression then this
        is substituted into the equation.
    dim : int
        Dimension of the wave equation (1, 2, or 3). Default is 2.
    time : bool
        If time-dependent equations or not. Default is True.
    mixed_form: bool
        If True, use the mixed formulation of the wave equation.

    Examples
    ========
    >>> we = WaveEquation(c=0.8, dim=3)
    >>> we.pprint()
      wave_equation: u__t__t - 0.64*u__x__x - 0.64*u__y__y - 0.64*u__z__z
    >>> we = WaveEquation(c='c', dim=2, time=False)
    >>> we.pprint()
      wave_equation: -c**2*u__x__x - c**2*u__y__y - 2*c*c__x*u__x - 2*c*c__y*u__y
    """

    name = "WaveEquation"

    def __init__(self, u="u", c="c", s_func="s_func", dim=3, time=True, mixed_form=False):
        # set params
        self.u = u
        self.dim = dim
        self.time = time
        self.mixed_form = mixed_form

        # coordinates
        x, y, z = Symbol("x"), Symbol("y"), Symbol("z")

        # time
        t = Symbol("t")

        # make input variables
        input_variables = {"x": x, "y": y, "z": z, "t": t}
        if self.dim == 1:
            input_variables.pop("y")
            input_variables.pop("z")
        elif self.dim == 2:
            input_variables.pop("z")
        if not self.time:
            input_variables.pop("t")

        # Scalar function
        assert type(u) == str, "u needs to be string"
        u = Function(u)(*input_variables)

        # wave speed coefficient
        if type(c) is str:
            c = Function(c)(*input_variables)
        elif type(c) in [float, int]:
            c = Number(c)
        
        # wave source function
        if type(s_func) is str:
            s_func = Function(s_func)(*input_variables)
        elif type(s_func) in [float, int]:
            s_func = Number(s_func)

        # set equations
        self.equations = {}

        if not self.mixed_form:
            self.equations["wave_equation"] = (
                u.diff(t, 2)
                - c**2 * u.diff(x, 2)
                - c**2 * u.diff(y, 2)
                - c**2 * u.diff(z, 2)
                - s_func
            )
        


@modulus.main(config_path="conf", config_name="config_source")
def run(cfg: ModulusConfig) -> None:

    # make list of nodes to unroll graph on
    we = Wave2DEquation(u="u", c=1.0, s_func="s_func", dim=2, time=True)
    wave_net = instantiate_arch(
        input_keys=[Key("x"), Key("y"), Key("t")],
        output_keys=[Key("u")],
        cfg=cfg.arch.fully_connected,
    )
    source_net = instantiate_arch(
        input_keys=[Key("x"), Key("y"), Key("t")],
        output_keys=[Key("s_func")],
        cfg=cfg.arch.fully_connected,
    )
   

    #creating nodes
    nodes = (
        we.make_nodes(detach_names=["s_func"]) 
        + [
            wave_net.make_node(name="wave_network"),
            source_net.make_node(name="source_network"),
        ]
    )
    


    # add constraints to solver
    # make geometry
    x, y, t_symbol = Symbol("x"), Symbol("y"), Symbol("t")
    L = 1.0
    geo = Rectangle((0,0),(L,L))
    time_range = {t_symbol: (0, L)}

    # make domain
    domain = Domain()

    # Add source points
    
    #Source function parameters
    locx1, locx2, locx3, locx4, locx5 = 1.2, 0.7, 0.5, 0.7, 1.2
    locy1, locy2, locy3, locy4, locy5 = 4.3, 3.4, 2.5, 1.6, 0.7
    Amp = 2                                                                         #Amplitude of the source
    sigma = 0.2                                                                     #Width of the Gaussian impulse
    freq = 1 

    mesh_x, mesh_y, mesh_t = np.meshgrid(
        np.linspace(0, 5, 128), np.linspace(0, 5, 128), np.linspace(0, 5, 128), indexing="ij"
    )
    wave_source_invar = {}
    wave_source_invar["x"] = np.expand_dims(mesh_x.flatten(), axis=-1)
    wave_source_invar["y"] = np.expand_dims(mesh_y.flatten(), axis=-1)
    wave_source_invar["t"] = np.expand_dims(mesh_t.flatten(), axis=-1)
    
    wave_source_outvar = {}
    wave_source_outvar["s_func"] = Amp*sin(2*np.pi*freq*wave_source_invar["t"])*(\
        exp(-0.5*( ((wave_source_invar["x"]-locx1)/sigma)**2 + ((wave_source_invar["y"]-locy1)/sigma)**2 ))\
                +exp(-0.5*( ((wave_source_invar["x"]-locx2)/sigma)**2 + ((wave_source_invar["y"]-locy2)/sigma)**2 ))\
                    +exp(-0.5*( ((wave_source_invar["x"]-locx3)/sigma)**2 + ((wave_source_invar["y"]-locy3)/sigma)**2 ))\
                        +exp(-0.5*( ((wave_source_invar["x"]-locx4)/sigma)**2 + ((wave_source_invar["y"]-locy4)/sigma)**2 ))\
                            +exp(-0.5*( ((wave_source_invar["x"]-locx5)/sigma)**2 + ((wave_source_invar["y"]-locy5)/sigma)**2 ))\
    )
                


    # add source points constraint
    source = PointwiseConstraint.from_numpy(
        nodes=nodes, invar=wave_source_invar, outvar=wave_source_outvar, batch_size=1024
    )
    domain.add_constraint(source, "Source")

    # initial condition
    IC = PointwiseInteriorConstraint(
        nodes=nodes,
        geometry=geo,
        outvar={"u": 0, "u__t": 0},
        batch_size=cfg.batch_size.IC,
        lambda_weighting={"u": 1.0, "u__t": 1.0},
        parameterization={t_symbol: 0.0},
    )
    domain.add_constraint(IC, "IC")


    # boundary condition (source)
    BC = PointwiseBoundaryConstraint(
        nodes=nodes,
        geometry=geo,
        outvar={"u": 0},
        batch_size=cfg.batch_size.BC,
        parameterization=time_range,
    )
    domain.add_constraint(BC, "BC")


    # interior
    interior = PointwiseInteriorConstraint(
        nodes=nodes,
        geometry=geo,
        outvar={"wave_equation": 0},
        batch_size=cfg.batch_size.interior,
        parameterization=time_range,
    )
    domain.add_constraint(interior, "interior")


    # Read in npz files generated using finite difference simulator
    # add validation data
    wf_filename = to_absolute_path(f"FDM_data/wave_2d_source.npz")
    wave = np.load(wf_filename)
    t = np.array(wave["t_input"]).flatten()[:, None].astype(np.float32)
    y = np.array(wave["y_input"]).flatten()[:, None].astype(np.float32)
    x = np.array(wave["x_input"]).flatten()[:, None].astype(np.float32)
    u = np.array(wave["u_output"]).astype(np.float32)
    X, Y, T = np.meshgrid(x, y, t)
    X = np.expand_dims(X.flatten(), axis=-1)
    Y = np.expand_dims(Y.flatten(), axis=-1)
    T = np.expand_dims(T.flatten(), axis=-1)
    u = np.expand_dims(u.flatten(), axis=-1)
    invar_numpy = {"x": X, "y": Y, "t": T}
    outvar_numpy = {"u": u}

    validator = PointwiseValidator(
        nodes=nodes,
        invar=invar_numpy,
        true_outvar=outvar_numpy,
        batch_size=1024,
        plotter=ValidatorPlotter(),
    )
    domain.add_validator(validator)
    
    # add inferencer data
    grid_inference = PointwiseInferencer(
        nodes=nodes,
        invar=invar_numpy,
        output_names=["u"],
        batch_size=1024,
        plotter=InferencerPlotter(),
    )
    domain.add_inferencer(grid_inference, "inf_data")

    # make solver
    slv = Solver(cfg, domain)

    # start solver
    slv.solve()


if __name__ == "__main__":
    run()

I’m suspecting it has to do with sympy. Any suggestions on how to fix this?

8

Hi @cpe.sk

Yes, this is a sympy error at the light where you’re defining wave_source_outvar["s_func"]. Is there are reason you’re using sympy math operations here, could you change things to numpy operations like np.sin(*)? I’m assuming you’re calculating an array of target values here. Best to avoid mixing numpy and sympy together when possible, as they can lead to complications.

1 Like
9

It worked! Thank you!
The 2D validator plots have a significant error value. I’m trying to change the loss function type and the learning rate to fix this. Do you have any suggestions for that?

As for the parameterization I was trying to attempt, I thought about testing it for 1D first. As I added the parameterization for frequency. I got this error:

Error executing job with overrides: []
Traceback (most recent call last):
  File "wave_1d_param.py", line 213, in run
    parameterization={time_range, freq_range},
TypeError: unhashable type: 'dict'

Here’s my trial code:

import numpy as np
from sympy import Symbol, sin, exp, Eq, Number, Function

import modulus
from modulus.hydra import to_absolute_path, ModulusConfig, instantiate_arch
from modulus.solver import Solver
from modulus.domain import Domain
from modulus.geometry import Parameterization, Parameter
from modulus.geometry.primitives_1d import Line1D
from modulus.domain.constraint import (
    PointwiseBoundaryConstraint,
    PointwiseInteriorConstraint,
    PointwiseConstraint,
)

from modulus.domain.validator import PointwiseValidator
from modulus.domain.inferencer import PointwiseInferencer
from modulus.key import Key
from modulus.node import Node
#from wave_equation import WaveEquation1D
from modulus.eq.pdes.wave_equation import WaveEquation
from modulus.utils.io import (
    ValidatorPlotter,
    InferencerPlotter,
)
from modulus.eq.pde import PDE

class Wave1DEquation(PDE):
    """
    Wave equation

    Parameters
    ==========
    u : str
        The dependent variable.
    c : float, Sympy Symbol/Expr, str
        Wave speed coefficient. If `c` is a str then it is
        converted to Sympy Function of form 'c(x,y,z,t)'.
        If 'c' is a Sympy Symbol or Expression then this
        is substituted into the equation.
    dim : int
        Dimension of the wave equation (1, 2, or 3). Default is 2.
    time : bool
        If time-dependent equations or not. Default is True.
    mixed_form: bool
        If True, use the mixed formulation of the wave equation.

    Examples
    ========
    >>> we = WaveEquation(c=0.8, dim=3)
    >>> we.pprint()
      wave_equation: u__t__t - 0.64*u__x__x - 0.64*u__y__y - 0.64*u__z__z
    >>> we = WaveEquation(c='c', dim=2, time=False)
    >>> we.pprint()
      wave_equation: -c**2*u__x__x - c**2*u__y__y - 2*c*c__x*u__x - 2*c*c__y*u__y
    """

    name = "WaveEquation"

    def __init__(self, u="u", c="c", s_func="s_func", dim=3, time=True, mixed_form=False):
        # set params
        self.u = u
        self.dim = dim
        self.time = time
        self.mixed_form = mixed_form

        # coordinates
        x, y, z = Symbol("x"), Symbol("y"), Symbol("z")

        # time
        t = Symbol("t")

        # make input variables
        input_variables = {"x": x, "y": y, "z": z, "t": t}
        if self.dim == 1:
            input_variables.pop("y")
            input_variables.pop("z")
        elif self.dim == 2:
            input_variables.pop("z")
        if not self.time:
            input_variables.pop("t")

        # Scalar function
        assert type(u) == str, "u needs to be string"
        u = Function(u)(*input_variables)

        # wave speed coefficient
        if type(c) is str:
            c = Function(c)(*input_variables)
        elif type(c) in [float, int]:
            c = Number(c)
        
        # wave source function
        if type(s_func) is str:
            s_func = Function(s_func)(*input_variables)
        elif type(s_func) in [float, int]:
            s_func = Number(s_func)

        # set equations
        self.equations = {}

        if not self.mixed_form:
            self.equations["wave_equation"] = (
                u.diff(t, 2)
                - c**2 * u.diff(x, 2)
                - c**2 * u.diff(y, 2)
                - c**2 * u.diff(z, 2)
                - s_func
            )


@modulus.main(config_path="conf", config_name="config_param")
def run(cfg: ModulusConfig) -> None:

    # make geometry
    x, t_symbol = Symbol("x"), Symbol("t")
    L = 1.0
    geo = Line1D(0, L)
    time_range = {t_symbol: (0, L)}
    freq = Symbol("freq")
    freq_range = {freq: (4, 16)}
    param_ranges = {
        t_symbol: time_range,
        freq: freq_range,
    }

    pr = Parameterization(param_ranges)
    
    # make list of nodes to unroll graph on
    we = Wave1DEquation(u="u", c=1.0, s_func="s_func", dim=1, time=True)
    wave_net = instantiate_arch(
        input_keys=[Key("x"), Key("t")],
        output_keys=[Key("u")],
        cfg=cfg.arch.fully_connected,
    )
    
    source_net = instantiate_arch(
        input_keys=[Key("x"), Key("t"), Key("freq"),],
        output_keys=[Key("s_func")],
        cfg=cfg.arch.fully_connected,
    )
    nodes = we.make_nodes(detach_names=["s_func"]) + [
        wave_net.make_node(name="wave_network"),
        source_net.make_node(name="source_network"),
    ]

    # add constraints to solver
    

    # make domain
    domain = Domain()

    #Source function parameters
    locx1 = 0.0
    Amp = 1                                                                         #Amplitude of the source
    sigma = 0.2                                                                     #Width of the Gaussian impulse
    freq = 1 

    mesh_x, mesh_t, mesh_f = np.meshgrid(
        np.linspace(0, 5, 512), np.linspace(0, 5, 512), np.linspace(4, 16, 6), indexing="ij"
    )
    wave_source_invar = {}
    wave_source_invar["x"] = np.expand_dims(mesh_x.flatten(), axis=-1)
    wave_source_invar["t"] = np.expand_dims(mesh_t.flatten(), axis=-1)
    wave_source_invar["freq"] = np.expand_dims(mesh_t.flatten(), axis=-1)
    
    wave_source_outvar = {}
    wave_source_outvar["s_func"] = Amp * np.sin(4*np.pi*wave_source_invar["freq"]*wave_source_invar["t"]) * (\
        np.exp(-0.5*(((wave_source_invar["x"]-locx1)/sigma)**2))\
    )
        
    # add source points constraint
    source = PointwiseConstraint.from_numpy(
        nodes=nodes,
        invar=wave_source_invar,
        outvar=wave_source_outvar,
        batch_size=1024,
    )
    domain.add_constraint(source, "Source")
        
        
    # initial condition
    IC = PointwiseInteriorConstraint(
        nodes=nodes,
        geometry=geo,
        outvar={"u": 0, "u__t": 0},
        batch_size=cfg.batch_size.IC,
        lambda_weighting={"u": 1.0, "u__t": 1.0},
        parameterization={t_symbol: 0}, 
    )
    domain.add_constraint(IC, "IC")


    # right boundary condition (source)
#     BCr = PointwiseBoundaryConstraint(
#         nodes=nodes,
#         geometry=geo,
#         #outvar={"u": sin(4*freq*np.pi*t_symbol)},
#         outvar={"u": 0},
#         batch_size=cfg.batch_size.BC,
#         parameterization=time_range,
#         criteria=Eq(x, 0)
#     )
#     domain.add_constraint(BCr, "BCr")


    # left boundary condition
    BCl = PointwiseBoundaryConstraint(
        nodes=nodes,
        geometry=geo,
        outvar={"u": 0},
        batch_size=cfg.batch_size.BC,
        parameterization={time_range, freq_range},
        criteria=Eq(x, L)
    )
    domain.add_constraint(BCl, "BCl")

    # interior
    interior = PointwiseInteriorConstraint(
        nodes=nodes,
        geometry=geo,
        outvar={"wave_equation": 0},
        batch_size=cfg.batch_size.interior,
        parameterization={time_range, freq_range},
    )
    domain.add_constraint(interior, "interior")

    # Read in npz files generated using finite difference simulator
    # add validation data
    wf_filename = to_absolute_path(f"FDM_data/wave_1d_source.npz")
    wave = np.load(wf_filename)
    t = np.array(wave["t_input"]).flatten()[:, None].astype(np.float32)
    x = np.array(wave["x_input"]).flatten()[:, None].astype(np.float32)
    u = np.array(wave["u_output"]).astype(np.float32)
    X, T = np.meshgrid(x, t)
    X = np.expand_dims(X.flatten(), axis=-1)
    T = np.expand_dims(T.flatten(), axis=-1)
    u = np.expand_dims(u.flatten(), axis=-1)
    invar_numpy = {"x": X, "t": T}
    outvar_numpy = {"u": u}

    validator = PointwiseValidator(
        nodes=nodes,
        invar=invar_numpy,
        true_outvar=outvar_numpy,
        batch_size=128,
        plotter=ValidatorPlotter(),
    )
    domain.add_validator(validator)
    
    # add inferencer data
    grid_inference = PointwiseInferencer(
        nodes=nodes,
        invar=invar_numpy,
        output_names=["u"],
        batch_size=1024,
        plotter=InferencerPlotter(),
    )
    domain.add_inferencer(grid_inference, "inf_data")

    # make solver
    slv = Solver(cfg, domain)

    # start solver
    slv.solve()


if __name__ == "__main__":
    run()

Any ideas on how to get around this?

UPDATE: I managed to run the code successfully by changing the parameterization definition to:

#     time_range = {t_symbol: (0, L)}
    time_range =  (0, L)
    freq = Parameter("freq")
    freq_range = (4, 16)
    parameterization = Parameterization(
        {freq: freq_range, t_symbol: time_range}
    )

And then using parameterization in the interior and boundary constraints. This got the code running without errors, however, I’m facing two other challenges:

  • The results got a lot worse. So, I’m not sure if I’m implementing the parameterization correctly.
    validator_u
    validator_u1481Ă—412 47.5 KB
  • How to plot the results for different parameterization values separately? Is there a good example that I can start from?

Please advise.

10

@cpe.sk

Great glad you’ve got things running.

The results got a lot worse. So, I’m not sure if I’m implementing the parameterization correctly.

Yes, looks like the model is really having a tough time learning. Start with lowering your parameterization range and increasing the size of your model / training points. Get the training working for this simpler set up then start increasing your parameterization range.

How to plot the results for different parameterization values separately? Is there a good example that I can start from?

I’m not sure off the top of my head. If you want to plot more results during training consider adding a plotter object for your validation or inference constraint. The default plotter classes can be found here. This should allow you to create multiple plots to your tensorboard report. Worse case you can create multiple constraints for evaluation each with a fixed parameterization.

11

Start with lowering your parameterization range and increasing the size of your model / training points. Get the training working for this simpler set up then start increasing your parameterization range.

I lowered the freq_range from (4,16) to (4,5) and increased the training points. Unfortunately, the results did not improve. I tried other solutions like testing another architecture model, fused_fully_connected and I got:

File "C:\ProgramData\Anaconda3\envs\modulus_env\lib\site-packages\modulus-22.9-py3.8.egg\modulus\models\fused_mlp.py", line 24, in <module>
    import tinycudann as tcnn
ModuleNotFoundError: No module named 'tinycudann'

So I used this command from Issues · NVlabs/nvdiffrec · GitHub

pip install --global-option="--no-networks" git+https://github.com/NVlabs/tiny-cuda-nn#subdirectory=bindings/torch

and tried running again:

 File "C:\ProgramData\Anaconda3\envs\modulus_env\lib\site-packages\tinycudann\modules.py", line 207, in __init__
    raise RuntimeError(f"Cannot create `NetworkWithInputEncoding` because tiny-cuda-nn was not compiled with neural network support.")
RuntimeError: Cannot create `NetworkWithInputEncoding` because tiny-cuda-nn was not compiled with neural network support.

UPDATE:
I also tried using another arch model, the multiscale_fourier and got this error:

File "C:\ProgramData\Anaconda3\envs\modulus_env\lib\site-packages\modulus-22.9-py3.8.egg\modulus\models\multiscale_fourier_net.py", line 114, in __init__
    layers.FourierLayer(
  File "C:\ProgramData\Anaconda3\envs\modulus_env\lib\site-packages\modulus-22.9-py3.8.egg\modulus\models\layers\fourier_layers.py", line 18, in __init__
    if isinstance(frequencies[0], str):
TypeError: 'int' object is not subscriptable

Can you help me get around this?
Thank you so much for your support.

12

Hi @cpe.sk

The Fused Neural networks are highly optimized but not as robust (from an installation stand point) as the other Pytorch models. It will require TinyCUDADNN to be built on your machine for your specific CUDA version. I would advise against using this if possible when just prototyping / testing.

For multiscale_fourier error, what does your code look like for constructing the model? Seems its complaining about the frequencies parameter being an interger when it needs to be something like a list or tuple. Please have a look at the API documentation for guidance on model parameters.

13

This topic was automatically closed 14 days after the last reply. New replies are no longer allowed.