hyteg issueshttps://i10git.cs.fau.de/hyteg/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/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/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/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/195Sort out usage of namespaces2023-10-06T18:35:43+02:00Marcus MohrSort out usage of namespacesHi,
since I plan to close #158 now, I am adding here an issue that Nils, I guess, brought up then:
> ccd5f369 inspired me in adding a row here. Namespaces are all over the place, we may want to have that sorted out (i.e. only single na...Hi,
since I plan to close #158 now, I am adding here an issue that Nils, I guess, brought up then:
> ccd5f369 inspired me in adding a row here. Namespaces are all over the place, we may want to have that sorted out (i.e. only single namespace hyteg:: or add some more structure by directories, or something else ...). Our doxygen documentation and IDE code completion would also really benefit.
Cheers
MarcusNils KohlNils Kohlhttps://i10git.cs.fau.de/hyteg/hyteg/-/issues/178Shared forms for DG and CG2022-04-11T16:03:04+02:00Nils KohlShared forms for DG and CGDG forms are currently carrying multiple methods for the integration of
* volumes
* inner facets (coupling of element with itself)
* outer facets (coupling of element with neighbor element)
* domain boundary facets
* integration for fac...DG forms are currently carrying multiple methods for the integration of
* volumes
* inner facets (coupling of element with itself)
* outer facets (coupling of element with neighbor element)
* domain boundary facets
* integration for facet contributions required for Dirichlet boundaries
The volume integrals are **identical** for conforming elements if the same basis functions are chosen (which is currently always the case as we are using Lagrangian polynomials for both discretizations).
Thus it may be convenient to either
* wrap the forms for conforming discretizations as a "special case" of the non-conforming version
* replace the old interface altogether
Something similar to a wrapper is already available for DG forms that only involve volume integrals (e.g. mass):
https://i10git.cs.fau.de/hyteg/hyteg/-/blob/df450dd9fbfdcb98396ae62af6c1e7194b5e489d/src/hyteg/dgfunctionspace/DGFormVolume.hpp
Note that the DGOperator uses the virtual method `onlyVolumeIntegrals()` to check if facet integrals can be skipped, see:
https://i10git.cs.fau.de/hyteg/hyteg/-/blob/df450dd9fbfdcb98396ae62af6c1e7194b5e489d/src/hyteg/dgfunctionspace/DGFormVolume.hpp#L39
So even the DGOperator already checks if facet integration is necessary during run time.
This refactoring could unify the form interface and remove the "special" DG forms.
At the same time we could get rid of some of the `Point3D` related code.https://i10git.cs.fau.de/hyteg/hyteg/-/issues/172Function copy mechanism2024-02-16T10:39:55+01:00Nils KohlFunction copy mechanism# Problem
Especially in solvers (but probably also elsewhere) it is necessary to create copies of existing functions.
Consider for example the temporary variables in solvers. Those should have the same properties as the
function that r...# Problem
Especially in solvers (but probably also elsewhere) it is necessary to create copies of existing functions.
Consider for example the temporary variables in solvers. Those should have the same properties as the
function that represents the solution variable (e.g. boundary conditions).
Somehow such properties need to be copied/duplicated - but the function type might not be known when the copy
is requested.
A simple copy-constructor is probably not what we want either, since the duplicate may have another name, different
allocated min and max-level etc. So certain properties should not be copied.
There are two approaches: ctor or a method that offers a duplicate instance.
An orthogonal issue is the following: consider current solver implementation:
```c++
// pseudo solver
SomeSolver(
const std::shared_ptr< PrimitiveStorage >& storage,
uint_t minLevel,
uint_t maxLevel,
uint_t maxIter = std::numeric_limits< uint_t >::max(),
real_t tolerance = 1e-16,
)
// When the temporary functions are constructed, there is no access yet to an actual function.
// So the properties have to be copied later on - OR we make member functions pointers that are initialized in a second step.
: tmp_( "tmp", storage, minLevel, maxLevel )
// ...
void solve( const OperatorType& A, const FunctionType& x, const FunctionType& b, const uint_t level ) override
{
// Only at this point we even _know_ the required properties of the tmp variable.
tmp_.copyBoundaryConditionFromFunction( x );
// ...
```
# Proposal
a) ctor - cannot be virtual so we need to force implementation in derived functions via virtual function anyway(?)
```c++
// Function
Function( const FunctionType & other, std::string name, uint_t minLevel, uint_t maxLevel )
{
ctor_impl( other, name, minLevel, maxLevel );
}
virtual void Function::ctor_impl( ... )
{
WALBERLA_ABORT( "Not implemented." )
}
// Derived
void ctor_impl( ... )
{
// here we go
}
// usage
// pseudo solver
SomeSolver(
const std::shared_ptr< PrimitiveStorage >& storage,
uint_t minLevel,
uint_t maxLevel,
uint_t maxIter = std::numeric_limits< uint_t >::max(),
real_t tolerance = 1e-16,
)
// ...
std::shared_ptr< FunctionType > tmp_;
void solve( const OperatorType& A, const FunctionType& x, const FunctionType& b, const uint_t level ) override
{
// not even sure if this works
tmp_ = std::make_shared< Function >( x, "tmp", 2, 2 );
// ...
```
b) virtual copy/duplicate methods
```c++
// Function
virtual std::shared_ptr< FunctionType > duplicate( std::string name, uint_t minLevel, uint_t maxLevel, bool copyValues )
{
WALBERLA_ABORT( "Not implemented." )
}
// Derived
virtual std::shared_ptr< FunctionType > duplicate( std::string name, uint_t minLevel, uint_t maxLevel, bool copyValues )
{
// copy stuff
}
// usage
// pseudo solver
SomeSolver(
const std::shared_ptr< PrimitiveStorage >& storage,
uint_t minLevel,
uint_t maxLevel,
uint_t maxIter = std::numeric_limits< uint_t >::max(),
real_t tolerance = 1e-16,
)
// ...
std::shared_ptr< FunctionType > tmp_;
void solve( const OperatorType& A, const FunctionType& x, const FunctionType& b, const uint_t level ) override
{
// not even sure if this works
tmp_ = x.duplicate( "tmp", 2, 2, false );
// ...
```
Writing this I feel that the duplicate method approach is simpler and straightforward.
Input appreciated!Nils KohlNils Kohlhttps://i10git.cs.fau.de/hyteg/hyteg/-/issues/136Inconsistent use of PrimitiveStorage inside VTKOutput class2020-09-28T09:35:31+02:00Marcus MohrInconsistent use of PrimitiveStorage inside VTKOutput classHi,
at some point we should make our use of the PrimitiveStorage inside VTKOutput more consistent. Currently we pass a PrimitiveStorage to the object in its constructor. However, internally that get's used sometimes, but e.g. not in wri...Hi,
at some point we should make our use of the PrimitiveStorage inside VTKOutput more consistent. Currently we pass a PrimitiveStorage to the object in its constructor. However, internally that get's used sometimes, but e.g. not in writeP1() or writeP2(). The latter ask the **first** of the functions to export for its PrimitiveStorage. This can of course lead to consistency issues, see e.g. this snippet from ```EdgeDoFFunctionTest``` [cff70c1ce6c7b1cbaf78bf722774e61a661d586f]
```
auto p1 = std::make_shared< P1Function< real_t > >( "p1", storage2, minLevel, maxLevel );
VTKOutput vtkOutput( "../../output", "interpolate_test", storage );
vtkOutput.add( *p1 );
```
Also we implicitely assume that all functions that we add in order to write to the same VTK file use the same PrimitiveStorage.
Maybe we should make ```add()``` check whether the functions have the same storage as the one passed in the constructor?
Cheers
MarcusMarcus MohrMarcus Mohrhttps://i10git.cs.fau.de/hyteg/hyteg/-/issues/115Remove old stencil assembly code for P1 in 3D2021-07-27T17:52:24+02:00Nils KohlRemove old stencil assembly code for P1 in 3DThe functions that assemble the 3D stencils for P1 elements
that still operate directly on the face, edge and vertex primitives
and return stencils in linearized arrays should be removed and
the respective routines using these stencils s...The functions that assemble the 3D stencils for P1 elements
that still operate directly on the face, edge and vertex primitives
and return stencils in linearized arrays should be removed and
the respective routines using these stencils should be replaced
with the std::map variants.Nils KohlNils Kohlhttps://i10git.cs.fau.de/hyteg/hyteg/-/issues/84Align 2D and 3D stencil assembly routines2023-03-21T11:11:01+01:00Nils KohlAlign 2D and 3D stencil assembly routinesNils KohlNils Kohl