# Could not unroll graph!

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,
)

# 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,
)

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

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 )
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(
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?

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

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!

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

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.

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.

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

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
)

# 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},
)

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

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

# Read in npz files generated using finite difference simulator
# add validation data
wf_filename = to_absolute_path(f"FDM_data/wave_2d_source.npz")
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(),
)

# add inferencer data
grid_inference = PointwiseInferencer(
nodes=nodes,
invar=invar_numpy,
output_names=["u"],
batch_size=1024,
plotter=InferencerPlotter(),
)

# 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?

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

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,
)

# 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},
)

# 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)
#     )

# 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)
)

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

# Read in npz files generated using finite difference simulator
# add validation data
wf_filename = to_absolute_path(f"FDM_data/wave_1d_source.npz")
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(),
)

# add inferencer data
grid_inference = PointwiseInferencer(
nodes=nodes,
invar=invar_numpy,
output_names=["u"],
batch_size=1024,
plotter=InferencerPlotter(),
)

# 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.
• How to plot the results for different parameterization values separately? Is there a good example that I can start from?

@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.

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.

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.

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