Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Add Deep ritz to NeuralPDE.jl #857

Draft
wants to merge 5 commits into
base: master
Choose a base branch
from
Draft

Conversation

ayushinav
Copy link
Contributor

Checklist

  • Appropriate tests were added
  • Any code changes were done in a way that does not break public API
  • All documentation related to code changes were updated
  • The new code follows the
    contributor guidelines, in particular the SciML Style Guide and
    COLPRAC.
  • Any new documentation only uses public API

Additional context

Add any other context about the problem here.

@ayushinav
Copy link
Contributor Author

For an eqn as $\Delta u = f(u,x,y)$ with appropriate BCs, we'd want to minimize $\int_{\Omega} \frac{1}{2} |\nabla u|^2 - f(u,x,y) u(x,y) d\Omega$ for Deep-Ritz algorithm. I think my implementation minimizes the square of the integrand above. I guess the way to proceed would be to dispatch on SciMLBase.symbolic_discretize, though most of the code would essentially remain the same until the loss function. Just wanted to make sure this would be a good way to proceed and know if there are any other ideas.

@ayushinav
Copy link
Contributor Author

# using Pkg

# Pkg.activate("./NeuralPDE.jl/.")

using NeuralPDE

using ModelingToolkit, Optimization, OptimizationOptimisers, Distributions, MethodOfLines,
      OrdinaryDiffEq
import ModelingToolkit: Interval, infimum, supremum
using Lux #: tanh, identity

# To make a similar architecture as theirs
function block(inp, out, hid, n, act)
    chain =  Chain(Dense(inp, hid, act), [Dense(hid, hid, act) for _ in 2:n-1]..., Dense(hid, out, act))
    SkipConnection(chain, +)
end

# the authors use the commented activation fn. Here I'm using tanh.
drm_act(x) = tanh(x); #max(x^3, 0)

@parameters x y
@variables u(..)
Dx= Differential(x)
Dy= Differential(y)
Dxx = Differential(x)^2
Dyy = Differential(y)^2

# 2D PDE
#=
The first commented out eq is the one we should be getting. 
The second is modified form of the first one so that we get the right integrand, and not the square of it.
=#
# eq = Dxx(u(x, y)) + Dyy(u(x, y)) ~ -sin(pi * x) * sin(pi * y)
# eq = sqrt(abs(0.5 * (Dx(u(x,y)))^2 + 0.5 * (Dy(u(x,y)))^2 + sin(pi * x) * sin(pi * y) * u(x,y))) ~ 0

# The current implementation returns the following PDE 
(except that I have a bug and do not have the 0.5 multiplied to the first term.
eq = 0.5 * (Dx(u(x,y)))^2 + 0.5 * (Dy(u(x,y)))^2 + sin(pi * x) * sin(pi * y) * u(x,y) ~ 0

# Initial and boundary conditions
bcs = [u(0, y) ~ 0.0, u(1, y) ~ -sin(pi * 1) * sin(pi * y),
    u(x, 0) ~ 0.0, u(x, 1) ~ -sin(pi * x) * sin(pi * 1)]

# Space and time domains
domains = [x ∈ Interval(0.0, 1.0), y ∈ Interval(0.0, 1.0)]
nonadaptive_loss = NeuralPDE.NonAdaptiveLoss(pde_loss_weights = 10, bc_loss_weights = 1000)

strategy = QuasiRandomTraining(512, minibatch = 32)
chain_ = Chain(Dense(2,10), # similar architecture as the one in the paper 
    block(10, 10, 10, 2, drm_act),
    block(10, 10, 10, 2, drm_act),
    block(10, 10, 10, 2, drm_act),
    block(10, 10, 10, 2, drm_act),
    Dense(10,1))

discretization = PhysicsInformedNN(chain_, strategy)

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

global iter = 0
callback = function (p, l)
    global iter += 1
    if iter % 50 == 0
        println("$iter => $l")
    end
    return false
end

res = Optimization.solve(prob, Adam(1e-3); callback = callback, maxiters = 500)
prob = remake(prob, u0 = res.u)
res = Optimization.solve(prob, Adam(1e-5); callback = callback, maxiters = 500)
prob = remake(prob, u0 = res.u)
res = Optimization.solve(prob, Adam(1e-7); callback = callback, maxiters = 1000)
phi = discretization.phi

xs, ys = [infimum(d.domain):0.01:supremum(d.domain) for d in domains]
analytic_sol_func(x, y) = (sin(pi * x) * sin(pi * y)) / (2pi^2)

u_predict = reshape([first(phi([x, y], res.u)) for x in xs for y in ys],
    (length(xs), length(ys)));
u_real = reshape([analytic_sol_func(x, y) for x in xs for y in ys],
    (length(xs), length(ys)));

    
    # @test u_predict≈u_real atol=0.1
# end

diff_u = abs.(u_real .- u_predict);
using Plots
p1 = plot(xs, ys, u_real, linetype = :contourf, title = "analytic");
p2 = plot(xs, ys, u_predict, linetype = :contourf, title = "predict");
p3 = plot(xs, ys, diff_u, linetype = :contourf, title = "error");
plot(p1, p2, p3)

Since the above code does not give good results, I think the issue if with the quadrature but not completely sure because even the second eqn which is commented out does not give good results.

@sathvikbhagavan
Copy link
Member

The loss function generation can be dispatched but it might require refactoring. Possibly, #678 might also be needed.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

2 participants