Write in the handbook the steps and implementation details to take for the solution of the plain vanilla ROM.

Find the right way to compute random vectors in the parameter space within the line defined by the Dirichlet boundary condition coefficient, omega * delta / a0.

• Algorithm.
• Sketch.

Once we have #2 and #5 we can build a ROM of our problem.
This is done by projecting the FOM operators at each step and solving the resultant ROM system of equations.
Once this is done, the reduced solution needs to be projected back into the original space to asses the error.

Warning: How does it work with boundary conditions?

• Build ROM system.
• Dense solver.
• Reconstruct solution in the FOM mesh.

• PB: Product Backlog
• FOM: Full Order Model
• ROM: Reduced Order Model
• DEIM: Discrete Empirical Interpolation Method
• MDEIM: Matrix Discrete Empirical Interpolation Method
• TH: Theoretical Analysis

See Chapters 3 and 5 from the FEniCS tutorial, volume I.
https://fenicsproject.org/tutorial/

The goal of this task is to touch base with the FEniCS solver with an specific application of the problem we aim to reduce.

• #7 Set up reproducible environment steps for FEniCS use.
• Plain vanilla equation.
• Introduce physical parameters.
• Wrap into OOC.

Useful references

• hplgit/fenics-tutorial#74

• https://www2.karlin.mff.cuni.cz/~hron/fenics-tutorial/convection_diffusion/doc.html

Given the current sinusoidal movement of the piston, its initial velocity is not zero!

L(t)  = L0 * (1 - delta * h(t))
L'(t) = -L0 * delta * h'(t)

h(t)  = sin(wt)
h'(t) = w * cos(wt)

L'(t) = -L0 * delta * w * cos(wt) 
L'(0) = -L0 * delta * w

This is not compatible with a fluid departing from rest.

The fix easy,

L(t)  = L0 * (1 - delta * h(t))
L'(t) = -L0 * delta * h'(t)

h(t)  = 1 - cos(wt)
h'(t) = w * sin(wt)

L'(t) = -L0 * delta * w * sin(wt) 
L'(0) = 0

This has an impact in the definition of the boundary conditions.

• FOM Solver:
  • ¿SUPG required? ¿P2-P1 Taylor Hood elements?
  • ALE Formulation.
• ROM Solver:
  • Local assembly for vector fields.

Bare in mind ... KikeM/romtime#14 (comment)

https://people.math.sc.edu/Burkardt/fenics_src/fenics_src.html

Use #22 over the vector form of a matrix.

• Implement.

Instead of using a parabola, I think it would be better to use a sinusoidal function (still taking of from rest).

This is necessary since the previous results were obtained before KikeM/romtime#25 was fixed.

Additionally, it was done with the heat equation, should be try with the piston equations?

Define a function whose solution is exact so that we can test our implementation and further ROM.

• Define test function.
• Implement and certify.

References

• http://www.math.lmu.de/~michel/lax-milgram

father: #2

How do we know how many elements are necessary to detect the basis.

This is actually a nice discussion, because it is not true that we need as many basis elements as nodal points.

In fact, this helps support our decision to rule out the linear operators from our experiments on how many modes can be removed without affecting too much the solution.

Related to KikeM/romtime#23

This is an good post-processing variable whose difference to compute between the FOM and the ROM models.

Use local assemblers instead of generic solve method.

I will need to invoke and obtain the FOM operator to project it.

Reference: https://fenicsproject.org/pub/tutorial/html/._ftut1018.html#ch:poisson0:linalg

Hyper-reduction levels

• None (Naive implementation)
• Only vectors (DEIM)
• All operators (DEIM & MDEIM)

See KikeM/msc-literature-review#23

If we solve a fixed domain problem with an interior moving mesh, we have simplified model equations, and we can verify our use of the convective terms due to mesh movement.

• Compute mass conservation.

Technical details

At the present moment (0.1.0) I am checking the ROM with the L2 norm against the exact solution.
This is fine for the moment, since I am working with a manufactured problem.
However, in the ROM community there are standard ways to measure how well a model is doing.

• Figure them out.
• Document them.

We remove the mesh velocity terms from the convection matrix and the lifting operator to appreciate its effects.

When such term is not included mass is no longer preserved, proving the necessity of a term that accounts for mesh deformation.

