Playground: Residual components in DAE case when using SDC-DAE sweeper arised from epsilon-embedding #455
Add this suggestion to a batch that can be applied as a single commit.
This suggestion is invalid because no changes were made to the code.
Suggestions cannot be applied while the pull request is closed.
Suggestions cannot be applied while viewing a subset of changes.
Only one suggestion per line can be applied in a batch.
Add this suggestion to a batch that can be applied as a single commit.
Applying suggestions on deleted lines is not supported.
You must change the existing code in this line in order to create a valid suggestion.
Outdated suggestions cannot be applied.
This suggestion has been applied or marked resolved.
Suggestions cannot be applied from pending reviews.
Suggestions cannot be applied on multi-line comments.
Suggestions cannot be applied while the pull request is queued to merge.
Suggestion cannot be applied right now. Please check back later.
The$\varepsilon$ -embedding can be used to derive reasonable schemes for (semi-explicit) DAEs that are close to stiff ODEs. When we consider a scalar problem of the form $$y'=f(y, z),$$ $$\varepsilon z' = g(y,z),$$ for $0 < \varepsilon \ll 1$ , then setting $\varepsilon=0$ results in the semi-explicit DAE $$y'=f(y, z),$$ $$0 = g(y,z).$$ Thus, when SDC is applied to the ODE we get the scheme $$\mathbf{y}^{k+1} = \mathbf{y}_0 + \Delta t (\mathbf{Q}-\mathbf{Q}_d) f(\mathbf{y}^k,\mathbf{z}^k) + \Delta t \mathbf{Q}_d f(\mathbf{y}^{k+1},\mathbf{z}^{k+1}),$$ $$\varepsilon\mathbf{z}^{k+1} = \varepsilon\mathbf{z}_0 + \Delta t (\mathbf{Q}-\mathbf{Q}_d) g(\mathbf{y}^k,\mathbf{z}^k) + \Delta t \mathbf{Q}_d g(\mathbf{y}^{k+1},\mathbf{z}^{k+1}).$$ Again, setting $\varepsilon=0$ results in a first possible SDC scheme to solve the DAE: $$\mathbf{y}^{k+1} = \mathbf{y}_0 + \Delta t (\mathbf{Q}-\mathbf{Q}_d) f(\mathbf{y}^k,\mathbf{z}^k) + \Delta t \mathbf{Q}_d f(\mathbf{y}^{k+1},\mathbf{z}^{k+1}),$$ $$0 = \Delta t (\mathbf{Q}-\mathbf{Q}_d) g(\mathbf{y}^k,\mathbf{z}^k) + \Delta t \mathbf{Q}_d g(\mathbf{y}^{k+1},\mathbf{z}^{k+1}).$$ which is called $\mathtt{SDC-E}$ or $$\mathbf{y}^{k+1} = \mathbf{y}_0 + \Delta t (\mathbf{Q}-\mathbf{Q}_d) f(\mathbf{y}^k,\mathbf{z}^k) + \Delta t \mathbf{Q}_d f(\mathbf{y}^{k+1},\mathbf{z}^{k+1}),$$ $$0 = g(\mathbf{y}^{k+1},\mathbf{z}^{k+1}),$$ called $\mathtt{SDC-C}$ or
genericImplicitEmbedded
. [E. Hairer, G. Wanner, eq. (1.12) on p. 375] suggested treating the algebraic constraints directly. Thus, another reasonable SDC scheme isgenericImplicitConstrained
. Here, no quadrature is applied to the algebraic parts.This PR contains these both sweepers to solve semi-explicit DAEs. In the playground the residual components are plotted for both sweepers for a new implemented problem$\Delta t$ where the residual components on last collocation node $\tau_M$ are plotted along iterations. I recommend playing around here to see how it does behave. And then in the plots corresponding to
LinearTestDAEMinion
which is a stiff linear one. Moreover, the original SDC-ODE sweepergeneric_implicit
sweeper can also be applied to the DAE problem (where we ignore the unnecessary computations). This can be done as long as thesolve_system
in the problem class does the correct things. Then, when you run the scriptresidual_components.py
it generates plots for different time step sizesgeneric_implicit
's solves you see that the residual in the algebraic component is always large. And the reason here for that is that in compute_residual the residual is not computed correctly. To check this consider thecompute_residual
methods of the two new sweepers in this PR.However, words above are not set in stone. This is only my assumption, since @brownbaerchen already mentioned issues in convergence when using the residual and you also mentioned that @Ouardghi has issues with the residual. This playground is thus for ..playing around (Who would have thought it?) and see how the residual behaves in the different treatments. Moreover, it could serve as basis for discussion when we three will meet for that.
Since this is only a playground there are no tests here, but I'll submit them later for sure. I already tried to make everything clean. I know it is not perfect (as usual) but when transferring the sweepers + problems to the
projects
folder we can discuss how to make it cleaner. For playing around this should be enough I hope.Oh, I forgot to mention that I added
genericImplicitOriginal
which is in principlegeneric_implicit
only overload thecompute_residual
method and having an attributeresidual_components
which is used for logging and plotting in the playground. I didn't want to clutter upgeneric_implicit
..