gridap/Gridap.jl

Not able to run the ODE-solvers properly

isakhammer opened this issue · 3 comments

isakhammer commented

Problem statement

Hi Gridap! Thanks for making a absolutely great FEM software library.

I am working on time-dependent problems, but unfortunately am I not able to get your solvers to work properly. Seems like ThetaMethod and RungeKutta_BE_1_0_1 works as expected, but RungeKutta_SDIRK_2_1_2 starts to diverge and RungeKutta_TRBDF2_3_3_2 is crashing. For Newmark and GeneralizedAlpha am I not even able to compile. I am not sure if this is a user error or a bug.

Any idea about this?

Minimalistic Example

I constructed a constructed a minimalistic example of the heat equation with Neumann condtions to test the solvers invidually.

using Dates
using Gridap

function run(;n, u_ex,  dt, solver_choice="ThetaMethod", vtkdirname=nothing)
    # Manufactured solution
    f(t) = x ->  ∂t(u_ex)(x,t) - Δ(u_ex(t))(x)
    ∇u_ex(t) = x ->  ∇(u_ex(t))(x)


    # Setup
    L= 1.0
    domain = (0,L,0,L)
    partition = (n,n)
    model = CartesianDiscreteModel(domain, partition)
    order = 2
    reffe = ReferenceFE(lagrangian, Float64, order)
    V = FESpace( model, reffe, conformity=:H1)
    U = TrialFESpace(V)

    Ω = Triangulation(model)
    Γ = BoundaryTriangulation(model)
    degree = 2*order
    dΩ = Measure(Ω,degree)
    dΓ = Measure(Γ,degree)
    n_Γ = get_normal_vector(Γ)

    g(t) = ∇u_ex(t)⋅n_Γ
    a(u,v) = ∫(∇(v)⋅∇(u))dΩ
    b(v,t) = ∫(v*f(t))dΩ + ∫((g(t)⋅v))dΓ

    res(t,u,v) = a(u,v) + ∫(∂t(u)*v)dΩ - b(v,t)
    jac(t,u,du,v) = a(du,v)
    jac_t(t,u,dut,v) = ∫(dut*v)dΩ

    op = TransientFEOperator(res,jac,jac_t,U,V)

    # Nonlinear/Linear solver
    s =LUSolver()
    # s = NLSolver(LUSolver();show_trace=true,method=:newton) #linesearch=BackTracking())

    function choose_solver(solver_choice, solver_method, dt)
        if solver_choice == "ThetaMethod"
            return ThetaMethod(solver_method, dt, 1)
        elseif solver_choice == "RungeKutta_BE_1_0_1"
            return RungeKutta(solver_method, dt, :BE_1_0_1)
        elseif solver_choice == "RungeKutta_SDIRK_2_1_2"
            return RungeKutta(solver_method, dt, :SDIRK_2_1_2)
        elseif solver_choice == "RungeKutta_TRBDF2_3_3_2"
            return RungeKutta(solver_method, dt, :TRBDF2_3_3_2)
        elseif solver_choice == "Newmark"
            γ, β  = 0.5, 0.25
            return Newmark(solver_method,dt,γ,β)
        elseif solver_choice == "GeneralizedAlpha"
            ρ∞ = 1.0 # Equivalent to Newmark(0.5, 0.25)
            return GeneralizedAlpha(solver_method, dt, ρ∞) # Does not compile
        else
            error("Invalid solver_choice: $solver_choice")
        end
    end

    ode_solver = choose_solver(solver_choice, s, dt)

    # Inital condition
    t_0 = 0.0 # Edit: fixed floating error 
    T = 1.0
    U_0 = interpolate_everywhere(0, U(0.0))

    #################

    U_h_t = solve(ode_solver, op, U_0, t_0, T)

    ts = Float64[]
    el2_ts = Float64[]
    eh1_ts = Float64[]
    eh_energy_ts = Float64[]

    solname = vtkdirname*"/sol_dt_$dt"*"_n_$n"
    mkpath(solname)
    createpvd(solname*".pvd") do pvd
        for (U_h, t) in U_h_t
            e = u_ex(t) - U_h
            pvd[t] = createvtk(Ω, solname*"_$t"*".vtu",cellfields=["u_h"=>U_h,"e"=>e])
            el2_t = sqrt(sum( ∫(e*e)dΩ ))
            eh1_t = sqrt(sum( ∫( e*e + ∇(e)⋅∇(e) )*dΩ ))
            eh_energy = sqrt(sum( ∫(∇(e)⋅∇(e) )*dΩ ))
            push!( ts, t)
            push!( el2_ts, el2_t )
            push!( eh1_ts, eh1_t )
            push!( eh_energy_ts, eh_energy)
        end
    end

    el2s_L2 = sqrt(sum(dt* e_ti^2 for e_ti in el2_ts))
    eh1s_L2 = sqrt(sum(dt* e_ti^2 for e_ti in eh1_ts))
    ehs_energy_L2 = sqrt(sum(dt* e_ti^2 for e_ti in eh_energy_ts))

    el2s_inf = maximum(abs.(el2_ts))
    eh1s_inf = maximum(abs.(eh1_ts))
    ehs_energy_inf = maximum(abs.(eh_energy_ts))

    println("\nTesting solver choice $solver_choice")
    println("L2 Norms:")
    println("el2s_L2: $(round(el2s_L2, digits=4))")
    println("eh1s_L2: $(round(eh1s_L2, digits=4))")
    println("ehs_energy_L2: $(round(ehs_energy_L2, digits=4))")

    println("Infinity Norms:")
    println("el2s_inf: $(round(el2s_inf, digits=4))")
    println("eh1s_inf: $(round(eh1s_inf, digits=4))")
    println("ehs_energy_inf: $(round(ehs_energy_inf, digits=4))")
end


function main()
    dirname= "figures/heat_eq/sim_"*string(Dates.now())
    println(dirname)
    mkpath(dirname)

    u_ex(x,t) = sin(t)*cos(2π*x[1])*sin(2π*x[2])
    u_ex(t::Real) = x -> u_ex(x,t)
    run(n=2^5, dt=2^-4, u_ex=u_ex, vtkdirname=dirname, solver_choice="ThetaMethod")
    run(n=2^5, dt=2^-4, u_ex=u_ex, vtkdirname=dirname, solver_choice="RungeKutta_BE_1_0_1")
    run(n=2^5, dt=2^-4, u_ex=u_ex, vtkdirname=dirname, solver_choice="RungeKutta_SDIRK_2_1_2")
    #run(n=2^5, dt=2^-4, u_ex=u_ex, vtkdirname=dirname, solver_choice="RungeKutta_TRBDF2_3_3_2")
    #run(n=2^5, dt=2^-4, u_ex=u_ex, vtkdirname=dirname, solver_choice="Newmark")
    #run(n=2^5, dt=2^-4, u_ex=u_ex, vtkdirname=dirname, solver_choice="GeneralizedAlpha")
end

main()

Where I get the following output

Testing solver choice ThetaMethod
L2 Norms:
el2s_L2: 0.0001
eh1s_L2: 0.0036
ehs_energy_L2: 0.0036
Infinity Norms:
el2s_inf: 0.0002
eh1s_inf: 0.0056
ehs_energy_inf: 0.0056

Testing solver choice RungeKutta_BE_1_0_1
L2 Norms:
el2s_L2: 0.0001
eh1s_L2: 0.0036
ehs_energy_L2: 0.0036
Infinity Norms:
el2s_inf: 0.0002
eh1s_inf: 0.0056
ehs_energy_inf: 0.0056

Testing solver choice RungeKutta_SDIRK_2_1_2
L2 Norms:
el2s_L2: 3.591747234711817e40
eh1s_L2: 1.2390496633141898e41
ehs_energy_L2: 1.1858488884248314e41
Infinity Norms:
el2s_inf: 1.4366961322166212e41
eh1s_inf: 4.95618880074577e41
ehs_energy_inf: 4.7433860956400755e41

As you can see ThetaMethod and RungeKutta_BE_1_0_1 works properly and RungeKutta_SDIRK_2_1_2 is diverging.

Other Error Messages

If I run RungeKutta_TRBDF2_3_3_2 I get the following crash report.

ERROR: LoadError: MethodError: no method matching update!(::Gridap.ODEs.ODETools.RungeKuttaNonlinearOperator, ::Int64, ::Vector{Vector{Float64}}, ::Int64)
Closest candidates are:
  update!(::Gridap.ODEs.ODETools.RungeKuttaNonlinearOperator, ::Float64, ::AbstractVector, ::Int64) at ~/.julia/packages/Gridap/971dU/src/ODEs/ODETools/RungeKutta.jl:220
Stacktrace:
  [1] solve_step!(uf::Vector{Float64}, solver::RungeKutta, op::Gridap.ODEs.TransientFETools.ODEOpFromFEOp{Gridap.ODEs.ODETools.Nonlinear}, u0::Vector{Float64}, t0::Int64, cache::Nothing)
    @ Gridap.ODEs.ODETools ~/.julia/packages/Gridap/971dU/src/ODEs/ODETools/RungeKutta.jl:114
  [2] solve_step!(uF::Vector{Float64}, solver::RungeKutta, op::Gridap.ODEs.TransientFETools.ODEOpFromFEOp{Gridap.ODEs.ODETools.Nonlinear}, u0::Vector{Float64}, t0::Int64)
    @ Gridap.ODEs.ODETools ~/.julia/packages/Gridap/971dU/src/ODEs/ODETools/ODESolvers.jl:27
  [3] iterate(sol::Gridap.ODEs.ODETools.GenericODESolution{Vector{Float64}})
    @ Gridap.ODEs.ODETools ~/.julia/packages/Gridap/971dU/src/ODEs/ODETools/ODESolutions.jl:47
  [4] iterate(sol::Gridap.ODEs.TransientFETools.TransientFESolution)
    @ Gridap.ODEs.TransientFETools ~/.julia/packages/Gridap/971dU/src/ODEs/TransientFETools/TransientFESolutions.jl:65
  [5] (::var"#131#146"{var"#u_ex#151", String, Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Gridap.ODEs.TransientFETools.TransientFESolution, Gridap.CellData.GenericMeasure, Gridap.Geometry.BodyFittedTriangulation{2, 2, CartesianDiscreteModel{2, Float64, typeof(identity)}, CartesianGrid{2, Float64, typeof(identity)}, Gridap.Arrays.IdentityVector{Int64}}})(pvd::WriteVTK.CollectionFile)
    @ Main ~/code/heat_eq_minimalistic.jl:80
  [6] paraview_collection(f::var"#131#146"{var"#u_ex#151", String, Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Gridap.ODEs.TransientFETools.TransientFESolution, Gridap.CellData.GenericMeasure, Gridap.Geometry.BodyFittedTriangulation{2, 2, CartesianDiscreteModel{2, Float64, typeof(identity)}, CartesianGrid{2, Float64, typeof(identity)}, Gridap.Arrays.IdentityVector{Int64}}}, args::String; kwargs::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
    @ WriteVTK ~/.julia/packages/WriteVTK/UoUwH/src/WriteVTK.jl:168
  [7] paraview_collection(f::Function, args::String)
    @ WriteVTK ~/.julia/packages/WriteVTK/UoUwH/src/WriteVTK.jl:164
  [8] #createpvd#26
    @ ~/.julia/packages/Gridap/971dU/src/Visualization/Vtk.jl:23 [inlined]
  [9] createpvd
    @ ~/.julia/packages/Gridap/971dU/src/Visualization/Vtk.jl:22 [inlined]
 [10] run(; n::Int64, u_ex::var"#u_ex#151", dt::Float64, solver_choice::String, vtkdirname::String)
    @ Main ~/code/heat_eq_minimalistic.jl:79
 [11] main()
    @ Main ~/code/heat_eq_minimalistic.jl:124
 [12] top-level scope
    @ ~/code/heat_eq_minimalistic.jl:130
 [13] include(fname::String)
    @ Base.MainInclude ./client.jl:476
 [14] top-level scope
    @ REPL[1]:1
in expression starting at /root/code/heat_eq_minimalistic.jl:130

And for both Newman and GeneralizedAlpha I simply get this error.

ERROR: LoadError: This function belongs to an interface definition and cannot be used.
ericneiva commented

Hi, @isakhammer,

ERROR: LoadError: MethodError: no method matching update!(::Gridap.ODEs.ODETools.RungeKuttaNonlinearOperator, ::Int64, ::Vector{Vector{Float64}}, ::Int64)
Closest candidates are:
update!(::Gridap.ODEs.ODETools.RungeKuttaNonlinearOperator, ::Float64, ::AbstractVector, ::Int64) at ~/.julia/packages/Gridap/971dU/src/ODEs/ODETools/RungeKutta.jl:220

Concerning the RungeKutta_TRBDF2_3_3_2 error, have you tried setting the initial time in your driver as t_0 = 0.0 instead of t_0 = 0?

isakhammer commented

Good, I changed to t_0 = 0.0 and was now able to now run RungeKutta_TRBDF2_3_3_2. Unfortunately it is diverging, just like RungeKutta_SDIRK_2_1_2.

Testing solver choice RungeKutta_BE_1_0_1
L2 Norms:
el2s_L2: 0.0001
eh1s_L2: 0.0036
ehs_energy_L2: 0.0036
Infinity Norms:
el2s_inf: 0.0002
eh1s_inf: 0.0056
ehs_energy_inf: 0.0056

Testing solver choice RungeKutta_SDIRK_2_1_2
L2 Norms:
el2s_L2: 3.591747234711834e40
eh1s_L2: 1.2390496633141958e41
ehs_energy_L2: 1.185848888424837e41
Infinity Norms:
el2s_inf: 1.4366961322166274e41
eh1s_inf: 4.956188800745793e41
ehs_energy_inf: 4.743386095640098e41

Testing solver choice RungeKutta_TRBDF2_3_3_2
L2 Norms:
el2s_L2: 2.224420149260568e52
eh1s_L2: 1.599167000005805e53
ehs_energy_L2: 1.583620738656578e53
Infinity Norms:
el2s_inf: 8.897679971017945e52
eh1s_inf: 6.39666755076324e53
ehs_energy_inf: 6.334482509749388e53
oriolcg commented

Hi @isakhammer,

Thanks for reporting these errors. The Newmark and GeneralizedAlpha require a tuple with the solution, first time derivative and second time derivative as initial value (displacement, velocity and acceleration), something like this:

U_h_t = solve(ode_solver, op, (U_0,V_0,A_0), t_0, T)

Note that the Newmark only works for 2nd order ODEs, so your operator should also incorporate jac_tt: op = TransientFEOperator(res,jac,jac_t,jac_tt,U,V) (or op = TransientFEOperator(res,U,V,order=2) if using automatic differentiation)

The Runge-Kutta methods have not been tested and there might have errors in the implementation. We would be very grateful if you want to investigate in more detail what is the cause of the problem.

  • 👍1