Failed Unrolling Graph

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")]

    # add constraints to solver
    # 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,
    )
    #domain.add_validator(validator)




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.

Hi @shubhamsp2195

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.