# Second Fick Law -- NeuralPDE

**URL:** <https://discourse.julialang.org/t/second-fick-law-neuralpde/99497>\
**Category:** Machine Learning\
**Tags:** question\
**Created:** [May 27, 2023, 11:00pm UTC](https://discourse.julialang.org/t/second-fick-law-neuralpde/99497 "2023-05-27T23:00:01Z")\
**Posts on this page:** 10\
**Page:** 1

<div class="post-metadata">

**Author:** ![aquarelleX332](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/aquarellex332/32/21048_2.png) [@aquarelleX332](https://discourse.julialang.org/u/aquarelleX332)\
**Post date:** [May 27, 2023, 11:00pm UTC](https://discourse.julialang.org/t/second-fick-law-neuralpde/99497/1 "2023-05-27T23:00:02Z")

</div>

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](https://neuralpde.sciml.ai/stable/examples/wave/)) so I tried this:

```julia
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:

 ![Fick](https://global.discourse-cdn.com/julialang/original/3X/4/7/47154581b7c5522cc0e4b55360b5bd076f0d9424.jpeg)  
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](https://neuralpde.sciml.ai/stable/tutorials/gpu/) and obtained the following error:

```julia
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:

 ![Fick_PINN](https://global.discourse-cdn.com/julialang/original/3X/9/3/930d8e6c5ff39128f8c16ad8e22fea19e7331bcb.jpeg)

---

<div class="post-metadata">

**Author:** ![avikpal](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/avikpal/32/6550_2.png) [@avikpal](https://discourse.julialang.org/u/avikpal)\
**Post date:** [May 28, 2023, 1:54pm UTC](https://discourse.julialang.org/t/second-fick-law-neuralpde/99497/2 "2023-05-28T13:54:01Z")

</div>

> [@aquarelleX332](#):
>
> Q3: As I have an Nvidia 3070 Ti graphic card, I tried to use CUDA. But before doing so, I tested the [following code](https://neuralpde.sciml.ai/stable/tutorials/gpu/) and obtained the following error:

Use [https://github.com/SciML/NeuralPDE.jl/blob/master/docs/src/tutorials/gpu.md](https://github.com/SciML/NeuralPDE.jl/blob/master/docs/src/tutorials/gpu.md). I think the latest version of the code hasn’t been deployed to the documentation website.

---

<div class="post-metadata">

**Author:** ![MilkshakeForReal](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/milkshakeforreal/32/32861_2.png) [@MilkshakeForReal](https://discourse.julialang.org/u/MilkshakeForReal)\
**Post date:** [May 28, 2023, 6:44pm UTC](https://discourse.julialang.org/t/second-fick-law-neuralpde/99497/3 "2023-05-28T18:44:37Z")

</div>

> [@aquarelleX332](#):
>
> Q1: Do you think the way the initial condition is presented is fair (i.e. in 2 lines)?

```julia
u₀(x) = ifelse(x < 0.5, 2x, 2(1-x))
@register_symbolic u₀(x)

bcs = [u(t,0) ~ 0.,       
       u(t,1) ~ 0.,
       u(0,x) ~u₀(x) ]

```

---

<div class="post-metadata">

**Author:** ![aquarelleX332](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/aquarellex332/32/21048_2.png) [@aquarelleX332](https://discourse.julialang.org/u/aquarelleX332)\
**Post date:** [May 28, 2023, 9:04pm UTC](https://discourse.julialang.org/t/second-fick-law-neuralpde/99497/4 "2023-05-28T21:04:02Z")

</div>

Thank you very much, Yicheng Wu. 👍 😉  
The way you have expressed the initial condition is perfect and elegant. The loss function runs much better (loss \< about 2.04e-4).  
So, this resolves Q1 and Q2.

---

<div class="post-metadata">

**Author:** ![aquarelleX332](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/aquarellex332/32/21048_2.png) [@aquarelleX332](https://discourse.julialang.org/u/aquarelleX332)\
**Post date:** [May 28, 2023, 9:13pm UTC](https://discourse.julialang.org/t/second-fick-law-neuralpde/99497/5 "2023-05-28T21:13:12Z")

</div>

Thank you very much avikpal. 👍 😉  
Everything works perfectly as well concerning the [code you mentionned](https://github.com/SciML/NeuralPDE.jl/blob/master/docs/src/tutorials/gpu.md) and for the modification I incorporated in my code (second Fick Law).

---

<div class="post-metadata">

**Author:** ![MilkshakeForReal](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/milkshakeforreal/32/32861_2.png) [@MilkshakeForReal](https://discourse.julialang.org/u/MilkshakeForReal)\
**Post date:** [May 28, 2023, 9:19pm UTC](https://discourse.julialang.org/t/second-fick-law-neuralpde/99497/6 "2023-05-28T21:19:43Z")

</div>

Not an ad but you can try using [Sophon.jl,](https://github.com/YichengDWu/Sophon.jl) it should be faster and more accurate. NeuralPDE is being rewritten, its old implementation has some shortcomings. Sophon is a rewritten version of NeuralPDE.

---

<div class="post-metadata">

**Author:** ![aquarelleX332](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/aquarellex332/32/21048_2.png) [@aquarelleX332](https://discourse.julialang.org/u/aquarelleX332)\
**Post date:** [May 28, 2023, 9:49pm UTC](https://discourse.julialang.org/t/second-fick-law-neuralpde/99497/7 "2023-05-28T21:49:13Z")

</div>

I’m already looking at the GitHub pages.  
Thanks for getting my attention. Everything is going so fast with Julia now.

---

<div class="post-metadata">

**Author:** ![aquarelleX332](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/aquarellex332/32/21048_2.png) [@aquarelleX332](https://discourse.julialang.org/u/aquarelleX332)\
**Post date:** [May 29, 2023, 6:50pm UTC](https://discourse.julialang.org/t/second-fick-law-neuralpde/99497/9 "2023-05-29T18:50:16Z")

</div>

After trying sucessfully the code recommanded by Avik Pal:  
[NeuralPDE.jl/gpu.md at master · SciML/NeuralPDE.jl · GitHu](https://github.com/SciML/NeuralPDE.jl/blob/master/docs/src/tutorials/gpu.md)  
I modified my program as this:

```julia
using NeuralPDE, Lux, CUDA, Random, ComponentArrays
using Optimization, OptimizationOptimisers
using Plots
import ModelingToolkit: Interval

@parameters t, x
@variables u(..)
Dxx = Differential(x)^2
Dt = Differential(t)

#1D-PDE
κ = 0.75 # Fick diffusion coefficient
eq = Dt(u(t,x)) ~ κ*Dxx(u(t,x)) # Second Fick law

# Initial and boundary conditions
u₀(x) = ifelse(x < 0.5, 2x, 2(1-x))
@register_symbolic u₀(x)
bcs = [u(t,0) ~ 0., # for all t > 0 -- BC   
       u(t,1) ~ 0., # for all t > 0 -- BC
       u(0,x) ~ u₀(x)] # for all 0 ≤ x ≤ 1/2 and for t > 0 -- Init.Cond
                                       # 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)
ps = Lux.setup(Random.default_rng(), chain)[1]
ps = ps |> ComponentArray |> gpu .|> Float64
discretization = PhysicsInformedNN(chain, strategy, init_params = ps)

@named pde_system = PDESystem(eq, bcs, domains, [t,x], [u(t,x)])
prob = discretize(pde_system, discretization)
symprob = symbolic_discretize(pde_system, discretization)

callback = function (p,l)
    println("Loss: $l")
    return false
end

# optimizer
#opt = OptimizationOptimJL.BFGS()
res = Optimization.solve(prob, Adam(0.001); callback = callback, maxiters=3000)
prob = remake(prob, u0 = res.u)
res = Optimization.solve(prob, Adam(0.001); 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, xaxis = "t (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)

```

However I get this error:

```julia
ERROR: CuArray only supports element types that are allocated inline.
Real is not allocated inline

```

I cannot detect my error, unfortunately.  
Thank you in advance for any help.

---

<div class="post-metadata">

**Author:** ![aquarelleX332](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/aquarellex332/32/21048_2.png) [@aquarelleX332](https://discourse.julialang.org/u/aquarelleX332)\
**Post date:** [May 29, 2023, 9:42pm UTC](https://discourse.julialang.org/t/second-fick-law-neuralpde/99497/10 "2023-05-29T21:42:03Z")

</div>

The step where the ERROR is thrown is:

```julia
prob = discretize(pde_system, discretization)

```

Up to the line as follows included

```julia
@named pde_system = PDESystem(eq, bcs, domains, [t,x], [u(t,x)])

```

REPL indicates::

```julia
PDESystem
Equations: Equation[Differential(t)(u(t, x)) ~ 0.75Differential(x)(Differential(x)(u(t, x)))]
Boundary Conditions: Equation[u(t, 0.0) ~ 0.0, u(t, 1.0) ~ 1.2246467991473532e-16exp(-7.4022033008170185t), u(0.0, 0) ~ sin(πx)]
Domain: Symbolics.VarDomainPairing[Symbolics.VarDomainPairing(t, 0.0..1.0), Symbolics.VarDomainPairing(x, 0.0..1.0)]
Dependent Variables: Num[u(t, x)]
Independent Variables: Num[t, x]
Parameters: SciMLBase.NullParameters()
Default Parameter ValuesDict{Any, Any}()

```

Then?..

---

<div class="post-metadata">

**Author:** ![aquarelleX332](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/aquarellex332/32/21048_2.png) [@aquarelleX332](https://discourse.julialang.org/u/aquarelleX332)\
**Post date:** [June 3, 2023, 8:37pm UTC](https://discourse.julialang.org/t/second-fick-law-neuralpde/99497/11 "2023-06-03T20:37:34Z")

</div>

I checked the pilote of my Nvida card and select 530 instead of 525. This resolves the problem.