See parameters["form_compiler"]["quadrature_degree"] = 4

from fenics import *

mesh = UnitSquareMesh(32,32)

V = FunctionSpace(mesh,'CG',1)

u = TrialFunction(V)
v = TestFunction(V)

L = inner(grad(u),grad(v))*dx
local_matrices = [assemble_local(u*v*dx, c) for c in cells(mesh)]

References

• https://fenicsproject.org/qa/4511/local-assembly-of-bilinear-form/
• https://fenicsproject.org/olddocs/dolfin/2017.2.0/python/programmers-reference/fem/assembling_local/assemble_local.html

• Remove alpha from the parameter vector (it does not belong to the problem parametrization).
• Run several FOM simulations for different values of it and analyse what changes.
• Read KikeM/msc-literature-review#38

This could be studied together with #45.

Now that I have the naive ROM implementation put up for a fixed domain, the question arises on which task to tackle next.

We have the following unresolved issues, ordered from complex to simple:

• 1. (M)DEIM Implementation
  • #22
  • #23
• 2. Moving Domain: Mesh Deformation
  • #24
  • #46
  • KikeM/msc-literature-review#20
• 3. Moving Domain: Generalised Transformation.
• 4. A posteriori errors, ¿Riesz representation theorem?
• 5. POD Implementation (#5).
• 6. ROM-community error calculations (#20).

1. The very first issue is the one where I feel the most uncertain about.
2. Regarding the second issue, fenics has an ALE.move() method which allegedly displaces the mesh.
   However, I am skeptical it will be easy and straightforward.
3. No problems here, I just need to change the weak form definition (this will probably end up in a bit of refactoring, to deal with several problems elegantly).
4. There seems to a way to compute a posteriori errors via the Riesz representation theorem.
   I need to review how to do that.
5. POD should be easy, I only need to review the papers describing its implementation, and read a couple of posts online.
6. Easy peasy :)

Show the difference between BDF-1 and BDF-2.

• Convergence of outflow value in time.
• Net mass error.
• Error difference between successive discretizations.
• Test with and without nonlinear term.

For the ROM too?

• MFP-1
• MFP-2

References:

• asaithambi1992.pdf
• mathematics-09-00749.pdf
• article_Madrid_2_LV.pdf

Once #2 is done, we can generate a thin-layer around the FOM solver to obtain the FOM snapshots (solution, matrices, vectors, etc.).

The goal is to have a clean interface to our solver in order to extract the ingredients for the basis generation.

• Create parameter manager class. Given a parametric space, for the moment define a random generator of snapshots.
• Create snapshot generator class.

Set up a reproducible environment with FEniCS.

Implementation details with the boundary conditions arise when the basis is built.
Ideally, the basis should be built for homogeneous/natural boundary conditions.

• An "easy" workaround is to solve the heat equation for natural boundary conditions du/dn = 0 for the whole domain (although I am not sure if existence is guaranteed).
• Solve for homogeneous Dirichlet boundary conditions for the moment.
  There is no strict need to parametrize the boundary conditions for our purposes.
• The generic implementation is to introduce a lifting function.

• Look in the internet.
  • https://rtd.pymor.org/en/0.5.2/_modules/pymor/algorithms/ei.html#deim
• Implement.

Fixes #25.

Although in theory one can use the solution snapshots as a basis to project the operators, in practice these might define a basis with a large condition number.
Having it go through the POD, even if it is with the same number of basis vectors as snapshots, will produce an orthonormal basis.

Warning: watch-out for standardisations, this could introduce problems.

• Use SVD from scipy.
  • See https://github.com/KikeM/msc-thesis/issues/51#issuecomment-866179711
• Tests

Deliverable
Build POD space of solution snapshots.

References

• https://www.youtube.com/watch?v=axfUYYNd-4Y
• https://rtd.pymor.org/en/0.5.2/_modules/pymor/algorithms/pod.html#pod
• https://scikit-learn.org/stable/modules/generated/sklearn.decomposition.PCA.html
• https://github.com/users/KikeM/projects/2#card-54680317
• https://www.janheiland.de/post/fenics-burger-pod-dmd/
• https://deustotech.github.io/DyCon-Blog/tutorial/wp04/P0009

For those for which there is room for cutting off the basis, do so and show its effect.

Intuitive check that the implementation is correct.

Can I automate this in a test?

There are two degrees of freedom:

• Number of basis elements for RB space.
• Number of basis elements for nonlinear term.

Feedback from Steven:

• Log-log plots of the final error vs. timestep and number of nodes.
• Error vs. mesh velocity plots.
• #54

Originally posted by @KikeM in #48 (comment)

Once #2 is finished, we introduce a moving domain.

The goal of this task is reduce the problem with a non-linearity in the matrices and vectors.

• Reference solution.
• How is it different from the oscillating cylinder problem?

Bare in mind ... KikeM/romtime#14 (comment)

When I was implementing the DEIM procedure for the lifting term, it would not be able to build a basis for one of the tree branches for a certain combination of parameters, namely one where the stationary solution was quickly attained.

I proceeded to inspect the snapshots, and I found (with a certain horror and underlying feeling of a possible bug) that they were identically zero.
Yet, the tests for the FOM solver are passing, which means that the lifting term is being correctly assembled.

What is going on then?

Well, this behaviour is actually not buggy at all, but rather a consequence of the asymptotic stationary nature of MFP-1 (which is one of its desired characteristics).
Indeed, as the solution converges towards a stationary value,

• the forcing term tends towards a constant value (the stationary value is a parabola, hence constant second order derivative),
• and the lifting terms tends towards zero!

Why is this happening?
Well, if we have reached a stationary value, the lifting term only needs to enforce the constant boundary values, which are actually suppressed in the imposition of the homogeneous boundary conditions in the assembly and solution of the algebraic problem.

This means that the lifting term is identically zero!

How can we overcome this issue for asymptotic stationary problems?
By reducing together the forcing term and the lifting term.

• https://fenicsproject.org/qa/10386/impossible-operations-on-matrices-objects-a-b-a-t/

https://fenicsproject.org/qa/8104/generate-sparse-stiffness-matrix-associated-bilinear-choice/

• Search for analytical solutions with a moving domain.
• Define Manufactured Solution.
• Actual implementation.

References

• Stabilisation examples #33
• http://www.dolfin-adjoint.org/en/latest/documentation/time-dependent-wave/time-dependent-wave.html

Once the ROM is defined in terms of the weak formulations (mass, stiffness, forcings, etc.), boundary conditions, parameter space, and hyper reductions, it would be good to automate the "reporting".

Questions to answer in a ROM reporting:

• Which problem is solved.
• Efficiency rates: how much was the problem compressed.
• Errors: L2 (done), H1.
  • Time integration errors.
  • Take a look at the paper from Negri, he suggested a couple of metrics.
  • A posteriori errors.
• Error location: is it near a non-linearity?
• Spectral decay in the POD/SVD.
• How far away can we go from the trained set? (This actually proves if we have captured the dynamics or not).
  • Euclidean distance from the furthest point of the training set.
  • Euclidean distance from the "centre" of the cluster.
• Performance rates: how much time is saved up (this metric is implementation-dependent, be careful).

Goal: show the reduction works with an easy-to-interpret non-linear operator.

• FOM Solver:
  • Integration scheme.
  • Stabilisation terms.
  • #55
  • #56
• ROM Solver:
  • Hyper-reduction of a solution-dependent non-linearity.

Suggestions from Steven:

• Piston with a moving domain.
• Interacting waves (here the non-linearity is a bit more complex to "feel")

FeniCs:

• https://github.com/SourangshuGhosh/burgers_time_viscous/blob/master/burgers_time_viscous.py

The goal is to build a simple code base that solves the reduction problem completely for a fixed domain.

Is this naive approach, the Discrete Empirical Interpolation is not used, since it would only add extra layers of complexity.
The ROM operators are obtained by assembly and projection of the FOM operators.

I am aware this does not honor the offline-online split principle, but it does allow me to tackle all the implementation details.
Once these are dealt with, implementing the (M)DEIM algorithms should not be too cumbersome.

## Steps

• Definition of Manufactured Solution #10.
• Parametrized FOM Solver #7 #2 #9 #15.
• Basis construction: POD SVD orthogonalisation #4 #5.
• Parametrized ROM Solver #8 #6.