[WIP] Better staggered implementation #285

cmhyett · 2023-06-21T22:44:48Z

Working to solve #255. Previous work found a way to perform time staggering - now we want to perform space staggering as well. The first thing I did was encapsulate functionality to more easily identify component parts (and perhaps overwrite them if need be).

Currently, Iaxies & Igrid members of DiscreteSpace are not utilized (at least as best as I can tell), but hold the potential to provide for a mapping between independent and dependent grids.

By default, it's a bijective map, but in general it doesn't need to be. However, since they aren't currently used, it'll take a bit of work to identify correct locations to use them instead of the dependent variable axies & grids. As a basic example, suppose we have u(x,t) => u[1](t) = u(x[1],t), u[2](t) = u(x[3],t), u[3](t) = u(x[5],t) and v(x,t) => v[1](t) = v(x[2],t), v[2](t) = v(x[4],t), v[3](t) = v(x[6],t), and using centered in space finite differences to update ~ du[2] = (v[2] - v[1])/(2*dx) notice that this is actually second-order centered because of the staggering.

This would change the differential operator based on the grid type.....may or may not be viable.

@xtalax @ChrisRackauckas, let me know if you forsee issues here.

… it easier to extend/specialize on grid type

cmhyett · 2023-06-22T14:37:19Z

Having thought about it overnight, I think doing a simple variable transformation as described in the attached image achieves the same effect, with the benefit of having an explicit mapping between the variables.

Still open to utilization of Igrid and Iaxies (should be equivalent, but unsure of side-effects), if @xtalax has ideas.

This is particularly nice because I only have to map the second-order space discrete differential to the first order forward space differential.

…retize, but this shouldn't be necessary in general) to get spatial staggering working. very close now, just need to tease out the details of the discretized differential operators now.

…e available elsewhere

codecov · 2023-06-28T13:59:38Z

Codecov Report

Attention: 156 lines in your changes are missing coverage. Please review.

Comparison is base (a04e094) 83.02% compared to head (0694c0d) 79.62%.

Files	Patch %	Lines
...schemes/centered_difference/centered_difference.jl	0.00%	80 Missing ⚠️
src/discretization/generate_bc_eqs.jl	0.00%	32 Missing ⚠️
src/discretization/staggered_discretize.jl	0.00%	30 Missing ⚠️
...discretization/generate_finite_difference_rules.jl	0.00%	11 Missing ⚠️
src/interface/solution/timedep.jl	60.00%	2 Missing ⚠️
src/interface/MOLFiniteDifference.jl	50.00%	1 Missing ⚠️

Additional details and impacted files

@@            Coverage Diff             @@
##           master     #285      +/-   ##
==========================================
- Coverage   83.02%   79.62%   -3.40%     
==========================================
  Files          40       41       +1     
  Lines        1756     1944     +188     
==========================================
+ Hits         1458     1548      +90     
- Misses        298      396      +98

☔ View full report in Codecov by Sentry.
📢 Have feedback on the report? Share it here.

…ath forward specializing only low-level differential discretizer methods

…cmhyett/MethodOfLines.jl into better_staggered_implementation

cmhyett · 2023-07-13T00:17:01Z

Following up on the commits yesterday, we should be able to combine this and the time-staggering hack to get a stable wave equation method. VelocityVerlet complains unfortunately, so I used another method that I'm unfamiliar with, I don't think it is the issue..

The results are that we obtain some advection, (not the correct advection in both directions unfortunately), with instability. There seems to be an off-by-one error in the flux that I can't quite track down.

cmhyett · 2023-07-13T19:39:38Z

We have wave equation looking results! I implemented my own simple leapfrog update - there is a strange bit about how the symbolic tracing modifies the du variables - check this placeholder in calc_du!

There is also a preference to leftward advection; I think this must be due to an off-by-one error, and I suspect the boundary substitution to be the issue.

For reference, I put the code here including all the hacks necessary to do the time staggering.

