Conserve integral in a PINN

[ste]

I have the following code where I try to solve a little toy system of PDEs:

using ModelingToolkit, NeuralPDE, Flux, Optim, GalacticOptim, Plots, DiffEqSensitivity, OptimizationOptimisers, OptimizationOptimJL

import ModelingToolkit: Interval, infimum, supremum

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

eqn1 = Dt(ρ(t, x)) * Dxx(ρ(t, x) * Dx(u(t, x))) ~ 0 
eqn2 = sqrt(ρ(t, x)) ~ 0
eqns = [eqn1, eqn2]

L = 10.0
domains = [x ∈ Interval(-L/2, L/2),
           t ∈ Interval(0.0, 1.0)]

ρ_IC = ρ(0, x) ~ exp(-0.9 * (x - 1)^2)
u_IC = u(0, x) ~ 0.0
bcs = [ρ_IC, u_IC]

af = gelu
chain = [Chain(Dense(2, 16, af), Dense(16, 16, af), Dense(16, 1, softplus)),
            Chain(Dense(2, 16, af), Dense(16, 16, af), Dense(16, 1, softplus))]

discretization = PhysicsInformedNN(chain, QuadratureTraining())
@named pde_sys = PDESystem(eqns, bcs, domains, [t, x], [ρ(t, x), u(t, x)])
prob = NeuralPDE.discretize(pde_sys, discretization)

cb = function (p, l)
    println("Current loss is: $l")
    return false
end

res = Optimization.solve(prob, ADAM(0.1); callback = cb, maxiters = 100)

1. \rho should always be positive as in eqn2 a square-root is applied. I fixed this by using a softmax in the last layer instead of a linear activation function. It feels a bit strange to do it like this. Is there a way to have \rho positive as a constraint?
2. Similarly, I want \int_{L/2}^{L/2} \rho dx = 1 i.e. the integral of \rho, a probability density, should be normalized to 1 and then kept that way (thereby also automatically fulfilling 1)) at all times. How can I do that in SciML or something?

I am new to both Julia and PINNs, I hope my question makes sense.

Edit: I’m aware of this that seems quite close to what I want but the explanation of norm_loss_function is not very clear to me, especially since I have a system of 2 coupled PDEs.

[ChrisRackauckas]

@ste

ρ\rho should always be positive as in eqn2 a square-root is applied. I fixed this by using a softmax in the last layer instead of a linear activation function. It feels a bit strange to do it like this. Is there a way to have \rhoρ\rho positive as a constraint?

That’s a good way to do it. Or abs. Just make a positive function approximator and it’s fine.

@ste

Similarly, I want \int_{L/2}^{L/2} \rho dx = 1∫L/2L/2ρdx=1\int_{L/2}^{L/2} \rho dx = 1 i.e. the integral of \rhoρ\rho, a probability density, should be normalized to 1 and then kept that way (thereby also automatically fulfilling 1)) at all times. How can I do that in SciML or something?

Just add that as a 3rd equation.