# Failed Unrolling Graph with custom PDE

I am trying to solve a set of PDEs and algebraic equations with Modulus. Following is the class which defines the PDE node.

``````### P2D.py

from sympy import Symbol, Function, Number
from sympy import exp, log, Pow

from modulus.eq.pde import PDE
from modulus.node import Node

class Pseudo2D(PDE):
"""
Pseudo 2D model for battery

Parameters
==========
D : float, Sympy Symbol/Expr, str
Diffusion coefficient of electrolyte. If `D` is a str then it is
converted to Sympy Function of form `D(x,y,z,t)`.
If `D` is a Sympy Symbol or Expression then this
is substituted into the equation. This allows for
variable diffusion coefficient.

D_s : float, Sympy Symbol/Expr, str
Diffusion coefficient of lithium in solid matrix

ep : float, Sympy Symbol/Expr, str
Porosity of electrode

t_plus : float, Sympy Symbol/Expr, str
Transference number

a : float, Sympy Symbol/Expr, str
Specific interacial area (m^2/m^3)

sigma : float, Sympy Symbol/Expr, str
Conductivity of solid matrix

kappa : float, Sympy Symbol/Expr, str
Conductivity of electrolyte

f_A : float, Sympy Symbol/Expr, str
Activity coefficient of salt

U : float, Sympy Symbol/Expr, str
Open circuit voltage

c_t : float, Sympy Symbol/Expr, str
Total concentration

I : float, Sympy Symbol/Expr, str
Total current density

time : bool
If time-dependent equations or not. Default is True.

name = "Pseudo2D"

def __init__(self, ep =1, D=0.001, D_s=0.001, t_plus=1, a=1, sigma=1, kappa=1, f_A=1, U=1, c_t=1, U_0=1, I=1.0, T_0=273, M = 1, time=True):
# set params
#self.dim = dim
self.time = time
#self.mixed_form = mixed_form

#Constant parameters
F = Number(96500) # Faraday's constant
R = Number(8.314) # Universal gas constant
C_p = Number(4.18) # Specific heat capacity

# coordinates
x, r  = Symbol("x"), Symbol("r")

# time
t = Symbol("t")

# make input variables
input_variables = {"x": x, "r": r, "t": t}

if not self.time:
input_variables.pop("t")

# Concentration variables
c = Function("c")(*input_variables)
c_s = Function("c_s")(*input_variables)

# Current densities
i_1 = Function("i_1")(*input_variables)
i_2 = Function("i_2")(*input_variables)

# Potentials
Phi_1 = Function("Phi_1")(*input_variables)
Phi_2 = Function("Phi_2")(*input_variables)

# Pore-wall flux
j_n = Function("j_n")(*input_variables)

# Temperature
#T = Function("T")(*input_variables)
T = Number(300)  #Assume constant temperature

# Porosity
if isinstance(ep, str):
ep = Function(ep)(*input_variables)
elif isinstance(ep, (float, int)):
ep = Number(ep)

# Diffusion coefficient
if isinstance(D, str):
D = Function(D)(*input_variables)
if isinstance(D_s, str):
D_s = Function(D_s)(*input_variables)
elif isinstance(D_s, (float, int)):
D_s = Number(D_s)

# Transference coefficient
if isinstance(t_plus, str):
t_plus = Function(t_plus)(*input_variables)
elif isinstance(t_plus, (float, int)):
t_plus = Number(t_plus)

# Interfacial area
if isinstance(a, str):
a = Function(a)(*input_variables)
elif isinstance(a, (float, int)):
a = Number(a)

# Conductivitty coefficients
if isinstance(sigma, str):
sigma = Function(sigma)(*input_variables)   # Conductivitty of solid matrix
elif isinstance(sigma, (float, int)):
sigma = Number(sigma)
if isinstance(sigma, str):
kappa = Function(kappa)(*input_variables)   # Conductivitty of electrolyte
elif isinstance(sigma, (float, int)):
kappa = Number(kappa)

# Activity coefficient
if isinstance(f_A, str):
f_A = Function(f_A)(*input_variables)
elif isinstance(f_A, (float, int)):
f_A = Number(f_A)

# Open circuit voltage
if isinstance(U, str):
U = Function(U)(*input_variables)
elif isinstance(U, (float, int)):
U = Number(U)
if isinstance(U_0, str):
U_0 = Function(U_0)(*input_variables)
elif isinstance(U_0, (float, int)):
U_0 = Number(U_0)

if isinstance(c_t, str):
c_t = Function(c_t)(*input_variables)
elif isinstance(c_t, (float, int)):
c_t = Number(c_t)

if isinstance(I, str):
I = Function(I)(*input_variables)
elif isinstance(I, (float, int)):
I = Number(I)

if isinstance(T_0, str):
T_0 = Function(T_0)(*input_variables)

elif isinstance(T_0, (float, int)):
T_0 = Number(T_0)

if isinstance(M, str):
M = Function(M)(*input_variables)
elif isinstance(M, (float, int)):
M = Number(M)

e = Number(2.71828)

# set equations
self.equations = {}
self.equations["Material_balance"] = ep * c.diff(t) -(D * c.diff(x)).diff(x) + i_2 * t_plus.diff(x)/F + a * j_n * (1-t_plus)

self.equations["Insertion_material_current"] =  (
(i_1
+ sigma * Phi_1.diff(x))
)

self.equations["Total_current"] = I - i_1 - i_2

self.equations["Electrolyte_current"] = i_2 + kappa * Phi_2.diff(x) - kappa * R * T *(1+ c* (log(f_A)).diff(c) )* (1- t_plus) * log(c).diff(x)/F

self.equations["Pore_wall_flux"] = a * j_n - i_2.diff(x)/F

self.equations["Spherical_particles"] = c_s.diff(t) - D_s * (c_s.diff(r,2) + (2/r) * c_s.diff(r))

self.equations["Butler_Volmer_kinetics"] = j_n - kappa * c**(0.5) * (c_t - c_s) * c_s**(0.5) *(Pow(e,(Phi_1-Phi_2 - U)*F/(2*R*T)) - Pow(e,-(Phi_1-Phi_2 - U)*F/(2*R*T))   )

``````

The source code is

``````# Battery.py

import numpy as np
from sympy import Symbol, sin

import math

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

from modulus.domain.validator import PointwiseValidator
from modulus.key import Key
from modulus.node import Node
#from wave_equation import WaveEquation1D
from diffusion import Diffusion

from P2D import Pseudo2D

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

# make list of nodes to unroll graph on
p2d = Pseudo2D()
p2d.pprint()
battery_net = instantiate_arch(
input_keys=[Key("x"), Key("r"), Key("t")],
output_keys=[Key("c"), Key("i_1"),Key("i_2"),Key("j_n"),Key("Phi_1"),Key("Phi_2"),Key("c_s")],
cfg=cfg.arch.fully_connected,
)
nodes = p2d.make_nodes() + [battery_net.make_node(name="battery_network")]

# make geometry
x, r, t_symbol = Symbol("x"),Symbol("r"), Symbol("t")
L = float(np.pi)
time_range = {t_symbol: (0, 2 * L)}

height = 0.1
width = 0.1
rec = Rectangle((-width / 2, -height / 2), (width / 2, height / 2))

# make domain
domain = Domain()

# interior
interior = PointwiseInteriorConstraint(
nodes=nodes,
geometry=rec,
outvar={"Material_balance": 0, "Insertion_material_current": 0, "Total_current":0,"Electrolyte_current": 0, "Pore_wall_flux":0,"Spherical_particles":0,"Butler_Volmer_kinetics":0 },
#outvar={"Material_balance": 0},
batch_size=cfg.batch_size.Interior,
parameterization=time_range,
)

if __name__ = '__main__':
run()
``````

Running Battery.py results in the following error

``````[04:46:46] - JIT using the NVFuser TorchScript backend
[04:46:46] - JitManager: {'_enabled': True, '_arch_mode': <JitArchMode.ONLY_ACTIVATION: 1>, '_use_nvfuser': True, '_autograd_nodes': False}
[04:46:46] - GraphManager: {'_func_arch': False, '_debug': False, '_func_arch_allow_partial_hessian': True}
Material_balance: c__t - 0.001*c__x__x
Insertion_material_current: i_1 + Phi_1__x
Total_current: -i_1 - i_2 + 1.0
Electrolyte_current: i_2 + Phi_2__x
Pore_wall_flux: j_n - i_2__x/96500
Spherical_particles: -0.001*c_s__r__r + c_s__t - 0.002*c_s__r/r
Butler_Volmer_kinetics: -(1 - c_s)*(-2.71828**(-19.3448801218828*Phi_1 + 19.3448801218828*Phi_2 + 19.3448801218828) + 2.71828**(19.3448801218828*Phi_1 - 19.3448801218828*Phi_2 - 19.3448801218828))*c**0.5*c_s**0.5 + j_n
####################################
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: [Material_balance, Insertion_material_current, Total_current, Electrolyte_current, Pore_wall_flux, Spherical_particles, Butler_Volmer_kinetics]
computable var: [x, y, sdf, area, t, Material_balance]
####################################
Nodes in graph:
node: Sympy Node: Material_balance
evaluate: SympyToTorch
inputs: []
derivatives: [c__t, c__x__x]
outputs: [Material_balance]
optimize: False
node: Sympy Node: Insertion_material_current
evaluate: SympyToTorch
inputs: [i_1]
derivatives: [Phi_1__x]
outputs: [Insertion_material_current]
optimize: False
node: Sympy Node: Total_current
evaluate: SympyToTorch
inputs: [i_1, i_2]
derivatives: []
outputs: [Total_current]
optimize: False
node: Sympy Node: Electrolyte_current
evaluate: SympyToTorch
inputs: [i_2]
derivatives: [Phi_2__x]
outputs: [Electrolyte_current]
optimize: False
node: Sympy Node: Pore_wall_flux
evaluate: SympyToTorch
inputs: [j_n]
derivatives: [i_2__x]
outputs: [Pore_wall_flux]
optimize: False
node: Sympy Node: Spherical_particles
evaluate: SympyToTorch
inputs: [r]
derivatives: [c_s__r, c_s__r__r, c_s__t]
outputs: [Spherical_particles]
optimize: False
node: Sympy Node: Butler_Volmer_kinetics
evaluate: SympyToTorch
inputs: [Phi_1, Phi_2, c, c_s, j_n]
derivatives: []
outputs: [Butler_Volmer_kinetics]
optimize: False
node: Arch Node: battery_network
evaluate: FullyConnectedArch
inputs: [x, r, t]
derivatives: []
outputs: [c, i_1, i_2, j_n, Phi_1, Phi_2, c_s]
optimize: True
####################################
Error executing job with overrides: []
Traceback (most recent call last):
File "Battery.py", line 52, in run
interior = PointwiseInteriorConstraint(
File "/modulus/modulus/domain/constraint/continuous.py", line 493, in __init__
super().__init__(
File "/modulus/modulus/domain/constraint/constraint.py", line 55, in __init__
self.model = Graph(
File "/modulus/modulus/graph.py", line 94, in __init__
raise RuntimeError("Failed Unrolling Graph")
RuntimeError: Failed Unrolling Graph
``````

For instance, the second PDE (“Insertion_material_current”) expects ‘i_1’ as input. But ‘i_1’ is a dependent variable, the value of which is to be determined by solving the set of equations. Similar issue is observed for the “Material_balance” equation if the divergence term is made non-linear by multiplying variable ‘c’ to it.

Is there a way to get around this problem in Modulus without creating additional PDE nodes for each of the equations.

As you mentioned the issue is that you have no nodes that can compute many of the intermediate variables (`i_1`, `i_2`, and `j_n`) which causes the graph to fail. You need to create some node to compute these variables.

You could do another PDE, although this is likely over kill. If its a neural network, you can use the `make_node()` function. If these are just some analytical function, the `Node.from_sympy()` function is a good option. There’s come example use cases of this function in the turbulent channel example.

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