using OrdinaryDiffEq, ModelingToolkit, MethodOfLines, DomainSets, Plots;

@parameters t x
@variables ρ(..) ϕ(..)
Dt = Differential(t);
Dx = Differential(x);

a = 1.0;#1.0/2.0;
L = 8.0;
dx = 0.0625;
dt = dx/a;

initialFunction(x) = exp(-(x-dx)^2);
eq = [Dt(ρ(t,x)) + Dx(ϕ(t,x)) ~ 0,
      Dt(ϕ(t,x)) + a^2 * Dx(ρ(t,x)) ~ 0]
bcs = [ρ(0,x) ~ initialFunction(x),
       ϕ(0.0,x) ~ 0.0,
       ρ(t,-L) ~ ρ(t,L),
       ϕ(t,-L) ~ ϕ(t,L)];

domains = [t in Interval(0.0, 10.0),
           x in Interval(-L, L)];

@named pdesys = PDESystem(eq, bcs, domains, [t,x], [ρ(t,x), ϕ(t,x)]);

#discretization = MOLFiniteDifference([x=>dx], t, grid_align=MethodOfLines.CenterAlignedGrid());
discretization = MOLFiniteDifference([x=>dx], t, grid_align=MethodOfLines.StaggeredGrid());

prob = discretize(pdesys, discretization);

len = floor(Int, length(prob.u0)/2);
@variables rho[1:len] phi[1:len]
drho = (prob.f([collect(rho); collect(phi)], nothing, 0.0)[1:len]);
dphi = (prob.f([collect(rho); collect(phi)], nothing, 0.0)[len+1:end]);

gen_drho = eval(Symbolics.build_function(drho, collect(rho), collect(phi))[2]);
gen_dphi = eval(Symbolics.build_function(dphi, collect(rho), collect(phi))[2]);

dynamical_f1(_drho,u,p,t) = gen_drho(_drho, u[1:len], u[len+1:end]);
dynamical_f2(_dphi,u,p,t) = gen_dphi(_dphi, u[1:len], u[len+1:end]);
u0 = [prob.u0[1:len]; prob.u0[len+1:end]];
prob = DynamicalODEProblem(dynamical_f1, dynamical_f2, u0[1:len], u0[len+1:end], (0.0,1.0))
tsteps = 0.0:dt:10.0
#sol = solve(prob, IRKN3(), dt=dt, saveat=tsteps);

function calc_du!(du, u, p, t)
    placeholder = copy(du);
    dynamical_f1(du, u, p, t);
    dynamical_f2(placeholder, [(u[1:len] + dt*du[1:len]); u[len+1:end]], p, t);
    du[len+1:end] .= placeholder[1:len];
    return;
end

function update(u, p, t)
    du = zeros(length(u));
    calc_du!(du, u, p, t);
    return u .+ du*dt;
end

sol = zeros(length(u0), length(tsteps));
sol[:,1] .= u0;
for i in 2:length(tsteps)
    sol[:,i] = update(sol[:,i-1], 0, 0);
end

function plot_sol(sol)
    p_rho = plot();
    p_phi = plot();
    for i in 1:5:20;#length(sol)
        plot!(p_phi, sol[1:len,i], legend=false);
        plot!(p_rho, sol[len+1:end,i], legend=false)
    end
    return plot(p_rho, p_phi);
end

plot_sol(sol)

…to DynamicalODEProblem

…ample code

cmhyett · 2023-08-10T22:50:32Z

@xtalax , looking in the last commit, in staggered_discretize.jl symbolic_trace(), can you help me construct the findall statements?

I want to look for all the same variables, i.e., collect all the \rho in u1 and all the \phi in u2 but I can't figure out how to do the comparison

…correctly and staggered data is available in DiscreteSpace object

…iscretization kwarg

…ith OTS solvers

…to work with revised output from discretize

cmhyett · 2023-09-01T19:46:53Z

