hyteg issueshttps://i10git.cs.fau.de/hyteg/hyteg/-/issues2024-03-08T08:26:10+01:00https://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/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/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/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/201Update GMG tutorial2023-02-17T14:32:07+01:00Nils KohlUpdate GMG tutorialThe GMG tutorial (05) is a little outdated and could use some more in-depth description to provide a nice starting point for new users.The GMG tutorial (05) is a little outdated and could use some more in-depth description to provide a nice starting point for new users.Nils 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/84Align 2D and 3D stencil assembly routines2023-03-21T11:11:01+01:00Nils KohlAlign 2D and 3D stencil assembly routinesNils KohlNils Kohl