hyteg issueshttps://i10git.cs.fau.de/groups/hyteg/-/issues2024-03-28T10:34:38+01:00https://i10git.cs.fau.de/hyteg/hyteg/-/issues/257Potential tutorial folder structure2024-03-28T10:34:38+01:00Ponsuganth Ilangovan Ponkumar IlangoPotential tutorial folder structureThis issue is to discuss modifying the tutorials folder structure to make it well categorised, like basics, full apps, etc. As we would also like to have the benchmark apps like Blankenbach under the tutorials which would serve as an exa...This issue is to discuss modifying the tutorials folder structure to make it well categorised, like basics, full apps, etc. As we would also like to have the benchmark apps like Blankenbach under the tutorials which would serve as an example of community standard geophysical application.
A potential folder structure would be,
- Tutorials
- Basics
- 01_PrimitiveStorage
- 02_PrimitiveData
- 03_Communication
- 04_Indexing
- Benchmarks
- 13_Blankenbach
- Full Apps
- others
Any suggestions or improvement with this and maybe also ideas for more tutorials would be greatly helpful as we start to work on these now.Ponsuganth Ilangovan Ponkumar IlangoPonsuganth Ilangovan Ponkumar Ilangohttps://i10git.cs.fau.de/hyteg/hyteg/-/issues/256Doxygen issues2024-03-26T16:09:48+01:00Nils KohlDoxygen issuesA collection of current issues in our Doxygen documentation:
- [x] LaTeX does not seem to render, see https://hyteg.pages.i10git.cs.fau.de/hyteg/10_DGAMR.html (see https://i10git.cs.fau.de/hyteg/hyteg/-/issues/256#note_29522)
- [ ] CI ba...A collection of current issues in our Doxygen documentation:
- [x] LaTeX does not seem to render, see https://hyteg.pages.i10git.cs.fau.de/hyteg/10_DGAMR.html (see https://i10git.cs.fau.de/hyteg/hyteg/-/issues/256#note_29522)
- [ ] CI badges in README do not link to CI
- [ ] some links in the README do not work correctly when clicked in GitLabNils KohlNils Kohlhttps://i10git.cs.fau.de/hyteg/hyteg/-/issues/255Add Viscosity profile templates2024-03-14T10:58:25+01:00Eugenio D'AscoliAdd Viscosity profile templateshttps://i10git.cs.fau.de/hyteg/hyteg/-/issues/253Scaled Apply + GEMV2024-03-08T08:26:10+01:00Nils KohlScaled Apply + GEMVThe proposal of scaled operators in !703 initiated this issue.
Having the standard BLAS level 2 routine `gemv` (see e.g. [here](https://spec.oneapi.io/versions/latest/elements/oneMKL/source/domains/blas/gemv.html))
```math
y \gets \alph...The proposal of scaled operators in !703 initiated this issue.
Having the standard BLAS level 2 routine `gemv` (see e.g. [here](https://spec.oneapi.io/versions/latest/elements/oneMKL/source/domains/blas/gemv.html))
```math
y \gets \alpha A x + \beta y
```
as part of the `Operator` interface would be great for obvious reasons.
What we currently can do is
```math
y \gets A x + \gamma y
```
where $\gamma \in \{0, 1\}$ (a.k.a., `Replace` and `Add`).
Implementing the full `gemv` would be the optimal fix, however, since many of the matrix-free operator implementations work in an "additive" fashion (particularly all elementwise operators including the new generated versions) I do not see an obvious way around having to touch $y$ (i.e., the "dst" vector) *twice* to apply scaling with $\beta$. If there is no obvious solution, we are forced to first scale $y$ (if $\beta \neq 0$) and then perform a scaled matrix-vector multiplication after that. We can do this scaling of $y$ inside the operator, but I do not see a real reason to do that inside the operator and not just call `assign` beforehand.
MR !703 proposes an implementation that enables scaled operators
```math
y \gets \alpha A x + \gamma y
```
that is in my opinion not optimal since $\alpha$ becomes part of the operator's state.
To resolve the "state"-issue mentioned above, I propose to actually extend the `apply` routine by a parameter $\alpha$ instead of changing the operator's state:
```
virtual void apply( const ValueType& alpha // new parameter
const SourceFunction& src,
const DestinationFunction& dst,
size_t level,
DoFType flag,
UpdateType updateType = Replace ) const
```
Refactoring this requires some work, but should hopefully be doable. I'd keep the current apply, and just make it call the new one with $\alpha = 1$.
-----
EDIT: We basically already do this:
> We can do this scaling of $y$ inside the operator, but I do not see a real reason to do that inside the operator and not just call `assign` beforehand.
if `UpdateType == Replace` since we [interpolate `0`](https://i10git.cs.fau.de/hyteg/hyteg/-/blob/0a8460a3dd538d9d91a194b5d5588c1a4ef2b1ac/src/hyteg/elementwiseoperators/P1ElementwiseOperator.cpp#L96).
Therefore, for completeness I'll also propose an interface for the `gemv`:
```
virtual void gemv( const ValueType& alpha // new parameter
const SourceFunction& src,
const ValueType& beta // new parameter
const DestinationFunction& dst,
size_t level,
DoFType flag ) const
```
with the perk that it gets rid of the `UpdateType` parameter.Nils KohlNils Kohlhttps://i10git.cs.fau.de/hyteg/hyteg/-/issues/252Elementwise operator applied to parts of domain2024-03-11T07:16:38+01:00Maximilian DechantElementwise operator applied to parts of domainIn some cases it would be useful to apply an operator to only some part of the domain.
In particular, I need to apply a P1ElementwiseOperator to only a few faces that satisfy a certain condition. This is currently done by using a modifi...In some cases it would be useful to apply an operator to only some part of the domain.
In particular, I need to apply a P1ElementwiseOperator to only a few faces that satisfy a certain condition. This is currently done by using a modified apply function that takes such a condition as an argument and then checks the condition for every face.
Maybe this is useful for others as well in which case it could be implemented as a new feature or maybe there is already a method to do this which I am unaware of.
EDIT: A first simple version of the P1ElementwiseOperator that is only applied to some faces can be seen on the branch „dechant+wagneran/LaplaceBeltrami2“https://i10git.cs.fau.de/hyteg/hyteg/-/issues/251Document whether the pseudo-3D setting for 2D FullStokes is applied.2024-03-06T08:41:58+01:00Nils KohlDocument whether the pseudo-3D setting for 2D FullStokes is applied.From [this comment](https://i10git.cs.fau.de/hyteg/hyteg/-/issues/246#note_28573) in #246:
> just a comment on the `P2P1StokesFull` case in 2D. Though, I am not sure how important this will be for us at the moment. The current **form_st...From [this comment](https://i10git.cs.fau.de/hyteg/hyteg/-/issues/246#note_28573) in #246:
> just a comment on the `P2P1StokesFull` case in 2D. Though, I am not sure how important this will be for us at the moment. The current **form_stokes_full.ufl** for the 2D case uses 2/3 as factor, which is a *'pseudo 2D'* setting. That was motivated by people doing 2D simulations, with are actually 3D but with a cylindrical symmetry of sorts.
We should document what is implemented in the 2D case. Probably needs to be done through the HFG.Nils KohlNils Kohlhttps://i10git.cs.fau.de/hyteg/hyteg/-/issues/250Generate missing pieces for all relevant Stokes preconditioners2024-03-06T08:38:03+01:00Nils KohlGenerate missing pieces for all relevant Stokes preconditionersEmerged as a remainder from #246.
Mainly we need the diagonal operators and some lumping as well as [some smoke tests](https://en.wikipedia.org/wiki/Smoke_testing_(software)) in HyTeG.Emerged as a remainder from #246.
Mainly we need the diagonal operators and some lumping as well as [some smoke tests](https://en.wikipedia.org/wiki/Smoke_testing_(software)) in HyTeG.Nils KohlNils Kohlhttps://i10git.cs.fau.de/hyteg/hyteg/-/issues/249Gradient of FE functions2024-03-11T11:26:50+01:00Ponsuganth Ilangovan Ponkumar IlangoGradient of FE functionsThis issue was to discuss if there are already some implementation to evaluate the gradient of an FE function, one related function I could find was [evaluateGradient](https://i10git.cs.fau.de/hyteg/hyteg/-/blob/master/src/hyteg/p2functi...This issue was to discuss if there are already some implementation to evaluate the gradient of an FE function, one related function I could find was [evaluateGradient](https://i10git.cs.fau.de/hyteg/hyteg/-/blob/master/src/hyteg/p2functionspace/P2Function.hpp#L136) but this is only for 2D and I am not sure if it is tested well also.
For example, if we could evaluate the gradient of a scalar FE function and get the result as a vector FE function. There would be multiple issues to deal with like, what should we do on (faces, ) edges and nodes and how efficiently we perform the evaluation.
I feel this would be useful when we go towards nonlinear rheologies (viscosity) and also maybe for some post-processing in geophysics, like studying the stress field from the velocity field.Ponsuganth Ilangovan Ponkumar IlangoPonsuganth Ilangovan Ponkumar Ilangohttps://i10git.cs.fau.de/hyteg/hyteg/-/issues/247Add template viscosity profiles as input data for convection simulations2024-03-18T11:24:26+01:00Eugenio D'AscoliAdd template viscosity profiles as input data for convection simulationsFor convection simulations with non-constant viscosities (in radial direction or temperature dependent) a viscosity profile must be read in and stored in the corresponding data container.
Currently, no "standard" template viscosity prof...For convection simulations with non-constant viscosities (in radial direction or temperature dependent) a viscosity profile must be read in and stored in the corresponding data container.
Currently, no "standard" template viscosity profiles are available.
The idea is to provide a set of template viscosity profiles for simulation runs. These template profiles could be stored under ```hyteg/data/terraneo/viscosityProfiles```.
@HamishBrown @IsabelPapanagnou @kohlHamish BrownHamish Brownhttps://i10git.cs.fau.de/hyteg/hyteg/-/issues/245Discussion: Handling Boundary Conditions in Vector-Valued Functions2024-02-26T18:32:00+01:00Marcus MohrDiscussion: Handling Boundary Conditions in Vector-Valued FunctionsHi,
!697 triggered a discussion of how to handle boundary conditions for vector-valued functions, in this case for the `CSFVectorFunction` family.
The primary questions are:
1. Which of the following two options do we choose:
a. T...Hi,
!697 triggered a discussion of how to handle boundary conditions for vector-valued functions, in this case for the `CSFVectorFunction` family.
The primary questions are:
1. Which of the following two options do we choose:
a. The boundary condition for a `CSFVectorFunction` is a property of that function alone, i.e. we enforce that the individual component functions always have the same attribute.
b. We allow that the individual component functions carry different attributes.
2. If we go for the second option (b), this induces the question what the method `CSFVectorFunction::getBoundaryCondition()` should return. The suggestion was that in this case we return an array of size dimension with the boundary condition attributes of the scalar sub-functions.
**Note:** The current implementation is (blame on me) inconsistent in this respect. It started out assuming option (a), but does not enforce this and even offers a method to set the attribute on individual components.
In order to better understand the pros and cons of the approaches let us look at the *free-slip* boundary condition as an example:
```math
\mathbf{u} \cdot \mathbf{n} = 0
\enspace,\quad
\frac{\partial (\mathbf{u} \cdot \mathbf{t})}{\partial \mathbf{n} } = 0
```
with normal and tangential vectors $\mathbf{n}$ and $\mathbf{t}$.
**PRO for (b):**
- In the case of a cartesian box, the free-slip condition on a face simply translates into having a homogeneous Dirichlet boundary condition for the scalar sub-function of the normal component and homogeneous Neumann boundary conditions for the two other sub-functions. If we go for (b), this can nicely be used to handle the case.
- Always setting a well-defined boundary condition on each sub-functions is IMHO good. And, without the need to enforce uniqueness over the sub-functions, is simple in the current implementation.
**CONS for (b):**
- In any case where the domain boundaries are not axis-aligned, the splitting cannot be exploited. Especially in the case of an annulus or a thick spherical shell, we then need to set the boundary condition of the scalar sub-functions to `FreeslipBoundary`. Mathematically this feels a little unnatural for a scalar function.
- We would handle the `CSFVectorFunction` case differently from any other vector-valued FEM function space.
An *outflow* conditions of the form
```math
\frac{\partial (\mathbf{u} \cdot \mathbf{n})}{\partial \mathbf{n} } = 0
\enspace,\quad
\frac{\partial (\mathbf{u} \cdot \mathbf{t})}{\partial \mathbf{n} } = 0
```
could be treated in either scenario:
- in (a) we could set the attribute of the `CSFVectorFunction` to `NeumannBoundary`
- in (b) we could do the same for the scalar sub-functions
Any input on the question is welcome.
Cheers
MarcusMarcus MohrMarcus Mohrhttps://i10git.cs.fau.de/hyteg/hyteg-operators/-/issues/2How do we make sure, that a new generation does not damage an old one?2024-02-13T14:45:16+01:00Michael ZikeliHow do we make sure, that a new generation does not damage an old one?While generating operators in different precisions, I came by an example that is currently not doable with the generate script.
I either have to generate almost all things that I want to use at once, or some files are simply not known by...While generating operators in different precisions, I came by an example that is currently not doable with the generate script.
I either have to generate almost all things that I want to use at once, or some files are simply not known by the compiler or linker.
This is not only a problem of different precisions but, for the CMakeLists of different function spaces as well.
The example can be found here, [HFG !31](https://i10git.cs.fau.de/terraneo/hyteg-form-generator/-/issues/31#note_28113).https://i10git.cs.fau.de/hyteg/hyteg/-/issues/244Shadowing makes pipeline fail with FOG operators in app directory2024-02-13T16:33:09+01:00Marcus MohrShadowing makes pipeline fail with FOG operators in app directoryHi,
I am encountering pipeline failure due to shadowing variable declarations in case of a generated operator placed inside an app's directory:
```
../../../apps/nlDiffusion/P1ElementwiseNonlinearDiffusionNewtonGalerkin.cpp:971:26: erro...Hi,
I am encountering pipeline failure due to shadowing variable declarations in case of a generated operator placed inside an app's directory:
```
../../../apps/nlDiffusion/P1ElementwiseNonlinearDiffusionNewtonGalerkin.cpp:971:26: error: declaration shadows a local variable [-Werror,-Wshadow]
const real_t tmp_7_WHITE_UP = tmp_6_WHITE_UP * 0.050086823222829389;
^
```
See e.g. pipeline [#62046](https://i10git.cs.fau.de/hyteg/hyteg/-/pipelines/62046)Fabian BöhmFabian Böhmhttps://i10git.cs.fau.de/hyteg/hyteg/-/issues/243Bringing TerraNeo app to master2024-03-27T14:58:43+01:00Nils KohlBringing TerraNeo app to masterThis issue collects subtasks that need to be tackled to step-by-step migrate the TerraNeo to the master branch.
**State:**
As of writing this issue there are multiple branches with different states of the `3d_convection` app.
**Agenda...This issue collects subtasks that need to be tackled to step-by-step migrate the TerraNeo to the master branch.
**State:**
As of writing this issue there are multiple branches with different states of the `3d_convection` app.
**Agenda:**
Rebuild `3d_convection` app (maybe rename) only from features that are in the master branch. If some feature is missing - add it to master first. This is just a proposal, but possibly the most straightforward way to accomplish this.
For example:
* branch away from master
* copy _single_ feature from older branch
* create a small test
* merge it back to master
* continue
**Subtasks**
- [ ] separate RHS assembly from operator classes (@dascoli told me there are some old functions from @MarkusWiedemann that could be reused)
- [ ] sort out how/where to apply all the projections (free-slip and pressure) in the solvers @burkhart @dascoli
- [ ] bring solvers to master
- [ ] generate all required operators and have them in master
- [ ] use new generated operators and build composites where necessary
- [ ] modularization of initialization functions for temperature, viscosity, etc. (could be directory in `src/terraneo`)
- [ ] modularization of I/O functions (readers for file formats etc.)
- [ ] reduce/streamline number of FE function in TerraNeo app, maybe add `std::map< std::string, std::shared_ptr< FunctionType > >` (see [here]( https://i10git.cs.fau.de/hyteg/hyteg/-/blob/0fd62a325c68ac8da01661356e4812a8baaf3796/apps/3d_convection/ConvectionSimulationEntropy.hpp#L446))
Please feel free to contribute to this task including discussion: @IsabelPapanagnou @burkhart @dascoli @pponkumar @mohrIsabel PapanagnouIsabel Papanagnouhttps://i10git.cs.fau.de/hyteg/hyteg/-/issues/241Uniform way to apply all GeometryMap subclasses to the SetupStorage2024-01-24T11:24:55+01:00Nils KohlUniform way to apply all GeometryMap subclasses to the SetupStorageThe application of a GeometryMap to all primitives currently works differently for each type/subclass.
Most maps implement a static method `setMap( ... )` that does that. However, the signature is always different which makes the select...The application of a GeometryMap to all primitives currently works differently for each type/subclass.
Most maps implement a static method `setMap( ... )` that does that. However, the signature is always different which makes the selection of a map during run time annoying. Some examples:
* [IcosahedralShellMap::setMap()](https://i10git.cs.fau.de/hyteg/hyteg/-/blob/7d4e5176463d7fa5e7c3391fc1d56f30a2aa83ce/src/hyteg/geometry/IcosahedralShellMap.hpp#L239)
* [TokamakMap::setMap()](https://i10git.cs.fau.de/hyteg/hyteg/-/blob/7d4e5176463d7fa5e7c3391fc1d56f30a2aa83ce/src/hyteg/geometry/TokamakMap.hpp#L357)
* [CircularMap](https://i10git.cs.fau.de/hyteg/hyteg/-/blob/7d4e5176463d7fa5e7c3391fc1d56f30a2aa83ce/src/hyteg/geometry/CircularMap.hpp) does not even have that method.
The procedure then usually goes as follows:
```
SetupPrimitiveStorage setupStorage( meshInfo, numProcesses );
// Different call for all maps
AnnulusMap::setMap( setupStorage );
```
Would be nice if we could do something like this:
```
void createStorage( std::shared_ptr< GeometryMap > someMap )
{
SetupPrimitiveStorage setupStorage( meshInfo, numProcesses );
// Same call for all maps.
someMap->setMap( setupStorage );
}
```
Is there any reason `setMap()` cannot be a virtual method in the `GeometryMap` class? Maybe I am missing something.Nils KohlNils Kohlhttps://i10git.cs.fau.de/hyteg/hyteg/-/issues/240Non double precision operators2024-01-17T11:36:53+01:00Daniel BauerNon double precision operatorsRemaining issue from https://i10git.cs.fau.de/hyteg/hyteg/-/merge_requests/673#note_27144.
The operators are only generated for double precision.
In single precision builds, vectorization and float conversion warnings are disabled to ma...Remaining issue from https://i10git.cs.fau.de/hyteg/hyteg/-/merge_requests/673#note_27144.
The operators are only generated for double precision.
In single precision builds, vectorization and float conversion warnings are disabled to make the CI pass.
We have to decide whether we want specialized implementations for all precisions and if yes, integrate them.https://i10git.cs.fau.de/hyteg/hyteg/-/issues/239Reducing metadata for ADIOS2 checkpointing2024-01-19T19:01:37+01:00Ponsuganth Ilangovan Ponkumar IlangoReducing metadata for ADIOS2 checkpointingCurrently with ADIOS2, for checkpointing, we write out data for all primitives separately, for this we need to create adios local variables for each primitive types. As recently discussed, this growing number of local variables puts a he...Currently with ADIOS2, for checkpointing, we write out data for all primitives separately, for this we need to create adios local variables for each primitive types. As recently discussed, this growing number of local variables puts a heavy load on the ADIOS2 metadata management and thus affecting IO performance.
The current idea is to only write out data for the highest dimensional primitive (Cell in 3D, Face in 2D), and communicate to other primitive types when reading in at a checkpoint.Ponsuganth Ilangovan Ponkumar IlangoPonsuganth Ilangovan Ponkumar Ilangohttps://i10git.cs.fau.de/hyteg/hyteg/-/issues/236Numeric instability2023-11-20T15:15:43+01:00Daniel BauerNumeric instabilityIn !665 I introduced some very small changes.
As a result, the `P2P1StokesElementwiseAnnulusBlendingTest` fails but only when compiled with the icx compiler in single precision (e.g. https://i10git.cs.fau.de/hyteg/hyteg/-/jobs/1139667).
...In !665 I introduced some very small changes.
As a result, the `P2P1StokesElementwiseAnnulusBlendingTest` fails but only when compiled with the icx compiler in single precision (e.g. https://i10git.cs.fau.de/hyteg/hyteg/-/jobs/1139667).
Nevertheless, once it succeeded with the same code and configuration (https://i10git.cs.fau.de/hyteg/hyteg/-/jobs/1139562).
Other compilers do not seem to be affected by the changes at all.
After playing around with this compiler on one of our workstations, I have the suspicion that numerical instability is at play here.
The test solves the stokes system on a 2D annulus with GMG (Uzawa smoother) and also the HyTeG MinRes solver.
After solving with MinRes, the discrete L2 error in the pressure is reduced to 3.6e-2 in both single and double precision.
Let us assume that this is the discretization error.
The GMG reduces the error to almost the same number.
Only when compiling with icx and using single precision the error (GMG) is not reduced below 1.1e-1.
I can reproduce these exact numbers on one of our workstations.
The ci uses the compile flag `-fp-model=precise` to force the compiler to not reorder floating point operations (`WALBERLA_BUILD_WITH_FASTMATH=OFF`).
If I change the compiler options to either of `-O3 -xHost -fp-model=precise` or `-O3`, the test succeeds.
With `-O3 -fp-model=precise` the test fails.
I think that either of these combinations of options allows the compiler to insert FMA instructions.
As I learnt today, those are not only faster but also more accurate because the result is only rounded once, after performing both FP operations.
Without the explicit `-fp-model=precise`, the compiler uses
> [...] more aggressive optimizations when implementing floating-point calculations. These optimizations increase speed, but may affect the accuracy or reproducibility of floating-point computations.
Although I am not sure that this is actually the explanation, I am very concerned that convergence is destroyed by setting `-fp-model=precise`.
To me this suggests that our implementation is numerically unstable.https://i10git.cs.fau.de/hyteg/hyteg/-/issues/235Version Requirements for Optional Dependencies2024-03-28T10:12:56+01:00Marcus MohrVersion Requirements for Optional DependenciesHi,
!661 made me think whether we should list in our [README](https://i10git.cs.fau.de/hyteg/hyteg/-/blob/master/README.md) minimal versions for our optional dependencies and, if yes, what these would be?
| Dependency | Requirement |
|...Hi,
!661 made me think whether we should list in our [README](https://i10git.cs.fau.de/hyteg/hyteg/-/blob/master/README.md) minimal versions for our optional dependencies and, if yes, what these would be?
| Dependency | Requirement |
| :--------- | :---------- |
| MPI | which standard must be supported? |
| ADIOS2 | >= 2.9.2 (this includes the bugfix we need for current checkpointing implementation) |
| PETSc | no idea how far back we could go with PETSc |
| Trilinos | see !661 |
| ParMETIS | no idea |
| Doxygen | >= 1.10.0 (for best results) |
Any thoughts?
Cheers
MarcusMarcus MohrMarcus Mohrhttps://i10git.cs.fau.de/hyteg/hyteg/-/issues/234Rewrite N1E1toP1LiftingTest and N1E1RestrictionTest using Eigen sparse matrix...2023-10-27T15:21:36+02:00Daniel BauerRewrite N1E1toP1LiftingTest and N1E1RestrictionTest using Eigen sparse matrix proxyOnce !658 lands.Once !658 lands.Daniel BauerDaniel Bauerhttps://i10git.cs.fau.de/hyteg/hyteg/-/issues/232Generic surface integral implementation2023-10-31T10:49:53+01:00Ponsuganth Ilangovan Ponkumar IlangoGeneric surface integral implementationI believe we currently don't have the functionality to perform generic surface integrals on the domain, the idea is to mark the boundaries of a macro primitive and perform surface integrals if needed on them. The marker can be the existi...I believe we currently don't have the functionality to perform generic surface integrals on the domain, the idea is to mark the boundaries of a macro primitive and perform surface integrals if needed on them. The marker can be the existing `NeumannBoundary` markers or some more generic markers for surfaces inside the domain. The forms needs to be modified a bit to accommodate the cell coordinates and the facet coordinated where the surface integral needs to be evaluated.Ponsuganth Ilangovan Ponkumar IlangoPonsuganth Ilangovan Ponkumar Ilango