Hi everyone,
I am trying to solve a PDE with the julia package Neural.jl.
I have 3 bcs. At two boundarys the pressure has to be zero. The last bsc ist a periodic boundary condition. The pressure at the one side has to be the same as the pressure at the other side.
My final loss is (depending on the training strategie) around 1e-07 for grid-training and 1e-11 for quadraturTraining.
But the periodic bcs is not fullfullid at all. Differences around (0.5 for quadraturTraining and 0.04 for gridtraining.
The code I am using:
using NeuralPDE, Lux, Optimization, OptimizationOptimJL, Plots
import ModelingToolkit: Interval
function main()
@parameters x y
@variables p(..)
Dx = Differential(x)
Dy = Differential(y)
Dxx = Differential(x)^2
# Bearing
ω = 600.0 # [rad/s]
D = 18e-03 # [m]
B = 10e-03
psi = 3.58e-03
eta = 6.4e-03
R = D/2
C = psi * D/2
u = sign(ω*D/2);
p_fak = eta * abs(ω*R) * R/(C^2)
# 2D PDE+
eta_nd = 1
x0 = -1e-05
y0 = 1e-05
phi = atan(x0,y0)
e = sqrt(x0^2 + y0^2)
H_spalt(x,y) = (C + e*cos(pi/2 + phi + x))/C
eq = [Dx(H_spalt(x,y)^3/12 * Dx(p(x,y))) + Dy(H_spalt(x,y)^3/12 * Dy(p(x,y))) ~ Dx(H_spalt(x,y))]
# Boundary conditions
bcs = [p(x,0) ~ 0, # pressure at the left and rigth side
p(x,B/R) ~ 0,
p(0,y) ~ p(2*pi,y)] # periodic boundary condition
# Space and time domains
domains = [x ∈ Interval(0.0,2*pi),
y ∈ Interval(0.0,B/R)]
# Neural network
dim = 2 # number of dimensions
n_neurons = 16
af = Lux.tanh
chain = Lux.Chain(Dense(dim,n_neurons,af),
Dense(n_neurons,n_neurons,af),
Dense(n_neurons,1))
# Discretization
nx = 60
ny = 40
dx = pi*D/nx/R
dy = B/(ny-1)/R
strategy = GridTraining([dx,dy])
strategy = QuadratureTraining(batch = 2000)
discretization = PhysicsInformedNN(chain,strategy,
adaptiv_loss = GradientScaleAdaptiveLoss(1))
@named pde_system = PDESystem(eq,bcs,domains,[x,y],[p(x, y)])
prob = discretize(pde_system,discretization)
#Optimizer
opt = OptimizationOptimJL.BFGS()
#Callback function
callback = function (p,l)
println("Current loss is: $l")
return false
end
@time res =
Optimization.solve(prob, opt, callback = callback, maxiters=1000)
phi = discretization.phi
return domains, res, phi, dx, dy, p_fak
end
function post_process(domains,res,phi,dx ,dy, p_fak)
# Bearing
D = 18e-03 # [m]
B = 10e-03
R = D/2
xs = collect(0:dx:2*pi)
ys = collect(0:dy:B/R)
u_predict = zeros(length(xs),length(ys))
for i = eachindex(xs)
for j = eachindex(ys)
u_predict[i,j] = first(phi([xs[i],ys[j]],res.u))*1
if u_predict[i,j] < 0
u_predict[i,j] = 0
end
end
end
plotly()
p = plot(ys,xs,u_predict,seriestype=:surface, camera=(70, 30),size = (800, 600))
display(p)
return u_predict
end
domains, res, phi, dx, dy, p_fak = main()
u_predict = post_process(domains,res,phi, dx, dy, p_fak)
And the result:
Thanks for your help and a happy new year:)