dorie issueshttps://ts-gitlab.iup.uni-heidelberg.de/dorie/dorie/-/issues2020-10-22T13:53:13+02:00https://ts-gitlab.iup.uni-heidelberg.de/dorie/dorie/-/issues/196flux reconstruction and transport issues in simplex grid2020-10-22T13:53:13+02:00Hannes Bauserflux reconstruction and transport issues in simplex grid
### Summary
I tried running richards+transport and could not get it working.
The cases I tried contain a rain event with solute concentration. DORiE is able to provide the Richards solution, but not the the transport solution and stops in the first time step of the rain event.
EDIT: *After Discussion with @lriedel, I could get the simulations running if I allow negative concentrations (`[transport.solutionCheck] negativeConcentration = warn`). However, the possible issue with the flux reconstruction remains.*
I have compiled two examples ('2D_transport' and '3D_transport') that show the issues I had.
It seems like part of the issue comes from the flux reconstruction.
In simplex grids, the reconstructed flux seems to depend on the form of the grid. There are systematic differences between the 'shallow' and the 'deep' part of a grid cell at the surface. This can possibly lead to issues with the transport.
EDIT: *I have attached an image showing the issue for the 3D case.![flux_reconstruction](/uploads/cb8ab645891c96bf618206026dcdda4a/flux_reconstruction.png)*
However, I also tried an example with a regular grid ('2d_transport_regular'). This case also fails (but the transport is at least able to make a few time steps). In this case also the resolution is sufficient and should not be an issue. Possibly I am making bad decision in the config file for the transport (however, they are based on the dorie create files). The only way I was able to get this working was by choosing FE order 0.
EDIT: *Allowing negative concentrations resolves the regular grid issues. No obvious flux reconstruction issues are visible here.*
### Steps to reproduce
Please run the provided cases. They all fail. Please have especially a look at the flux reconstruction in '2D_transport' and '3D_transport', which seems rather wrong to me.
### What is the current _bug_ behaviour?
Transport fails. Reconstructed flux at surface does not follow boundary condition.
EDIT: *Transport does not fail when allowing negative concentrations.*
### What is the expected _correct_ behaviour?
Transport works. Reconstructed flux is correct and follows boundary condition.
#### Reproducing input
[2D_transport.zip](/uploads/159c84be524081cf8b55089b5da468b4/2D_transport.zip)
[2D_transport_regular.zip](/uploads/3cb6e534e561ccb67ae095b4e332e464/2D_transport_regular.zip)
[3D_transport.zip](/uploads/2196645642a1aae5fa4e34f0305db5e1/3D_transport.zip)
### Ideas how to fix this?
_Add them here, if you have any._
### Summary
I tried running richards+transport and could not get it working.
The cases I tried contain a rain event with solute concentration. DORiE is able to provide the Richards solution, but not the the transport solution and stops in the first time step of the rain event.
EDIT: *After Discussion with @lriedel, I could get the simulations running if I allow negative concentrations (`[transport.solutionCheck] negativeConcentration = warn`). However, the possible issue with the flux reconstruction remains.*
I have compiled two examples ('2D_transport' and '3D_transport') that show the issues I had.
It seems like part of the issue comes from the flux reconstruction.
In simplex grids, the reconstructed flux seems to depend on the form of the grid. There are systematic differences between the 'shallow' and the 'deep' part of a grid cell at the surface. This can possibly lead to issues with the transport.
EDIT: *I have attached an image showing the issue for the 3D case.![flux_reconstruction](/uploads/cb8ab645891c96bf618206026dcdda4a/flux_reconstruction.png)*
However, I also tried an example with a regular grid ('2d_transport_regular'). This case also fails (but the transport is at least able to make a few time steps). In this case also the resolution is sufficient and should not be an issue. Possibly I am making bad decision in the config file for the transport (however, they are based on the dorie create files). The only way I was able to get this working was by choosing FE order 0.
EDIT: *Allowing negative concentrations resolves the regular grid issues. No obvious flux reconstruction issues are visible here.*
### Steps to reproduce
Please run the provided cases. They all fail. Please have especially a look at the flux reconstruction in '2D_transport' and '3D_transport', which seems rather wrong to me.
### What is the current _bug_ behaviour?
Transport fails. Reconstructed flux at surface does not follow boundary condition.
EDIT: *Transport does not fail when allowing negative concentrations.*
### What is the expected _correct_ behaviour?
Transport works. Reconstructed flux is correct and follows boundary condition.
#### Reproducing input
[2D_transport.zip](/uploads/159c84be524081cf8b55089b5da468b4/2D_transport.zip)
[2D_transport_regular.zip](/uploads/3cb6e534e561ccb67ae095b4e332e464/2D_transport_regular.zip)
[3D_transport.zip](/uploads/2196645642a1aae5fa4e34f0305db5e1/3D_transport.zip)
### Ideas how to fix this?
_Add them here, if you have any._Santiago Ospina De Los Ríossospinar@gmail.comSantiago Ospina De Los Ríossospinar@gmail.comhttps://ts-gitlab.iup.uni-heidelberg.de/dorie/dorie/-/issues/175Introduce Model::pre_step() method to set boundary condition times independen...2020-01-10T15:39:53+01:00Lukas Riedelmail@lukasriedel.comIntroduce Model::pre_step() method to set boundary condition times independently from step computationIn the coupled model, the flux reconstruction at the begin of the step uses the previous boundary conditions, but that currently does not matter because we apply the `NextStep` policy when evaluating it.
### Description
During the `Model::step()` method, the Richards model applies the current start time to its boundary condition manager to retrieve the correct BCs during the step. *Before* calling `step()`, the boundary condition manager therefore applies the boundary conditions of the *previous* step.
### Proposal
Extend the base `Model` interface with a virtual `pre_step()` method which does nothing by default. Override this method in the Richards model to apply the current time to the boundary condition manager. Additionally, the flux reconstruction cache can be reset in this method rather than at the end of `step()`.
Within the coupled model, call `pre_step()` on the Richards model *before* fetching the water flux reconstruction, such that the current boundary conditions are correctly applied. Then, call `step()` and retrieve the second reconstruction. See above for an
Additionally, raise a warning in the Richards model if `step()` is called without calling `pre_step()` before.
```c++
void ModelRichardsTransportCoupling<Traits>::step()
{
// ... //
// NOTE: Proposal. Make sure the boundary condition manager returns the BCs for this time step
//_richards->pre_step();
// set initial state of the water flux to container
auto gf_water_flux_begin = _richards->get_water_flux_reconstructed();
igf_water_flux->push(gf_water_flux_begin,time_begin);
// ^----- Currently applies the old BC for this time
// always do an step of richards
_richards->step();
// ^----- Currently enables the correct BCs for this step
// ... //
// set final state of the water flux to container
auto gf_water_flux_end = _richards->get_water_flux_reconstructed();
igf_water_flux->push(gf_water_flux_end,time_end);
// ^----- Applies the correct BC for this time
// ... //
```
### How to test the implementation?
* Extend `test-model-base` to check if `pre_step()` is called during `run()`.
### Related issuesIn the coupled model, the flux reconstruction at the begin of the step uses the previous boundary conditions, but that currently does not matter because we apply the `NextStep` policy when evaluating it.
### Description
During the `Model::step()` method, the Richards model applies the current start time to its boundary condition manager to retrieve the correct BCs during the step. *Before* calling `step()`, the boundary condition manager therefore applies the boundary conditions of the *previous* step.
### Proposal
Extend the base `Model` interface with a virtual `pre_step()` method which does nothing by default. Override this method in the Richards model to apply the current time to the boundary condition manager. Additionally, the flux reconstruction cache can be reset in this method rather than at the end of `step()`.
Within the coupled model, call `pre_step()` on the Richards model *before* fetching the water flux reconstruction, such that the current boundary conditions are correctly applied. Then, call `step()` and retrieve the second reconstruction. See above for an
Additionally, raise a warning in the Richards model if `step()` is called without calling `pre_step()` before.
```c++
void ModelRichardsTransportCoupling<Traits>::step()
{
// ... //
// NOTE: Proposal. Make sure the boundary condition manager returns the BCs for this time step
//_richards->pre_step();
// set initial state of the water flux to container
auto gf_water_flux_begin = _richards->get_water_flux_reconstructed();
igf_water_flux->push(gf_water_flux_begin,time_begin);
// ^----- Currently applies the old BC for this time
// always do an step of richards
_richards->step();
// ^----- Currently enables the correct BCs for this step
// ... //
// set final state of the water flux to container
auto gf_water_flux_end = _richards->get_water_flux_reconstructed();
igf_water_flux->push(gf_water_flux_end,time_end);
// ^----- Applies the correct BC for this time
// ... //
```
### How to test the implementation?
* Extend `test-model-base` to check if `pre_step()` is called during `run()`.
### Related issueshttps://ts-gitlab.iup.uni-heidelberg.de/dorie/dorie/-/issues/139Flux reconstruction yields large error on 3D simplices and 2nd order polynomi...2020-10-12T19:04:23+02:00Lukas Riedelmail@lukasriedel.comFlux reconstruction yields large error on 3D simplices and 2nd order polynomials (only)The following discussion from !105 should be addressed:
- [ ] @sospinar started a [discussion](https://ts-gitlab.iup.uni-heidelberg.de/dorie/dorie/merge_requests/105#note_20006): (+13 comments)
> I'm having troubles with the $`\mathcal{RT}_1`$ elements in 3D. For simplices, the jumps are quite high, while for cubes a local matrix is throwing an error saying that is singular. Not really sure how to proceed here...
### Collected info
A single testing case for the flux reconstruction shows strongly increased flux jump residuals:
| case | error |
| ------ | ------ |
| 2D_1_cube | < 1E-17 |
| 2D_2_cube | < 1E-17 |
| 2D_3_cube | < 1E-17 |
| 2D_1_simplex | < 1E-17 |
| 2D_2_simplex | < 1E-17 |
| 2D_3_simplex | < 1E-17 |
| 3D_1_cube | < 1E-17 |
| 3D_2_cube | < 1E-17 |
| 3D_3_cube | none |
| 3D_1_simplex | < 1E-17 |
| **3D_2_simplex** | **2.70125e-11** |
| 3D_3_simplex | < 1E-17 |
@sospinar:
>>>
We have several sources of numerical inaccuracies
* In the [DG cases](https://ts-gitlab.iup.uni-heidelberg.de/dorie/dorie/merge_requests/105#note_20066), they come because values that should represent the same value are generated by different local matrices; `Ax=b` problems.
* In the evaluation of the fluxes, they come from the Piola transformation, which have the following operations:
```c++
auto J = e.geometry().jacobianInverseTransposed(x); //! Geometry dependent
J.invert();
J.umtv(x,y); //! y += A^T x
y /= J.determinant();
```
>>>
### Proposal
@sospinar:
> Will forward this question to Prof. Bastian.
### How to test the implementation?
### Related issues
See #65, !105The following discussion from !105 should be addressed:
- [ ] @sospinar started a [discussion](https://ts-gitlab.iup.uni-heidelberg.de/dorie/dorie/merge_requests/105#note_20006): (+13 comments)
> I'm having troubles with the $`\mathcal{RT}_1`$ elements in 3D. For simplices, the jumps are quite high, while for cubes a local matrix is throwing an error saying that is singular. Not really sure how to proceed here...
### Collected info
A single testing case for the flux reconstruction shows strongly increased flux jump residuals:
| case | error |
| ------ | ------ |
| 2D_1_cube | < 1E-17 |
| 2D_2_cube | < 1E-17 |
| 2D_3_cube | < 1E-17 |
| 2D_1_simplex | < 1E-17 |
| 2D_2_simplex | < 1E-17 |
| 2D_3_simplex | < 1E-17 |
| 3D_1_cube | < 1E-17 |
| 3D_2_cube | < 1E-17 |
| 3D_3_cube | none |
| 3D_1_simplex | < 1E-17 |
| **3D_2_simplex** | **2.70125e-11** |
| 3D_3_simplex | < 1E-17 |
@sospinar:
>>>
We have several sources of numerical inaccuracies
* In the [DG cases](https://ts-gitlab.iup.uni-heidelberg.de/dorie/dorie/merge_requests/105#note_20066), they come because values that should represent the same value are generated by different local matrices; `Ax=b` problems.
* In the evaluation of the fluxes, they come from the Piola transformation, which have the following operations:
```c++
auto J = e.geometry().jacobianInverseTransposed(x); //! Geometry dependent
J.invert();
J.umtv(x,y); //! y += A^T x
y /= J.determinant();
```
>>>
### Proposal
@sospinar:
> Will forward this question to Prof. Bastian.
### How to test the implementation?
### Related issues
See #65, !105https://ts-gitlab.iup.uni-heidelberg.de/dorie/dorie/-/issues/138Allow flux reconstruction for non-conforming grids2019-01-09T05:08:38+01:00Santiago Ospina De Los Ríossospinar@gmail.comAllow flux reconstruction for non-conforming gridsThe following discussion from !105 should be addressed:
- [ ] @sospinar started a [discussion](https://ts-gitlab.iup.uni-heidelberg.de/dorie/dorie/merge_requests/105#note_19782): (+6 comments)
> I have been thinking about the non-conforming grids for flux reconstruction. It is possible, but it is absolutely needed to do something more than local flux reconstruction (that is what this MR is about):
>
> One has to solve the same problem but globally, and constraint the normal fluxes to be equal on the intersections. This procedure leads to the so-called hanging nodes. This is not needed for the conforming grids because normal fluxes are trivially equal for RT elements of the same order. In such a case, each grid entity has all the information to do the flux reconstruction, and computations become local. (Notice that p-adaptivity for the flux reconstruction also needs to solve the global problem).
>
> In that case, some of the things implemented here need to be modified:
>
> * Form a global linear system rather than a local linear system. It means that the Raviart Thomas local engine has to be modified so that it forms the right linear system. (Or just creating another local engine).
> * Local Function Spaces (LFS) cannot be the one I implemented here: It can be either that the LFS of PDELab receives intersections, or that the `MinimalLocalFunctionSpace` can construct the DOF and can form a composite space (yeah, it implies using `TypeTree`).
> * One has to constrain the degrees of freedom associated to the entities of codimension 1 in a very similar way as for hanging nodes in normal FEM (PDELab has a quite old code for it, not sure if it would work).
> * As well as for this MR, the lifting can be done separately (e.g. with this MR infrastructure), or directly together with the flux reconstruction RT elements.
> * Check again the computations in the local operator for the non-conforming Raviart Thomas finite elements (and its test functions).
>
> Since it is quite a lot of work I won't do what I just described and since local flux reconstruction in non-conforming grids is not useful at all for us, and I will disable local flux reconstruction computation for the cube-adaptative cases. @lriedel Are you fine with it?The following discussion from !105 should be addressed:
- [ ] @sospinar started a [discussion](https://ts-gitlab.iup.uni-heidelberg.de/dorie/dorie/merge_requests/105#note_19782): (+6 comments)
> I have been thinking about the non-conforming grids for flux reconstruction. It is possible, but it is absolutely needed to do something more than local flux reconstruction (that is what this MR is about):
>
> One has to solve the same problem but globally, and constraint the normal fluxes to be equal on the intersections. This procedure leads to the so-called hanging nodes. This is not needed for the conforming grids because normal fluxes are trivially equal for RT elements of the same order. In such a case, each grid entity has all the information to do the flux reconstruction, and computations become local. (Notice that p-adaptivity for the flux reconstruction also needs to solve the global problem).
>
> In that case, some of the things implemented here need to be modified:
>
> * Form a global linear system rather than a local linear system. It means that the Raviart Thomas local engine has to be modified so that it forms the right linear system. (Or just creating another local engine).
> * Local Function Spaces (LFS) cannot be the one I implemented here: It can be either that the LFS of PDELab receives intersections, or that the `MinimalLocalFunctionSpace` can construct the DOF and can form a composite space (yeah, it implies using `TypeTree`).
> * One has to constrain the degrees of freedom associated to the entities of codimension 1 in a very similar way as for hanging nodes in normal FEM (PDELab has a quite old code for it, not sure if it would work).
> * As well as for this MR, the lifting can be done separately (e.g. with this MR infrastructure), or directly together with the flux reconstruction RT elements.
> * Check again the computations in the local operator for the non-conforming Raviart Thomas finite elements (and its test functions).
>
> Since it is quite a lot of work I won't do what I just described and since local flux reconstruction in non-conforming grids is not useful at all for us, and I will disable local flux reconstruction computation for the cube-adaptative cases. @lriedel Are you fine with it?https://ts-gitlab.iup.uni-heidelberg.de/dorie/dorie/-/issues/132Support restarting of simulations2019-05-17T19:15:43+02:00Santiago Ospina De Los Ríossospinar@gmail.comSupport restarting of simulations### Description
From discussion at #122:
> @sospinar: **If we write `vtk`, it is logical that we must be able to read them.** It would, for example, allow us to set checkpoints on simulations or just use the end of some simulations to start a new one without any further postprocessing.
> @lriedel: I generally agree. Restarting from a written simulation would be nice to have. But there are two major issues:
>
>1. VTK is a really huge library. It does not only contain data I/O, but (mostly) functions for data visualization (!) and analysis. It is now available through most common package managers, but I can understand that the DUNE devs did not want to include it because it will be unavailable on most clusters.
>
>2. The written VTK file does not include the entire GFS information. This is a particular problem with local grid refinement, where the grid configuration cannot—at least not simply—be inferred from the VTK output. Additionally, one may choose cell or vertex data output, and different levels of subsampling. If this is evaluated as initial condition, one never has the guarantee that the assembled DOF vector is the same as in the simulation run which generated the VTK file.
> @sospinar: Well, that's better than nothing. Hopefully, subsampling can support part of the error one makes when writing into `vtk`s.
> @lriedel: But that's exactly my point. Then the restart option is no real restart, but a simulation that starts from an initial bcondition which is _somewhat_ similar to an old solution. I'm not sure if this is really useful.
>
> I see two general options for a program restart:
>
> 1. Use program checkpoints. At a checkpoint, the full precision DOF vector is written to a file. The grid is stored in the [DUNE grid format](https://www.dune-project.org/doxygen/2.6.0/group__DuneGridFormatParser.html). However, I'm unsure about its capabilities and have never worked with it.
>
> 2. Use a solution from a previous program run as initial condition. This is the cheap way. Grid configurations don't need to match and the old solution is interpolated onto the new grid function space.
But that's better discussed in another Issue.
### Proposal
### How to test the implementation?
### Related issues
See #...
<!--
PLEASE READ THIS
Briefly explain __what__ should be changed and __propose__ how this can happen.
Adding pseudo code or diagrams would be great!
Additionally, you can:
- add suitable labels
- assign a milestone
- mention other issues
-->### Description
From discussion at #122:
> @sospinar: **If we write `vtk`, it is logical that we must be able to read them.** It would, for example, allow us to set checkpoints on simulations or just use the end of some simulations to start a new one without any further postprocessing.
> @lriedel: I generally agree. Restarting from a written simulation would be nice to have. But there are two major issues:
>
>1. VTK is a really huge library. It does not only contain data I/O, but (mostly) functions for data visualization (!) and analysis. It is now available through most common package managers, but I can understand that the DUNE devs did not want to include it because it will be unavailable on most clusters.
>
>2. The written VTK file does not include the entire GFS information. This is a particular problem with local grid refinement, where the grid configuration cannot—at least not simply—be inferred from the VTK output. Additionally, one may choose cell or vertex data output, and different levels of subsampling. If this is evaluated as initial condition, one never has the guarantee that the assembled DOF vector is the same as in the simulation run which generated the VTK file.
> @sospinar: Well, that's better than nothing. Hopefully, subsampling can support part of the error one makes when writing into `vtk`s.
> @lriedel: But that's exactly my point. Then the restart option is no real restart, but a simulation that starts from an initial bcondition which is _somewhat_ similar to an old solution. I'm not sure if this is really useful.
>
> I see two general options for a program restart:
>
> 1. Use program checkpoints. At a checkpoint, the full precision DOF vector is written to a file. The grid is stored in the [DUNE grid format](https://www.dune-project.org/doxygen/2.6.0/group__DuneGridFormatParser.html). However, I'm unsure about its capabilities and have never worked with it.
>
> 2. Use a solution from a previous program run as initial condition. This is the cheap way. Grid configurations don't need to match and the old solution is interpolated onto the new grid function space.
But that's better discussed in another Issue.
### Proposal
### How to test the implementation?
### Related issues
See #...
<!--
PLEASE READ THIS
Briefly explain __what__ should be changed and __propose__ how this can happen.
Adding pseudo code or diagrams would be great!
Additionally, you can:
- add suitable labels
- assign a milestone
- mention other issues
-->https://ts-gitlab.iup.uni-heidelberg.de/dorie/dorie/-/issues/118use SubGrid to be able to have diferent grid views in each model2019-07-11T23:35:47+02:00Santiago Ospina De Los Ríossospinar@gmail.comuse SubGrid to be able to have diferent grid views in each model### Description
Use metagrid [`dune-subgrid`](http://numerik.mi.fu-berlin.de/dune-subgrid/index.php) so that we are able to have diferent grid views in each model.
### Proposal
The idea is that each model is able to adapt the grid to its own puposes while stille sharing the same host grid with the whole simulation. `dune-subgrid` would give us the possibility to create our custom grid views for each model and even giving us some support for the markin of the grid. It would remain to us to create the adapters that can transfer information between entities of the two sub grids.
### How to test the implementation?
Each model is still able to compute steps in the new grid views.
<!--
PLEASE READ THIS
Briefly explain __what__ should be changed and __propose__ how this can happen.
Adding pseudo code or diagrams would be great!
Additionally, you can:
- add suitable labels
- assign a milestone
- mention other issues
-->### Description
Use metagrid [`dune-subgrid`](http://numerik.mi.fu-berlin.de/dune-subgrid/index.php) so that we are able to have diferent grid views in each model.
### Proposal
The idea is that each model is able to adapt the grid to its own puposes while stille sharing the same host grid with the whole simulation. `dune-subgrid` would give us the possibility to create our custom grid views for each model and even giving us some support for the markin of the grid. It would remain to us to create the adapters that can transfer information between entities of the two sub grids.
### How to test the implementation?
Each model is still able to compute steps in the new grid views.
<!--
PLEASE READ THIS
Briefly explain __what__ should be changed and __propose__ how this can happen.
Adding pseudo code or diagrams would be great!
Additionally, you can:
- add suitable labels
- assign a milestone
- mention other issues
-->https://ts-gitlab.iup.uni-heidelberg.de/dorie/dorie/-/issues/102Specify output intervals2019-01-17T07:30:58+01:00Lukas Riedelmail@lukasriedel.comSpecify output intervals### Description
DORiE currently writes VTK output for every time step, or none at all. Users might want to specify output intervals for retrieving data and comparing them to measurements, which are typically conducted at certain time intervals, too.
### Proposal
Add a config key `output.interval`, specifying the output interval in seconds. This requires the `CalculationController` to adjust the time step to these intervals. `output.interval = 0` leads to the current output behavior.
### How to test the implementation?
Write a test which
* counts the number of VTK files if `output.interval` is used
* verifies that the VTK files have the expected time stamps
### Related issues
None.### Description
DORiE currently writes VTK output for every time step, or none at all. Users might want to specify output intervals for retrieving data and comparing them to measurements, which are typically conducted at certain time intervals, too.
### Proposal
Add a config key `output.interval`, specifying the output interval in seconds. This requires the `CalculationController` to adjust the time step to these intervals. `output.interval = 0` leads to the current output behavior.
### How to test the implementation?
Write a test which
* counts the number of VTK files if `output.interval` is used
* verifies that the VTK files have the expected time stamps
### Related issues
None.