Staggered interpolations

src/discretization/discretize_vars.jl

@@ -280,8 +280,28 @@
 
 
 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])
 

src/discretization/staggered_discretize.jl

@@ -27,17 +27,13 @@
     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

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