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 f5044634d78e8e8c36ed87877587fbde04dec7e8 Mon Sep 17 00:00:00 2001 From: chyett home Date: Sat, 13 Jan 2024 16:52:52 -0500 Subject: [PATCH 2/2] first attempt at variable interpolation. Needs revising, and does not work for periodic BCs, but otherwise seems to achieve goal --- src/discretization/discretize_vars.jl | 22 +++++++++++++++++++++- 1 file changed, 21 insertions(+), 1 deletion(-) diff --git a/src/discretization/discretize_vars.jl b/src/discretization/discretize_vars.jl index c02163df2..5cecf90e6 100644 --- a/src/discretization/discretize_vars.jl +++ b/src/discretization/discretize_vars.jl @@ -280,8 +280,28 @@ gridvals(s::DiscreteSpace{N}, u, I::CartesianIndex) where {N} = ndims(u, s) == 0 varmaps(s::DiscreteSpace, depvars, II, indexmap) = [u => s.discvars[u][Idx(II, s, u, indexmap)] for u in depvars] - valmaps(s::DiscreteSpace, u, depvars, II, indexmap) = length(II) == 0 ? [] : vcat(varmaps(s, depvars, II, indexmap), gridvals(s, u, II)) +function valmaps(s::DiscreteSpace{N,M,G}, u, depvars, II, indexmap) where {N,M,G<:StaggeredGrid} + if (length(II) == 0) + return 0; + else + interpolated_varmap = [] + if (length(depvars) > 1) + if (s.staggeredvars[operation(u)] == CenterAlignedVar) + ii = vcat(II, II+CartesianIndex(1)) + else + ii = vcat(II, II-CartesianIndex(1)) + end + interpolate(vars) = (vars[1] + vars[2])/2; #take average for now...can/should we generalize? + interpolated_var = depvars[findfirst(x->operation(x)!==operation(u), depvars)] + interpolated_varmap = [interpolated_var=>interpolate([s.discvars[interpolated_var][Idx(i, s, interpolated_var, indexmap)] for i in ii])] + end + + return vcat(varmaps(s, [depvars[findfirst(x->operation(x)==operation(u), depvars)]], II, indexmap), + interpolated_varmap, + gridvals(s, u, II)) + end +end valmaps(s, u, depvars, indexmap) = valmaps.([s], [u], [depvars], s.Igrid[u], [indexmap])