I want to couple solve an ODE and a PDE in a system in which the heat-diffusion equation is the PDE and the ODE is for concentration (c). My inputs are t, x and my outputs are T and c. My problem is that the loss for my convective boundary loss doesn’t decrease and if I give a learning rate higher than 1e-4, I get NaNs. I don’t know what my problem is. Thanks. Below is My full code for this problem.

```
using NeuralPDE, Lux, ModelingToolkit, Optimization, OptimizationOptimJL
import ModelingToolkit: Interval, infimum, supremum
@parameters t, x
@variables T(..) , c(..)
Dt = Differential(t)
#Dyy = Differential(y)^2
Dx = Differential(x)
#Dy = Differential(y)
Dxx = Differential(x)^2
########################
Le = 0.013 # m
He = 0.065 # m
kl = 0.7 # Thermal conductivity in radial direction (W/mk)
T0 = 273.15 + 25 # Initial Temperature (C)
ρ = 2231.2 # Cell density (kg/m^3)
Cp = 1100 # Heat capacity (J/kgK)
h = 12 # Heat transfer coefficient (W/m^2K)
Rc = 8.314 # Gas constant (J/molK)
E = 2.7e5 # Activation energy value of electrolyte [J/mol]
A = 5.14e25 # Reaction factor of electrolyte [1/s]
m = 1 # Reaction order of electrolyte
c0 = 1 # Initial dimensionless concentration
t_max = 133 * 60 # solution time
H = 6.2e5
W = 335
Tamb = 420
####################################################
eqs = [ρ*Cp*Dt(T(t,x)) ~ kl * Dxx(T(t, x)) - H*W*Dt(c(t,x)),
-Dt(c(t,x)) ~ A* c(t,x) * exp(-E/(Rc*T(t,x)))
]
bcs = [ T(0,x) ~ T0,
c(0,x) ~ c0,
Dx(T(t,0)) ~ 0,
h * (T(t, Le) - Tamb) - kl * Dx(T(t, 0)) ~ 0
]
# Space and time domains
domains = [t ∈ Interval(0, t_max),
x ∈ Interval(0, Le),]
# Neural network
input_ = length(domains)
n = 20
chain = [Lux.Chain(Dense(input_, n, Lux.relu), Dense(n, n, Lux.relu), Dense(n, n, Lux.relu), Dense(n, n, Lux.relu), Dense(n, n, Lux.relu), Dense(n, n, Lux.relu), Dense(n, n, Lux.relu), Dense(n, 1)) for _ in 1:2]
@named pdesystem = PDESystem(eqs, bcs, domains, [t,x], [T(t,x),c(t,x)])
strategy = NeuralPDE.QuadratureTraining()
discretization = PhysicsInformedNN(chain, strategy)
sym_prob = NeuralPDE.symbolic_discretize(pdesystem, discretization)
pde_loss_functions = sym_prob.loss_functions.pde_loss_functions
bc_loss_functions = sym_prob.loss_functions.bc_loss_functions
callback = function (p, l)
println("loss: ", l)
println("pde_losses: ", map(l_ -> l_(p), pde_loss_functions))
println("bcs_losses: ", map(l_ -> l_(p), bc_loss_functions))
return false
end
loss_functions = [pde_loss_functions; bc_loss_functions]
function loss_function(θ, p)
sum(map(l -> l(θ), loss_functions))
end
f_ = OptimizationFunction(loss_function, Optimization.AutoZygote())
prob = Optimization.OptimizationProblem(f_, sym_prob.flat_init_params)
using Optimization
using OptimizationOptimisers
res = Optimization.solve(prob, Adam(0.0001); callback = callback, maxiters = 100000)
```