Unable to handle non-linear convection terms

[v-chau]

I am trying to solve the equations for a fully developed laminar flow in a channel (continuity and x-momentum NSE equation) but my PINN code is unable to drive down the non-linear term ((rho/mu)udu/dx) to zero. The predicted solution matches the analytical when I replace (rho/mu)*u with a constant, say, K1 = 1, since the continuity equation is du/dx = 0. Does anyone know what seems to be the issue? Here is my code.

# SciML PINN technique for solving PDEs using NeuralPDE.jl package. https://neuralpde.sciml.ai/stable/
using NeuralPDE, Flux, ModelingToolkit, GalacticOptim, Optim, DiffEqFlux, Zygote, LaTeXStrings
import ModelingToolkit: Interval, infimum, supremum
using Plots, Printf, Quadrature,Cubature

# Define the independent and dependent variables, and the functional operators
@parameters x, y
@variables u(..)
Dx = Differential(x)
Dy = Differential(y)
Dxx = Differential(x)^2
Dyy = Differential(y)^2

# # Constants 
rho = 1.145
mu = 1.895e-5 
nu = mu/rho
k = 0.02625
alpha = 2.277e-5
H = 0.1 # channel height
L = 1 # channel length
ReH = 1000 # based on channel height
U = (ReH*nu)/H # Mean velocity
delP = (-12*U*mu)/(H^2) # pressure gradient
adam_max = 100
bfgs_max = 1000
K1 = 1

# Fully developed laminar flow in a channel
eqs = [Dx(u(x,y)) ~ 0, Dxx(u(x,y)) + Dyy(u(x,y)) ~ (delP/mu) + ((rho/mu)*u(x,y)*Dx(u(x,y)))] 
# eqs = [Dx(u(x,y)) ~ 0, Dxx(u(x,y)) + Dyy(u(x,y)) ~ (delP/mu) + (K1*Dx(u(x,y)))] 

# BCs 
bcs = [u(x,-H/2) ~ 0, u(x,H/2) ~ 0]

# Spatial domains
domains = [x ∈ Interval(0.0,L), y ∈ Interval(-H/2,H/2)]

# Set up the spatial arrays for solution 
dx = 0.1 
dy = 0.01
xs, ys = [infimum(d.domain):dx:supremum(d.domain) for (dx,d) in zip([dx,dy],domains)]

# Analytical solution
u_a = 1.5*U*(1.0.-((ys.^2)/((H/2)^2)))

# Discretization
dx_grid = dx/2 
dy_grid = dy/2 

# # Neural network
densein = 20 # no. of neurons
act_func = Flux.σ # activation function
chain = [FastChain(FastDense(2,densein,act_func,bias=true,initW=Flux.glorot_normal,initb = Flux.zeros32),
        FastDense(densein,densein,act_func,bias=true,initW=Flux.glorot_normal,initb = Flux.zeros32),
        FastDense(densein,densein,act_func,bias=true,initW=Flux.glorot_normal,initb = Flux.zeros32),
        FastDense(densein,1,bias=true,initW=Flux.glorot_normal,initb = Flux.zeros32)) for _ in 1:length(eqs)] # network of layers
initθ = map(c -> Float64.(c), DiffEqFlux.initial_params.(chain))
flat_initθ = reduce(vcat,initθ)
eltypeθ = eltype(flat_initθ) # determine the type of the elements
parameterless_type_θ = DiffEqBase.parameterless_type(flat_initθ) # generate a parameter-less version of a type
phi = NeuralPDE.get_phi.(chain,parameterless_type_θ) # trial solution of the NN - phi([t,x],θ)
strategy = NeuralPDE.GridTraining([dx_grid,dy_grid])
temp1 = map(phi_ -> phi_(rand(2,10), flat_initθ),phi)

# # Low-level framework begins from here
indvars = [x,y] # define independent variables
depvars = [u(x,y)] #define dependent variables
dim = length(domains)
varbls = NeuralPDE.get_vars(indvars, depvars) # get all the variables in a dictionary
pde_varbls = NeuralPDE.get_variables(eqs,indvars,depvars) # get the variables in the pde
bcs_varbls = NeuralPDE.get_variables(bcs,indvars,depvars) # get the variables in the BCs

derivative = NeuralPDE.get_numeric_derivative() # calculate the derivative (based from Zygote)
integral = NeuralPDE.get_numeric_integral(strategy, indvars, depvars, chain, derivative) # calculate the integral

