dorie issueshttps://ts-gitlab.iup.uni-heidelberg.de/dorie/dorie/-/issues2020-07-31T08:58:58Zhttps://ts-gitlab.iup.uni-heidelberg.de/dorie/dorie/-/issues/204endOfTransportStep output policy leads to desynchronization of models if exec...2020-07-31T08:58:58ZLukas Riedelmail@lukasriedel.comendOfTransportStep output policy leads to desynchronization of models if executed in parallel
### Summary
When selecting `transport.output.policy = endOfTransportStep`, and running the application in parallel, the models desynchronize. The transport model may then take time steps which advance further than the Richards model.
### Steps to reproduce
1. Model settings: `simulation.mode = richards+transport`, `transport.output.policy = endOfTransportStep`
2. Run in parallel: `dorie run -p 2 config.ini`
### What is the current _bug_ behaviour?
Transport model advances further in time than Richards model if it outputs every time step and the simulation is run in parallel.
### What is the expected _correct_ behaviour?
Transport model never advances further in time than Richards model.
### Relevant logs, screenshots, files...?
*Will be supplied later*
#### Reproducing input
| Input data | |
| - | - |
| Simulation Case | _Description goes here_ |
| PFG config file | _if any_ |
| Grid mapping file | _if any_ |
| GMSH grid file | _if any_ |
| Boundary Condition file | |
| Parameterization file | |
| Run config file | |
### Ideas how to fix this?
Check how `write_data()` method or `endOfTransportStep` output policy of transport model mess with time.
### Summary
When selecting `transport.output.policy = endOfTransportStep`, and running the application in parallel, the models desynchronize. The transport model may then take time steps which advance further than the Richards model.
### Steps to reproduce
1. Model settings: `simulation.mode = richards+transport`, `transport.output.policy = endOfTransportStep`
2. Run in parallel: `dorie run -p 2 config.ini`
### What is the current _bug_ behaviour?
Transport model advances further in time than Richards model if it outputs every time step and the simulation is run in parallel.
### What is the expected _correct_ behaviour?
Transport model never advances further in time than Richards model.
### Relevant logs, screenshots, files...?
*Will be supplied later*
#### Reproducing input
| Input data | |
| - | - |
| Simulation Case | _Description goes here_ |
| PFG config file | _if any_ |
| Grid mapping file | _if any_ |
| GMSH grid file | _if any_ |
| Boundary Condition file | |
| Parameterization file | |
| Run config file | |
### Ideas how to fix this?
Check how `write_data()` method or `endOfTransportStep` output policy of transport model mess with time.https://ts-gitlab.iup.uni-heidelberg.de/dorie/dorie/-/issues/196flux reconstruction and transport issues in simplex grid2020-08-11T04:47:15ZHannes 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/188Avoid dangling references in lambdas returned by parameterization interface2020-04-27T18:39:18ZLukas Riedelmail@lukasriedel.comAvoid dangling references in lambdas returned by parameterization interface### Overview
The following discussion from !187 should be addressed:
- [ ] @sospinar started a [discussion](https://ts-gitlab.iup.uni-heidelberg.de/dorie/dorie/-/merge_requests/187#note_48281): (+1 comment)
> I have to accept that I never realized of this `this` on the parameterization interface. And now that I do, I do not like it at all. This lambda is getting passed through at least two classes with these references and is shouting "I will hold a dangling reference whenever you don't expect it". It is _never_ advised to export lambdas with closures that hold references to other scopes (see [CppCoreGuidelines](https://github.com/isocpp/CppCoreGuidelines/blob/master/CppCoreGuidelines.md#Rf-value-capture) or [Effective Modern C++
> ](http://shop.oreilly.com/product/0636920033707.do?cmp=af-code-books-video-product_cj_0636920033707_7708709)).
>
> Now is perhaps too late to complain and ask for a redesign as we have several of them, so I would be pleased with just a big warning to use these functions and those in the parameter classes: "only if you know what you are doing" or with very precise instructions how to avoid the dangling references.
### Proposal
In the lambdas returned by `Dorie::Parameterization::Richards` and `::MualemVanGenuchten`, capture the required parameters by value. *Edit: Parameterizations in the Transport model are affected as well.*
### How to test the implementation?
Evaluate a parameterization function after destroying the parameterization interface which supplied it.
### Related issues### Overview
The following discussion from !187 should be addressed:
- [ ] @sospinar started a [discussion](https://ts-gitlab.iup.uni-heidelberg.de/dorie/dorie/-/merge_requests/187#note_48281): (+1 comment)
> I have to accept that I never realized of this `this` on the parameterization interface. And now that I do, I do not like it at all. This lambda is getting passed through at least two classes with these references and is shouting "I will hold a dangling reference whenever you don't expect it". It is _never_ advised to export lambdas with closures that hold references to other scopes (see [CppCoreGuidelines](https://github.com/isocpp/CppCoreGuidelines/blob/master/CppCoreGuidelines.md#Rf-value-capture) or [Effective Modern C++
> ](http://shop.oreilly.com/product/0636920033707.do?cmp=af-code-books-video-product_cj_0636920033707_7708709)).
>
> Now is perhaps too late to complain and ask for a redesign as we have several of them, so I would be pleased with just a big warning to use these functions and those in the parameter classes: "only if you know what you are doing" or with very precise instructions how to avoid the dangling references.
### Proposal
In the lambdas returned by `Dorie::Parameterization::Richards` and `::MualemVanGenuchten`, capture the required parameters by value. *Edit: Parameterizations in the Transport model are affected as well.*
### How to test the implementation?
Evaluate a parameterization function after destroying the parameterization interface which supplied it.
### Related issuesLukas Riedelmail@lukasriedel.comLukas Riedelmail@lukasriedel.comhttps://ts-gitlab.iup.uni-heidelberg.de/dorie/dorie/-/issues/175Introduce Model::pre_step() method to set boundary condition times independen...2020-01-10T14:39:53ZLukas 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-08-05T01:15:06ZLukas 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, !105