In addition to the included tests, a few words about the final commits:

@xtalax and I made the symbolic tracing more robust
We now output a SplitODEProblem which can be solved using standard (RK4) or specialized (CNLF2) solvers
Along with Dirichlet and Neumann BCs, time-derivatives in the boundary equations work - in our tests we implement transparent BCs.

In the attached pictures we show stability of transparent and reflecting boundary conditions - these cases are codified in the tests.

…cmhyett/MethodOfLines.jl into better_staggered_implementation

xtalax · 2023-09-04T15:04:54Z

2 things, We need to get tests passing, and we need a docs page detailing how this is used with an example

cmhyett · 2023-09-04T16:50:25Z

@xtalax , do you understand these errors? The offending lines of code (these get_index calls) are not code that I modified...Errors are in generate_grid for EdgeAlignedGrid

Perhaps I clobbered an include in src/MethodOfLines.jl?

xtalax · 2023-09-04T17:09:50Z

Somehow, you've broken the default construct_discrete_space, or changed the behaviour of MOLFiniteDifference, I suggest you look at master and compare

…known upstream issue

cmhyett · 2023-11-07T22:31:38Z

Integration tests failed due to an upstream failure addressed in:
SciML/SimpleNonlinearSolve.jl#96

@xtalax PR tests are ready to be re-run.

xtalax

src/interface/grid_types.jl

 const center_align=CenterAlignedGrid()
 const edge_align=EdgeAlignedGrid()
+const stagger_align=StaggeredGrid()


please make this staggered_grid

xtalax · 2023-06-28T13:34:09Z

examples/waveEquation.jl

@@ -0,0 +1,37 @@
+using OrdinaryDiffEq, ModelingToolkit, MethodOfLines, DomainSets, Plots;


you should add these to PDESystemLibrary.jl, and tag them appropriately so that you can retrieve them, rather than in a folder here

I think you commented on an old file? This has been moved to docs/staggered.md and it follows the same template as is found in e.g. heatss.md

Can you clarify?

src/discretization/staggered_discretize.jl

xtalax · 2023-08-30T14:23:46Z

src/discretization/generate_bc_eqs.jl

+    end
+    othervars = filter(v -> (length(arguments(v)) != 1) && any(isequal(x_), arguments(depvar(v, s))), othervars)
+
+    depvarderivbcmaps = [(Differential(x_)^d)(u_) => central_difference(derivweights, II, s, [], (x2i(s, u, x_), x_), u, ufunc, d) for d in derivweights.orders[x_]];


remember the depvar/derivbcmaps are different for variables on the other grid

xtalax · 2023-11-17T18:55:02Z

docs/src/tutorials/heatss.md

@@ -1,7 +1,7 @@
 # Steady State Heat Equation - No Time Dependence - NonlinearProblem

 Occasionally, it is desirable to solve an equation that has no time evolution, such as the steady state heat equation:
-```@example heatss
+```


need to change this back

It works now?

xtalax · 2023-11-17T18:55:26Z

test/pde_systems/MOL_NonlinearProblem.jl

@@ -23,7 +23,7 @@ using DomainSets
 end

 # Laplace's Equation, same as above but with MOL discretization
-@testset "1D Laplace - constant solution" begin
+@test_skip begin #@testset "1D Laplace - constant solution" begin


change all these back

cmhyett added 2 commits June 16, 2023 16:50

encapsulating basic functionality of construct_discrete_space to make…

b8ea3ad

… it easier to extend/specialize on grid type

adding staggered grid type

e88c34c

cmhyett changed the title ~~Better staggered implementation~~ [WIP] Better staggered implementation Jun 21, 2023

cmhyett added 2 commits June 22, 2023 11:03