train_sets = NeuralPDE.generate_training_sets(domains,[dx_grid,dy_grid],eqs,bcs,eltypeθ,indvars,depvars)
pde_train_set,bcs_train_set = train_sets
# build PDE loss function
build_pde_loss_functions = [NeuralPDE.build_loss_function(eq,indvars,depvars,phi,derivative,integral,
                        chain,initθ,strategy,bc_indvars = pde_varbls) for eq in eqs]
temp2 = map(loss_f -> loss_f(rand(2,10), flat_initθ),build_pde_loss_functions)
# Obtain the PDE loss function
get_pde_loss_functions = [NeuralPDE.get_loss_function(_loss_funcs,pde_trainset,eltypeθ,parameterless_type_θ,
                        strategy) for (_loss_funcs,pde_trainset) in zip(build_pde_loss_functions,pde_train_set)]
temp3 = map(l->l(flat_initθ) ,get_pde_loss_functions)

# build BCs loss function
build_bcs_loss_functions = [NeuralPDE.build_loss_function(bc,indvars,depvars,phi,derivative,integral,
                            chain,initθ,strategy,bc_indvars = bcs_vars) 
                            for (bc,bcs_vars) in zip(bcs,bcs_varbls)]
temp4 = map(loss_f -> loss_f(rand(2,10), flat_initθ),build_bcs_loss_functions)
# # Obtain the BCs loss function
get_bc_loss_functions = [NeuralPDE.get_loss_function(_loss_funcs,bcs_trainsets,eltypeθ,parameterless_type_θ,
                        strategy) for (_loss_funcs,bcs_trainsets) in zip(build_bcs_loss_functions,bcs_train_set)]
temp5 = map(l->l(flat_initθ) ,get_bc_loss_functions)

# # Put all the loss functions in a vector array and sum them together
loss_funcs = [get_pde_loss_functions; get_bc_loss_functions]

# # # The total loss function
function loss_func_(θ,p)
    sum(map(l->l(θ) ,loss_funcs))
end

# # # Frame the optimization problem
f = OptimizationFunction(loss_func_, GalacticOptim.AutoZygote())
prob = GalacticOptim.OptimizationProblem(f, flat_initθ)

# # Get the weights of the NN, and values of the total, pde, and BCs loss functions at each iteration
wghts = []
loss_list = []
pde_loss = []
bc_loss = []
cb_ = function (p,l)
       println("Current loss is: $l")
       push!(loss_list, l)
       push!(wghts,p)
       push!(pde_loss,[loss_funcs[i](p) for i in 1:length(eqs)])
       push!(bc_loss,[loss_funcs[i](p) for i in length(eqs)+1:length(eqs)+length(bcs)])
       return false
end

# # # # Choose the optimizer and solve the optimization problem by minimizing the total loss function
res = GalacticOptim.solve(prob, ADAM(0.01); cb = cb_, maxiters=adam_max)
prob = remake(prob,u0=res.minimizer)
res = GalacticOptim.solve(prob, BFGS(); cb = cb_, maxiters=bfgs_max)

acum =  [0;accumulate(+, length.(initθ))]
sep = [acum[i]+1 : acum[i+1] for i in 1:length(acum)-1]
minimizers_ = [res.minimizer[s] for s in sep]
u_p  = [[phi[i]([x,y],minimizers_[i])[1] for x in xs for y in ys] for i in 1:length(eqs)]
u_p1  = reshape(u_p[1], (length(ys),length(xs)))

plot(u_a,ys)
plot!(u_p1[:,5],ys)

[KirillZubov]

you need only one NN because you have only one unknown variable u

[KirillZubov]

so if Dx(u(x,y)) ~ 0 then simplifying Dxx(u(x,y)) + Dyy(u(x,y)) ~ (delP/mu) + ((rho/mu)*u(x,y)*Dx(u(x,y))) => Dxx(u(x,y)) + Dyy(u(x,y)) ~ (delP/mu)

[KirillZubov]

yeah, in this form Dxx(u(x,y)) + Dyy(u(x,y)) ~ (delP/mu) + ((rho/mu)*u(x,y)*Dx(u(x,y))) , optimisation converges to a trivial solution u(x,y) -> 0. A good test for convergence

[KirillZubov]

@v-chau try to use more points in train_data or use QuasiRandomTraining. it should works, I guess