Merge pull request #243 from SciML/update_quadrature_strategy

Update quadrature strategy

Author: ChrisRackauckas

src/pinns_pde_solve.jl

@@ -601,13 +601,13 @@ function get_loss_function(loss_functions, bounds, strategy::QuadratureTraining;
     end
 
     f = (lb,ub,loss_,θ) -> begin
-        _loss = (x,θ) -> sum(inner_loss(loss_, x, θ))
+        _loss = (x,θ) -> inner_loss(loss_, x, θ)
         prob = QuadratureProblem(_loss,lb,ub,θ;batch = strategy.batch)
-        solve(prob,
+        abs(solve(prob,
               strategy.quadrature_alg,
               reltol = strategy.reltol,
               abstol = strategy.abstol,
-              maxiters = strategy.maxiters)[1]
+              maxiters = strategy.maxiters)[1])
     end
     loss = (θ) -> τ*sum(f(lb,ub,loss_,θ) for (lb,ub,loss_) in zip(lbs,ubs,loss_functions))
     return loss
@@ -723,11 +723,6 @@ function DiffEqBase.discretize(pde_system::PDESystem, discretization::PhysicsInf
                                                        strategy)
          (pde_loss_function, bc_loss_function)
     elseif strategy isa QuadratureTraining
-        dim<=1 && error("QuadratureTraining works only with dimensionality more than 1")
-
-        (dim<=2 && strategy.quadrature_alg in [CubaCuhre(),CubaDivonne()]
-        && error("$(strategy.quadrature_alg) works only with dimensionality more than 2"))
-
         bounds = get_bounds(domains,bcs,dict_indvars,dict_depvars)
         pde_bounds, bcs_bounds = bounds
 
@@ -745,12 +740,39 @@ function DiffEqBase.discretize(pde_system::PDESystem, discretization::PhysicsInf
                                               strategy;
                                               τ=τp)
 
-        τb =  1.0f0 / (bsl*Int(round(τ_^(bl/pl))))
-
-        bc_loss_function = get_loss_function(_bc_loss_functions,
+        τb =  1.0f0 / (bsl * τ_^(bl/pl))
+
+        if bl == 0
+            bc_loss_function = get_loss_function(_bc_loss_functions,
+                                                 [[[]]],
+                                                 GridTraining(0.1))
+
+        elseif bl == 1 && strategy.quadrature_alg in [CubaCuhre(),CubaDivonne()]
+            @warn "$(strategy.quadrature_alg) does not work with one-dimensional
+            problems, so for the boundary conditions loss function,
+            the quadrature algorithm was replaced by HCubatureJL"
+
+            strategy = QuadratureTraining(quadrature_alg = HCubatureJL(),
+                                          reltol = bsl*(strategy.reltol)^(bl/pl),
+                                          abstol = bsl*(strategy.abstol)^(bl/pl),
+                                          maxiters = bsl*Int(round((strategy.maxiters)^(bl/pl))),
+                                          batch = strategy.batch)
+            bc_loss_function = get_loss_function(_bc_loss_functions,
+                                                 bcs_bounds,
+                                                 strategy;
+                                                 τ=τb)
+        else
+            strategy = QuadratureTraining(quadrature_alg = strategy.quadrature_alg,
+                                          reltol = bsl*(strategy.reltol)^(bl/pl),
+                                          abstol = bsl*(strategy.abstol)^(bl/pl),
+                                          maxiters = bsl*Int(round((strategy.maxiters)^(bl/pl))),
+                                          batch = strategy.batch)
+            bc_loss_function = get_loss_function(_bc_loss_functions,
                                              bcs_bounds,
                                              strategy;
                                              τ=τb)
+        end
+
         (pde_loss_function, bc_loss_function)
     end
 

test/NNPDE_tests.jl

@@ -37,7 +37,7 @@ dt = 0.1
 # Neural network
 chain = FastChain(FastDense(1,12,Flux.σ),FastDense(12,1))
 
-strategy = NeuralPDE.GridTraining(dt)
+strategy = NeuralPDE.QuadratureTraining()
 discretization = NeuralPDE.PhysicsInformedNN(chain,
                                              strategy;
                                              init_params = nothing,
@@ -160,8 +160,7 @@ function run_2d_poisson_equation(strategy)
     res = GalacticOptim.solve(prob, ADAM(0.01); cb = cb,  maxiters=2)
 end
 
-algs = [CubaVegas(), CubaSUAVE(),HCubatureJL(), CubatureJLh(), CubatureJLp()]
-#CubaDivonne(),CubaCuhre() doesn't work with dim = 2
+algs = [HCubatureJL(), CubatureJLh(), CubatureJLp(),CubaCuhre()]
 for alg in algs
     strategy =  NeuralPDE.QuadratureTraining(quadrature_alg = alg,reltol=1e-8,abstol=1e-8,maxiters=600)
     run_2d_poisson_equation(strategy)
@@ -232,12 +231,12 @@ domains = [x ∈ IntervalDomain(0.0,1.0), y ∈ IntervalDomain(0.0,1.0)]
 # Neural network
 chain1 = FastChain(FastDense(2,10,Flux.σ),FastDense(10,1))
 chain2 = FastChain(FastDense(2,10,Flux.σ),FastDense(10,1))
-
-discretization = NeuralPDE.PhysicsInformedNN([chain1,chain2],NeuralPDE.GridTraining([0.1, 0.1]))
+strategy = NeuralPDE.QuadratureTraining(quadrature_alg=CubaCuhre(),reltol= 1e-4,abstol= 1e-3,maxiters=20)
+discretization = NeuralPDE.PhysicsInformedNN([chain1,chain2],strategy)
 pde_system = PDESystem(eqs,bcs,domains,[x,y],[u1,u2])
 prob = NeuralPDE.discretize(pde_system,discretization)
 
-res = GalacticOptim.solve(prob,Optim.BFGS(); cb = cb, maxiters=600)
+res = GalacticOptim.solve(prob,Optim.BFGS(); cb = cb, maxiters=300)
 phi = discretization.phi
 
 analytic_sol_func(x,y) =[1/3*(6x - y), 1/2*(6x - y)]
@@ -250,7 +249,7 @@ sep = [acum[i]+1 : acum[i+1] for i in 1:length(acum)-1]
 minimizers = [res.minimizer[s] for s in sep]
 u_predict  = [[phi[i]([x,y],minimizers[i])[1] for x in xs  for y in ys] for i in 1:2]
 
-@test u_predict ≈ u_real atol = 10.0
+@test u_predict ≈ u_real atol = 2.0
 
 # p1 =plot(xs, ys, u_predict, st=:surface);
 # p2 = plot(xs, ys, u_real, st=:surface);
@@ -300,19 +299,18 @@ train_sets = NeuralPDE.generate_training_sets(domains,dx,bcs,indvars,depvars)
 pde_train_set,bcs_train_set,train_set = train_sets
 pde_bounds, bcs_bounds = NeuralPDE.get_bounds(domains,bcs,indvars,depvars)
 
-quadrature_strategy = NeuralPDE.QuadratureTraining(quadrature_alg=HCubatureJL(),reltol= 1e-3,abstol= 1e-3,maxiters=20)
-τp = 1/100
+quadrature_strategy = NeuralPDE.QuadratureTraining(quadrature_alg=CubaCuhre(),reltol= 1e-4,abstol= 1e-3,maxiters=20)
 pde_loss_function = NeuralPDE.get_loss_function(_pde_loss_function,
                                                 pde_bounds,
                                                 quadrature_strategy;
-                                                τ = τp)
+                                                τ = 1/100)
 
 
-grid_strategy = NeuralPDE.GridTraining(dx)
+quadrature_strategy = NeuralPDE.QuadratureTraining(quadrature_alg=HCubatureJL(),reltol= 1e-2,abstol= 1e-1,maxiters=5)
 bc_loss_function = NeuralPDE.get_loss_function(_bc_loss_functions,
-                                               bcs_train_set,
-                                               grid_strategy;
-                                               τ = nothing)
+                                               bcs_bounds,
+                                               quadrature_strategy;
+                                               τ = 1/20)
 
 function loss_function_(θ,p)
     return pde_loss_function(θ) + bc_loss_function(θ)
@@ -321,7 +319,12 @@ end
 f_ = OptimizationFunction(loss_function_, GalacticOptim.AutoZygote())
 prob = GalacticOptim.OptimizationProblem(f_, initθ)
 
-res = GalacticOptim.solve(prob,Optim.BFGS(); cb = cb, maxiters=400)
+cb_ = function (p,l)
+    println("Current losses are: ", pde_loss_function(p), " , ",  bc_loss_function(p))
+    return false
+end
+
+res = GalacticOptim.solve(prob,Optim.BFGS(); cb = cb_, maxiters=400)
 
 xs,ts = [domain.domain.lower:dx:domain.domain.upper for domain in domains]
 analytic_sol_func(x,t) =  sum([(8/(k^3*pi^3)) * sin(k*pi*x)*cos(C*k*pi*t) for k in 1:2:50000])