Hello everyone,
I’m trying to solve Fick’s second law (1D), which is analogous to the thermal diffusion equation (or heat equation). Fick’s equation explains how a solute diffuses in a solvent.
I was inspired by this code ( 1D Wave Equation with Dirichlet boundary conditions) so I tried this:
using NeuralPDE, Lux
using Optimization, OptimizationOptimJL, Plots
import ModelingToolkit: Interval
@parameters t, x
@variables u(..)
Dxx = Differential(x)^2
Dt = Differential(t)
#1D-PDE
κ = 1 # Fick diffusion coefficient
eq = Dt(u(t,x)) ~ κ*Dxx(u(t,x)) # Fick second law
# Initial and boundary conditions
bcs = [u(t,0) ~ 0., # for all t > 0 -- BC
u(t,1) ~ 0., # for all t > 0 -- BC
u(0,x) ~ 2x, # for all 0 ≤ x ≤ 1/2 and for t > 0 -- Init.Cond
u(0,x) ~ 2*(1. - x)] # for all 1/2 ≤ x ≤ 1 and for t > 0 -- Init.Cond
# Space and time domains
domains = [t ∈ Interval(0.0,1.0),
x ∈ Interval(0.0,1.0)]
# Discretization
dx = 0.01
# Neural network
inner = 20
chain = Lux.Chain(Dense(2,inner,Lux.σ),Dense(inner,inner,Lux.σ),Dense(inner,1))
strategy = GridTraining(dx)
discretization = PhysicsInformedNN(chain, strategy)
@named pde_system = PDESystem(eq,bcs,domains,[t,x],[u(t,x)])
prob = discretize(pde_system,discretization)
callback = function (p,l)
println("Loss: $l")
return false
end
# optimizer
opt = OptimizationOptimJL.BFGS()
res = Optimization.solve(prob,opt; callback = callback, maxiters=3000)
phi = discretization.phi
# Graphs
ts,xs = [infimum(d.domain):dx:supremum(d.domain) for d in domains]
analytic_sol_f(t,x) = sum([exp(-t*(n*π)^(2)) * sin(n*π*x) for n in 1:2:50000])
u_predict = reshape([first(phi([t,x],res.u)) for t in ts for x in xs],(length(ts),length(xs)))
u_analyt = reshape([analytic_sol_f(t,x) for t in ts for x in xs], (length(ts),length(xs)))
diff_u = abs.(u_predict .- u_analyt)
p1 = plot(ts, xs, u_analyt, linetype=:contourf, yaxis = "x (u.a.)", title = "Analytic");
p2 = plot(ts, xs, u_predict, linetype=:contourf, xaxis = "t (u.a.)", yaxis = "x (u.a.)", title = "Prediction");
p3 = plot(ts, xs, diff_u, linetype=:contourf, title = "Error");
plot(p1,p2,p3)
The analytical solution is:
where Cₙ and κ can be normalized to 1, L = xₘₐₓ = 1, 0.0 ≤ t ≤ 1.0, 0.0 ≤ x ≤ 1.0 and n is 1:2:50000
There are 3 questions:
- Q1: Do you think the way the initial condition is presented is fair (i.e. in 2 lines)?
- Q2: Does the comparison between the “Analytic” graph and the “Prediction” graph seems reasonable to you, considering the “Error” graph?
- Q3: As I have an Nvidia 3070 Ti graphic card, I tried to use CUDA. But before doing so, I tested the following code and obtained the following error:
ERROR: UndefVarError: `ComponentArray` not defined
Can you help me solve this error? Because when I adapt my code in order to use GPUs, I get the same error as reported above.
Thank you in advance for any help.
T.
Here the graphs obtained via the code as above:
