Thank you!
Yes - this looks more like wave.
What did you changed, @Knud_Sorensen ?
I change the network.
inner = 10 chain = Chain(Dense(3, inner, Lux.asinh), Dense(inner, inner, Lux.asinh), Dense(inner, inner, Lux.asinh), Dense(inner, 1))
Have you tried this with MethodOfLines? Second order wave is tested and working in v0.9.7
Thank you, @xtalax for looking at this!
Do you have an example which I can look at? I couldn’t find it at the website.
Thank you!
Indeed - that makes a bit difference - couldn’t get to the really low error, though… Will work around it… Need to setup some problem where I know the solution in advance - to check how it actually converges… Still have lot’s of questions… and no intuition about what is going on…
You can use the same PDESystem, just use the standard discretization that you can find all over the tutorials and examples,
as an example
...
disc = MOLFiniteDifference([x => 30, y => 30], t)
prob = discretize(sys, disc)
sol = solve(prob, Tsit5()) # Try different solvers for best results
udisc = sol[u(t, x, y)] # gives an array in t, x and y respectively
The trouble with NeuralPDE is that you have no guarantee that you will converge to something correct in a sensible amount of time although it is more flexible, with FD i.e. MOL you know at least that it is correct to within a margin of error, unless your discretization is unstable which is fairly obvious to see as it will blow up to infinity.
Also if you are interested in inhomogenous electromagnetics/acoustics, please reach out on the slack
Yes, I reduced the network to be able to run it on my Laptop. That might be why the loss is not so low.
1 Like
Thank you, @xtalax - I’ll follow up with my updated code (once I’ll test it).
In the meanwhile 2 questions:
- Can I turn it into GPU based solver?
- Is there any examples for solving inverse problems (given the wave fields reconstruct/update the velocity volume?)
I’ve installed 0.9.7 using add MethodOfLines#master.
For now, I’m getting this (with 0.9.7 for Brusselator example from the Getting Started tutorial):
Discretization:
┌ Warning: The system contains interface boundaries, which are not compatible with system transformation. The system will not be transformed. Please post an issue if you need this feature.
└ @ MethodOfLines C:\Users\USER\.julia\packages\MethodOfLines\BVnkB\src\system_parsing\pde_system_transformation.jl:39
24.688898 seconds (124.26 M allocations: 8.412 GiB, 8.73% gc time, 36.44% compilation time: 82% of which was recompilation)
Solve:
218.603521 seconds (205.71 M allocations: 10.883 GiB, 1.81% gc time, 47.73% compilation time)
So, I coded it up and it gives an issue:
using ModelingToolkit, MethodOfLines, OrdinaryDiffEq, DomainSets
println("Setting Parameters")
@parameters t x y
@variables u(..)
Dx = Differential(x)
Dy = Differential(y)
Dxx = Differential(x)^2
Dyy = Differential(y)^2
Dt = Differential(t)
Dtt = Differential(t)^2
t_min = 0.0
t_max = 2.0
x_min = -2.0
x_max = 2.0
x_center = (x_min + x_max)/2
y_min = -2.0
y_max = 2.0
y_center = (y_min + y_max)/2
c = 1.0
dx=0.05
dy=0.05
nx = 5
ny = 5
discOrder=2
first_maxiter = 250
println("Setting PDE")
# 2D PDE
eq = Dtt(u(t, x, y)) ~ c*Dx(c*Dx(u(t, x, y))) + c*Dy(c*Dy(u(t, x, y)))
println("Setting IC and BC")
open_bc = [
bc_scale*Dt(u(t,x_min,y)) ~ bc_scale*c*Dx(u(t,x_min,y)),
bc_scale*Dt(u(t,x_max,y)) ~ -bc_scale*c*Dx(u(t,x_max,y)),
bc_scale*Dt(u(t,x,y_min)) ~ bc_scale*c*Dx(u(t,x,y_min)),
bc_scale*Dt(u(t,x,y_max)) ~ -bc_scale*c*Dx(u(t,x,y_max))
]
# impulse source IC
ic_func(x,y) = (x-x_center+0.2)*exp(-12*((x-x_center)^2+(y-y_center)^2))
ic = [
u(t_min,x,y) ~ ic_func(x, y),
Dt(u(t_min,x,y)) ~ 0.0
]
all_conditions = append!(open_bc, ic)
println("Setting Space and Domain")
# Space and time domains
domains = [t ∈ Interval(t_min, t_max),
x ∈ Interval(x_min, x_max),
y ∈ Interval(y_min, y_max)]
println("Putting it all together into PDESystem")
@named pde_system = PDESystem(eq, all_conditions, domains, [t, x, y], [u(t, x, y)])
println("Discretization schema definitation:")
discretization = MOLFiniteDifference(
[x=>nx, y=>ny],
t,
approx_order=discOrder
)
println("Putting together PDESystem and Discretization schema:")
@time problem = discretize(pde_system, discretization)
println("Start to solve using TRBDF2 with saving at 0.1:")
callback = function (params, loss, other_args)
println("Current loss is: $loss")
return false
end
# using Tsit5() - would be another thing to try
@time res = solve(
problem,
TRBDF2();
callback = callback, maxiters = first_maxiter
)
Which gives me this:
Summary
The system of equations is:
Equation[4.0(u(t))[2, 2] + Differential(t)(Differential(t)((u(t))[2, 2])) - (u(t))[1, 2] - (u(t))[2, 1] - (u(t))[2, 3] - (u(t))[3, 2] ~ 0, 4.0(u(t))[3, 2] + Differential(t)(Differential(t)((u(t))[3, 2])) - (u(t))[2, 2] - (u(t))[3, 1] - (u(t))[3, 3] - (u(t))[4, 2] ~ 0, 4.0(u(t))[4, 2] + Differential(t)(Differential(t)((u(t))[4, 2])) - (u(t))[3, 2] - (u(t))[4, 1] - (u(t))[4, 3] - (u(t))[5, 2] ~ 0, 4.0(u(t))[2, 3] + Differential(t)(Differential(t)((u(t))[2, 3])) - (u(t))[1, 3] - (u(t))[2, 2] - (u(t))[2, 4] - (u(t))[3, 3] ~ 0, 4.0(u(t))[3, 3] + Differential(t)(Differential(t)((u(t))[3, 3])) - (u(t))[2, 3] - (u(t))[3, 2] - (u(t))[3, 4] - (u(t))[4, 3] ~ 0, 4.0(u(t))[4, 3] + Differential(t)(Differential(t)((u(t))[4, 3])) - (u(t))[3, 3] - (u(t))[4, 2] - (u(t))[4, 4] - (u(t))[5, 3] ~ 0, 4.0(u(t))[2, 4] + Differential(t)(Differential(t)((u(t))[2, 4])) - (u(t))[1, 4] - (u(t))[2, 3] - (u(t))[2, 5] - (u(t))[3, 4] ~ 0, 4.0(u(t))[3, 4] + Differential(t)(Differential(t)((u(t))[3, 4])) - (u(t))[2, 4] - (u(t))[3, 3] - (u(t))[3, 5] - (u(t))[4, 4] ~ 0, 4.0(u(t))[4, 4] + Differential(t)(Differential(t)((u(t))[4, 4])) - (u(t))[3, 4] - (u(t))[4, 3] - (u(t))[4, 5] - (u(t))[5, 4] ~ 0, Differential(t)((u(t))[1, 2]) ~ 2.0(u(t))[2, 2] - 0.5(u(t))[3, 2] - 1.5(u(t))[1, 2], Differential(t)((u(t))[1, 3]) ~ 2.0(u(t))[2, 3] - 0.5(u(t))[3, 3] - 1.5(u(t))[1, 3], Differential(t)((u(t))[1, 4]) ~ 2.0(u(t))[2, 4] - 0.5(u(t))[3, 4] - 1.5(u(t))[1, 4], Differential(t)((u(t))[5, 2]) ~ 2.0(u(t))[4, 2] - 0.5(u(t))[3, 2] - 1.5(u(t))[5, 2], Differential(t)((u(t))[5, 3]) ~ 2.0(u(t))[4, 3] - 0.5(u(t))[3, 3] - 1.5(u(t))[5, 3], Differential(t)((u(t))[5, 4]) ~ 2.0(u(t))[4, 4] - 0.5(u(t))[3, 4] - 1.5(u(t))[5, 4], Differential(t)((u(t))[2, 1]) ~ Differential(-1.0)((u(t))[2, 1]), Differential(t)((u(t))[3, 1]) ~ Differential(0.0)((u(t))[3, 1]), Differential(t)((u(t))[4, 1]) ~ Differential(1.0)((u(t))[4, 1]), Differential(t)((u(t))[2, 5]) ~ -Differential(-1.0)((u(t))[2, 5]), Differential(t)((u(t))[3, 5]) ~ -Differential(0.0)((u(t))[3, 5]), Differential(t)((u(t))[4, 5]) ~ -Differential(1.0)((u(t))[4, 5]), (u(t))[1, 1] ~ 0, (u(t))[5, 1] ~ 0, (u(t))[1, 5] ~ 0, (u(t))[5, 5] ~ 0]
Discretization failed, please post an issue on https://github.com/SciML/MethodOfLines.jl with the failing code and system at low point count.
ERROR: ArgumentError: Differential w.r.t. multiple variables Any[0.0, -1.0, t, 1.0] are not allowed.
Stacktrace:
[1] check_equations(eqs::Vector{Equation}, iv::SymbolicUtils.BasicSymbolic{Real})
@ ModelingToolkit C:\Users\USER\.julia\packages\ModelingToolkit\ua9Sp\src\utils.jl:202
[2] ODESystem(tag::UInt64, deqs::Vector{Equation}, iv::SymbolicUtils.BasicSymbolic{Real}, dvs::Vector{SymbolicUtils.BasicSymbolic{Real}}, ps::Vector{Any}, tspan::Nothing, var_to_name::Dict{Any, Any}, ctrls::Vector{Any}, observed::Vector{Equation}, tgrad::Base.RefValue{Vector{Num}}, jac::Base.RefValue{Any}, ctrl_jac::Base.RefValue{Any}, Wfact::Base.RefValue{Matrix{Num}}, Wfact_t::Base.RefValue{Matrix{Num}}, name::Symbol, systems::Vector{ODESystem}, defaults::Dict{Any, Any}, torn_matching::Nothing, connector_type::Nothing, preface::Nothing, cevents::Vector{ModelingToolkit.SymbolicContinuousCallback}, devents::Vector{ModelingToolkit.SymbolicDiscreteCallback}, metadata::MethodOfLines.MOLMetadata{Val{true}(), MethodOfLines.DiscreteSpace{2, 1, MethodOfLines.CenterAlignedGrid}, MOLFiniteDifference{MethodOfLines.CenterAlignedGrid, MethodOfLines.ScalarizedDiscretization}, PDESystem, Base.RefValue{Any}, MethodOfLines.ScalarizedDiscretization}, gui_metadata::Nothing, tearing_state::Nothing, substitutions::Nothing, complete::Bool, discrete_subsystems::Nothing, unknown_states::Nothing; checks::Bool)
@ ModelingToolkit C:\Users\USER\.julia\packages\ModelingToolkit\ua9Sp\src\systems\diffeqs\odesystem.jl:154
[3] ODESystem
@ C:\Users\USER\.julia\packages\ModelingToolkit\ua9Sp\src\systems\diffeqs\odesystem.jl:143 [inlined]
[4] ODESystem(deqs::Vector{Equation}, iv::Num, dvs::Vector{SymbolicUtils.BasicSymbolic{Real}}, ps::Vector{Num}; controls::Vector{Num}, observed::Vector{Equation}, systems::Vector{ODESystem}, tspan::Nothing, name::Symbol, default_u0::Dict{Any, Any}, default_p::Dict{Any, Any}, defaults::Dict{SymbolicUtils.BasicSymbolic{Real}, Float64}, connector_type::Nothing, preface::Nothing, continuous_events::Nothing, discrete_events::Nothing, checks::Bool, metadata::MethodOfLines.MOLMetadata{Val{true}(), MethodOfLines.DiscreteSpace{2, 1, MethodOfLines.CenterAlignedGrid}, MOLFiniteDifference{MethodOfLines.CenterAlignedGrid, MethodOfLines.ScalarizedDiscretization}, PDESystem, Base.RefValue{Any}, MethodOfLines.ScalarizedDiscretization}, gui_metadata::Nothing)
@ ModelingToolkit C:\Users\USER\.julia\packages\ModelingToolkit\ua9Sp\src\systems\diffeqs\odesystem.jl:218
[5] generate_system(disc_state::PDEBase.EquationState, s::MethodOfLines.DiscreteSpace{2, 1, MethodOfLines.CenterAlignedGrid}, u0::Vector{Pair{SymbolicUtils.BasicSymbolic{Real}, Float64}}, tspan::Tuple{Float64, Float64}, metadata::MethodOfLines.MOLMetadata{Val{true}(), MethodOfLines.DiscreteSpace{2, 1, MethodOfLines.CenterAlignedGrid}, MOLFiniteDifference{MethodOfLines.CenterAlignedGrid, MethodOfLines.ScalarizedDiscretization}, PDESystem, Base.RefValue{Any}, MethodOfLines.ScalarizedDiscretization}, disc::MOLFiniteDifference{MethodOfLines.CenterAlignedGrid, MethodOfLines.ScalarizedDiscretization})
@ PDEBase C:\Users\USER\.julia\packages\PDEBase\2OwkC\src\discretization_state.jl:41
[6] symbolic_discretize(pdesys::PDESystem, discretization::MOLFiniteDifference{MethodOfLines.CenterAlignedGrid, MethodOfLines.ScalarizedDiscretization})
@ PDEBase C:\Users\USER\.julia\packages\PDEBase\2OwkC\src\symbolic_discretize.jl:88
[7] discretize(pdesys::PDESystem, discretization::MOLFiniteDifference{MethodOfLines.CenterAlignedGrid, MethodOfLines.ScalarizedDiscretization}; analytic::Nothing,
kwargs::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
@ PDEBase C:\Users\USER\.julia\packages\PDEBase\2OwkC\src\discretization_state.jl:58
[8] discretize(pdesys::PDESystem, discretization::MOLFiniteDifference{MethodOfLines.CenterAlignedGrid, MethodOfLines.ScalarizedDiscretization})
@ PDEBase C:\Users\USER\.julia\packages\PDEBase\2OwkC\src\discretization_state.jl:55
[9] top-level scope
@ .\timing.jl:273
Any thoughts?
It’s not possible to discretize derivatives that aren’t orthogonal to the boundary in bcs with MOL, are you sure that the spatial derivs in the following:
bc_scale*Dt(u(t,x,y_min)) ~ bc_scale*c*Dx(u(t,x,y_min)),
bc_scale*Dt(u(t,x,y_max)) ~ -bc_scale*c*Dx(u(t,x,y_max))
shouldn’t be in the y direction?
Yes, @xtalax , you are absolutely right - another typo on my behalf. Indeed - this should be Dy, not Dx.
Here is updated code, which works for any nx=ny<=21 - but fails with the error below once I increase it to something bigger.
using ModelingToolkit, MethodOfLines, OrdinaryDiffEq, DomainSets
println("Setting Parameters")
@parameters t x y
@variables u(..)
Dx = Differential(x)
Dy = Differential(y)
Dxx = Differential(x)^2
Dyy = Differential(y)^2
Dt = Differential(t)
Dtt = Differential(t)^2
t_min = 0.0
t_max = 5.0
x_min = -2.0
x_max = 2.0
x_center = (x_min + x_max)/2
y_min = -2.0
y_max = 2.0
y_center = (y_min + y_max)/2
c = 1.0
dx=0.05
dy=0.05
nx = 21
ny = 21
discOrder=2
first_maxiter = 250
eq_scale=1.000
bc_scale=1.000
ic_scale=1.000
println("Setting PDE")
# 2D PDE
eq = eq_scale*Dtt(u(t, x, y)) ~ eq_scale*(c*Dx(c*Dx(u(t, x, y))) + c*Dy(c*Dy(u(t, x, y))))
println("Setting IC and BC")
open_bc = [
bc_scale*Dt(u(t,x_min,y)) ~ bc_scale*c*Dx(u(t,x_min,y)),
bc_scale*Dt(u(t,x_max,y)) ~ -bc_scale*c*Dx(u(t,x_max,y)),
bc_scale*Dt(u(t,x,y_min)) ~ bc_scale*c*Dy(u(t,x,y_min)),
bc_scale*Dt(u(t,x,y_max)) ~ -bc_scale*c*Dy(u(t,x,y_max))
]
# impulse source IC
ic_func(x,y) = (x-x_center+0.2)*exp(-12*((x-x_center)^2+(y-y_center)^2))
ic = [
ic_scale*u(t_min,x,y) ~ ic_scale*ic_func(x, y),
ic_scale*Dt(u(t_min,x,y)) ~ 0.0
]
all_conditions = append!(open_bc, ic)
println("Setting Space and Domain")
# Space and time domains
domains = [t ∈ Interval(t_min, t_max),
x ∈ Interval(x_min, x_max),
y ∈ Interval(y_min, y_max)]
println("Putting it all together into PDESystem")
@named pde_system = PDESystem(eq, all_conditions, domains, [t, x, y], [u(t, x, y)])
println("Discretization schema definitation:")
discretization = MOLFiniteDifference(
[x=>nx, y=>ny],
t,
approx_order=discOrder
)
println("Putting together PDESystem and Discretization schema:")
@time problem = discretize(pde_system, discretization)
println("Start to solve using TRBDF2 with saving at 0.1:")
# using Tsit5() - would be another thing to try
@time res = solve(
problem,
TRBDF2();
# maxiters = first_maxiter
)
println("Inspect the solution")
ts, xs, ys = [infimum(d.domain):0.1:supremum(d.domain) for d in domains]
using Plots
using Printf
function plot_heatmap(res)
discrete_x = res[x]
discrete_y = res[y]
discrete_t = res[t]
solu = res[u(t, x, y)]
anim = @animate for k in 1:length(discrete_t)
heatmap(solu[k, 2:end, 2:end], title="$(discrete_t[k])") # 2:end since end = 1, periodic condition
end
gif(anim, "heatmap_wave.gif", fps = 10)
end
println("Third plot")
plot_heatmap(res)
Error:
ERROR: ArgumentError: Equations (880), states (872), and initial conditions (872) are of different lengths.
To allow a different number of equations than states use kwarg check_length=false.
Full stack error
Setting Parameters
Setting PDE
Setting IC and BC
Setting Space and Domain
Putting it all together into PDESystem
Discretization schema definitation:
Putting together PDESystem and Discretization schema:
The system of equations is:
>>>> I REMOVED IT SINCE IT WAS TOO LONG TO PASTE HERE <<<<
[1:484]
ERROR: ArgumentError: Equations (880), states (872), and initial conditions (872) are of different lengths. To allow a different number of equations than states use kwarg check_length=false.
Stacktrace:
[1] check_eqs_u0(eqs::Vector{Equation}, dvs::Vector{Any}, u0::Vector{Float64}; check_length::Bool, kwargs::Base.Pairs{Symbol, Real, Tuple{Symbol, Symbol}, NamedTuple{(:t, :has_difference), Tuple{Float64, Bool}}})
@ ModelingToolkit C:\Users\USER\.julia\packages\ModelingToolkit\ua9Sp\src\systems\abstractsystem.jl:1714
[2] process_DEProblem(constructor::Type, sys::ODESystem, u0map::Vector{Pair}, parammap::SciMLBase.NullParameters; implicit_dae::Bool, du0map::Nothing, version::Nothing, tgrad::Bool, jac::Bool, checkbounds::Bool, sparse::Bool, simplify::Bool, linenumbers::Bool, parallel::Symbolics.SerialForm, eval_expression::Bool, use_union::Bool, tofloat::Bool, symbolic_u0::Bool, kwargs::Base.Pairs{Symbol, Real, Tuple{Symbol, Symbol, Symbol}, NamedTuple{(:t, :has_difference, :check_length), Tuple{Float64, Bool, Bool}}})
@ ModelingToolkit C:\Users\USER\.julia\packages\ModelingToolkit\ua9Sp\src\systems\diffeqs\abstractodesystem.jl:736
[3] (ODEProblem{true, SciMLBase.AutoSpecialize})(sys::ODESystem, u0map::Vector{Pair}, tspan::Tuple{Float64, Float64}, parammap::SciMLBase.NullParameters; callback::Nothing, check_length::Bool, kwargs::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
@ ModelingToolkit C:\Users\USER\.julia\packages\ModelingToolkit\ua9Sp\src\systems\diffeqs\abstractodesystem.jl:834
[4] ODEProblem
@ C:\Users\USER\.julia\packages\ModelingToolkit\ua9Sp\src\systems\diffeqs\abstractodesystem.jl:827 [inlined]
[5] (ODEProblem{true, SciMLBase.AutoSpecialize})(sys::ODESystem, u0map::Vector{Pair}, tspan::Tuple{Float64, Float64})
@ ModelingToolkit C:\Users\USER\.julia\packages\ModelingToolkit\ua9Sp\src\systems\diffeqs\abstractodesystem.jl:827
[6] (ODEProblem{true})(::ODESystem, ::Vector{Pair}, ::Vararg{Any}; kwargs::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
@ ModelingToolkit C:\Users\USER\.julia\packages\ModelingToolkit\ua9Sp\src\systems\diffeqs\abstractodesystem.jl:814
[7] (ODEProblem{true})(::ODESystem, ::Vector{Pair}, ::Vararg{Any})
@ ModelingToolkit C:\Users\USER\.julia\packages\ModelingToolkit\ua9Sp\src\systems\diffeqs\abstractodesystem.jl:813
[8] ODEProblem(::ODESystem, ::Vector{Pair}, ::Vararg{Any}; kwargs::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
@ ModelingToolkit C:\Users\USER\.julia\packages\ModelingToolkit\ua9Sp\src\systems\diffeqs\abstractodesystem.jl:810
[9] ODEProblem(::ODESystem, ::Vector{Pair}, ::Vararg{Any})
@ ModelingToolkit C:\Users\USER\.julia\packages\ModelingToolkit\ua9Sp\src\systems\diffeqs\abstractodesystem.jl:809
[10] discretize(pdesys::PDESystem, discretization::MOLFiniteDifference{MethodOfLines.CenterAlignedGrid, MethodOfLines.ScalarizedDiscretization}; analytic::Nothing, kwargs::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
@ PDEBase C:\Users\USER\.julia\packages\PDEBase\2OwkC\src\discretization_state.jl:74
[11] discretize(pdesys::PDESystem, discretization::MOLFiniteDifference{MethodOfLines.CenterAlignedGrid, MethodOfLines.ScalarizedDiscretization})
@ PDEBase C:\Users\USER\.julia\packages\PDEBase\2OwkC\src\discretization_state.jl:55
[12] top-level scope
@ .\timing.jl:273
Try using latest versions of all packages, MOL at 0.9.6