overloading transform_pdesys (and unfortunately for now symbolic_disc…

63d9303

…retize, but this shouldn't be necessary in general) to get spatial staggering working. very close now, just need to tease out the details of the discretized differential operators now.

this commit should be reverted, this is debug code that I want to hav…

e4416d1

…e available elsewhere

cmhyett and others added 3 commits July 11, 2023 12:27

removing symbolic discretize specialization, think we have a better p…

5875ddc

…ath forward specializing only low-level differential discretizer methods

implementing centered finite difference for staggered grid\!

d642f1a

Merge branch 'better_staggered_implementation' of https://github.com/…

22d428d

…cmhyett/MethodOfLines.jl into better_staggered_implementation

cmhyett added 6 commits July 27, 2023 16:29

fixes to resolve off-by-one error in staggered solve of wave equation

9d03ffe

working staggered boundary conditions

5ec416e

specializing SciMLBase.discretize to put staggered discretizations in…

a48c627

…to DynamicalODEProblem

tucking DynamicalODEProblem inside discretize, so removing it from ex…

d5090b9

…ample code

successfully returning DynamicalODEProblem from discretize

5829d72

a bit of massaging, simplyfying, poking at generality

a761517

cmhyett and others added 11 commits August 11, 2023 10:52

removing debug code

abc2053

formulating guiding examples as tests

fd114f1

generalizing symbolic tracing operation

795efc7

framing hack fixing changes

49fdb98

upstream pipelining just about done, need to check constructor works …

c8382de

…correctly and staggered data is available in DiscreteSpace object

DiscreteSpace now holds staggered dictionary created from staggered d…

34dfc1b

…iscretization kwarg

finished removing hacks, put in proper logic, beginning to add tests

c898630

adding staggered tests to runtests

7d72efb

adding additional tests...currently doesn't work

6e81969

updating staggered discretize to output split ode system compatible w…

7ea157e

…ith OTS solvers

removing duplicate/ambiguous code

40ee388

cmhyett added 2 commits August 31, 2023 20:15

correcting transparent boundary conditions, updating solution format …

48ccabe

…to work with revised output from discretize

updating test to use splitode solution interface

9335503

cmhyett added 4 commits September 1, 2023 12:48

Merge branch 'master' into better_staggered_implementation

49009ca

no change, trying to trigger github actions

26a14c5

Merge branch 'better_staggered_implementation' of https://github.com/…

b4a833d

…cmhyett/MethodOfLines.jl into better_staggered_implementation

adding staggered wave equation tests to CI

b68861e

cmhyett and others added 5 commits October 21, 2023 21:53

fixing precompilation bug, vars->vars.intervals

974873d

remedying includes

5cacaf9

adding required includes, correcting test bug

4611d1f

adding docs file for staggered grids

27a15d0

marking nonlinear laplacian tests as broken - directed by Alex for a …

c36019c

…known upstream issue

cmhyett mentioned this pull request Oct 30, 2023

Fix NonlinearProblem tests #334

Closed

removing from heat equation example due to upstream failing code

07f1a05

xtalax closed this Oct 31, 2023

xtalax reopened this Oct 31, 2023

dummy commit to trigger CI

4cd0e51

xtalax closed this Nov 8, 2023

xtalax reopened this Nov 8, 2023

cmhyett and others added 2 commits November 8, 2023 14:21

dummy commit to trigger CI

488cae2

merging with main branch updates

0694c0d

xtalax reviewed Nov 17, 2023

View reviewed changes

xtalax merged commit a3810a0 into SciML:master Nov 28, 2023

cmhyett deleted the better_staggered_implementation branch December 4, 2023 16:06

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Uh oh!

Uh oh!

[WIP] Better staggered implementation #285

[WIP] Better staggered implementation #285

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

		@@ -0,0 +1,37 @@
		using OrdinaryDiffEq, ModelingToolkit, MethodOfLines, DomainSets, Plots;

Uh oh!

[WIP] Better staggered implementation #285

[WIP] Better staggered implementation #285

Uh oh!

Conversation

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Codecov Report

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Uh oh!