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.