From 90ad3ea4d35fff341ce03e4641b8e1d0e89da93c Mon Sep 17 00:00:00 2001 From: chyett home Date: Sat, 13 Jan 2024 14:33:56 -0500 Subject: [PATCH 1/2] fixing tracing logic and adding inhomogenous, nonlinear wave eq test --- src/discretization/staggered_discretize.jl | 10 ++---- test/pde_systems/wave_eq_staggered.jl | 40 +++++++++++++++++++++- 2 files changed, 42 insertions(+), 8 deletions(-) diff --git a/src/discretization/staggered_discretize.jl b/src/discretization/staggered_discretize.jl index 29d65f7b1..36a77f6c0 100644 --- a/src/discretization/staggered_discretize.jl +++ b/src/discretization/staggered_discretize.jl @@ -27,17 +27,13 @@ function symbolic_trace(prob, sys) u2_var = get_var_from_state(states[findfirst(x->get_var_from_state(x)!=u1_var, states)]); u1inds = findall(x->get_var_from_state(x)===u1_var, states); u2inds = findall(x->get_var_from_state(x)===u2_var, states); - u1 = states[u1inds] - u2 = states[u2inds] - tracevec_1 = [i in u2inds ? states[i] : Num(0.0) for i in 1:length(states)] # maybe just 0.0 - tracevec_2 = [i in u1inds ? states[i] : Num(0.0) for i in 1:length(states)] - du1 = prob.f(tracevec_1, nothing, 0.0); - du2 = prob.f(tracevec_2, nothing, 0.0); + tmp = prob.f(states, nothing, 0.0) + du1 = [i in u1inds ? tmp[i] : Num(0.0) for i in 1:length(states)]; + du2 = [i in u2inds ? tmp[i] : Num(0.0) for i in 1:length(states)]; gen_du1 = eval(Symbolics.build_function(du1, states)[2]); gen_du2 = eval(Symbolics.build_function(du2, states)[2]); dynamical_f1(_du1,u,p,t) = gen_du1(_du1, u); dynamical_f2(_du2,u,p,t) = gen_du2(_du2, u); u0 = prob.u0;#[prob.u0[u1inds]; prob.u0[u2inds]]; return SplitODEProblem(dynamical_f1, dynamical_f2, u0, prob.tspan); -#return DynamicalODEProblem{false}(dynamical_f1, dynamical_f2, u0[1:length(u1)], u0[length(u1)+1:end], prob.tspan) end diff --git a/test/pde_systems/wave_eq_staggered.jl b/test/pde_systems/wave_eq_staggered.jl index ed442ff7e..1591b703f 100644 --- a/test/pde_systems/wave_eq_staggered.jl +++ b/test/pde_systems/wave_eq_staggered.jl @@ -19,7 +19,7 @@ using OrdinaryDiffEq bcs = [ρ(0,x) ~ initialFunction(x), ϕ(0.0,x) ~ 0.0, Dx(ρ(t,L)) ~ 0.0, - ϕ(t,-L) ~ 0.0];#-a^2*Dx(ρ(t,L))]; + ϕ(t,-L) ~ 0.0]; domains = [t in Interval(0.0, tmax), x in Interval(-L, L)]; @@ -106,3 +106,41 @@ end @test maximum(sol[1:128,1] .- sol[1:128,test_ind]) < max(dx^2, dt) @test maximum(sol[1:128,1] .- sol[1:128,2*test_ind]) < 10*max(dx^2, dt) end + +@testset "1D inhomogenous nonlinear wave equation, staggered grid, Mixed BC" begin + @parameters x + @variables t, ρ(..), ϕ(..) + Dt = Differential(t); + Dx = Differential(x); + + β = 0.11/(2*0.01); + a = 0.5; + L = 1.0; + dx = 1/2^7; + tspan = (0.0, 6.0) + + function RHS(ρ, ϕ) + return -β*(ϕ^2)*sign(ϕ)/ρ + end + + eqs = [Dt(ρ(t,x)) + Dx(ϕ(t,x)) ~ 0, + Dt(ϕ(t,x)) + a^2*Dx(ρ(t,x)) ~ RHS(ρ(t, x), ϕ(t,x))] + bcs = [ρ(0,x) ~ 50.0*exp(-((x-0.5)*20)^2) + 50.0, + ϕ(0,x) ~ 0.0, + Dx(ρ(t,L)) ~ 0.0, + ϕ(t,0.0) ~ 0.0] + domains = [t in Interval(tspan[1], tspan[end]), + x in Interval(0.0, L)] + + @named pdesys = PDESystem(eqs, bcs, domains, [t,x], [ρ(t,x), ϕ(t,x)]) + discretization = MOLFiniteDifference([x=>dx], t, grid_align=MethodOfLines.StaggeredGrid(), edge_aligned_var=ϕ(t,x)); + prob = discretize(pdesys, discretization); + @time sol = solve(prob, SplitEuler(), dt=(dx/a)^2, adaptive=false) + @test sol.retcode == ReturnCode.Success + # p = plot() + # for i in 1:tspan[end]*((a/dx)^2)/10:length(sol) + # ind = floor(Int, i) + # plot!(p, sol.u[ind], label="t = $(@sprintf("%.2f", sol.t[ind]))") + # end + # p +end From f60648653292ff7ae84eb13918479e727b6c65fb Mon Sep 17 00:00:00 2001 From: chyett home Date: Sat, 13 Jan 2024 19:49:03 -0500 Subject: [PATCH 2/2] dummy commit for CI --- src/discretization/staggered_discretize.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/discretization/staggered_discretize.jl b/src/discretization/staggered_discretize.jl index 36a77f6c0..a0491532c 100644 --- a/src/discretization/staggered_discretize.jl +++ b/src/discretization/staggered_discretize.jl @@ -4,7 +4,7 @@ function SciMLBase.discretize(pdesys::PDESystem, sys, tspan = SciMLBase.symbolic_discretize(pdesys, discretization) try simpsys = structural_simplify(sys) - if tspan === nothing + if tspan === nothing add_metadata!(get_metadata(sys), sys) return prob = NonlinearProblem(simpsys, ones(length(simpsys.states)); discretization.kwargs..., kwargs...)