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/202Add missing member functions to DGFunction and P0Function2023-06-29T08:21:04+02:00Marcus MohrAdd missing member functions to DGFunction and P0FunctionOur scalar function classes (like `P1Function` and `P2Function`) offer the three methods:
- getMaxValue()
- getMinValue()
- getMaxMagnitude()
Currently `DGFunction` implements only the latter. For consistency we should implement also th...Our scalar function classes (like `P1Function` and `P2Function`) offer the three methods:
- getMaxValue()
- getMinValue()
- getMaxMagnitude()
Currently `DGFunction` implements only the latter. For consistency we should implement also the former two and have `P0Function` delegate to them, like it already does for `getMaxMagnitude()`.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/186Rule-of-three for Functions2023-07-31T14:02:31+02:00Nils KohlRule-of-three for FunctionsFunctions only carry handles to access the allocated memory that resides on the Primitives.
In some sense the memory is owned by the Primitives. Currently destroying functions does not deallocate all of the data that was allocated before...Functions only carry handles to access the allocated memory that resides on the Primitives.
In some sense the memory is owned by the Primitives. Currently destroying functions does not deallocate all of the data that was allocated before.
This requires the implementation of custom destructors. So we should also implement copy-constructors and copy-assignment operators for all functions (https://en.cppreference.com/w/cpp/language/rule_of_three).
Copy constructors are relevant anyway (see #172 - although copy-constructors might not fix all problems related to solvers) and it is often misleading that Functions are not destroyed when leaving scope.
It is not clear to me whether introducing will break stuff - but on the other hand that would be a bad sign anyway.
What is missing?
* methods to remove data from Primitives
* destructors in all Function classes
* copy-constructors in all Function classes
* copy-assignment operators in all Function classesNils KohlNils Kohlhttps://i10git.cs.fau.de/hyteg/hyteg/-/issues/184Make evaluate() method of FE Function Classes Blending-Aware2023-05-26T16:53:56+02:00Marcus MohrMake evaluate() method of FE Function Classes Blending-AwareHi,
as recently discussed on [MatterMost](https://i10chat.cs.fau.de/lssall/pl/r1oucstaupdr7kucfpqf6144wr) the current implementation of the `evaluate()` methods of our FE function classes is implicitly assuming that the point coordinate...Hi,
as recently discussed on [MatterMost](https://i10chat.cs.fau.de/lssall/pl/r1oucstaupdr7kucfpqf6144wr) the current implementation of the `evaluate()` methods of our FE function classes is implicitly assuming that the point coordinates it receives are w.r.t. to the computational domain. If the coordinates are "blended", i.e. w.r.t. the physical domain, we will either get incorrect values for the point, or none at all, since no enclosing primitive will be found.
The outcome of our discussion IMHO is that we want to change this behaviour. Technically the requirement that the blending function must be globally a homeomorphism allows us to map the coordinates from the physical domain back to the computational domain, using a primitive's local blending map, and then check, whether the mapped point belongs to the primitive in the computational domain.
Note that this change will, of course, affect existing code that passes points with coordinates in the computational domain, such as in the current [MMOC implementation](https://i10git.cs.fau.de/hyteg/hyteg/-/blob/master/src/coupling_hyteg_convection_particles/MMOCTransport.hpp#L384).
As an aside I suggest to also remove in the process the specialisation of `evaluate()` for `real_t` using `constexpr if`, to get cleaner code.
Cheers
MarcusMarcus MohrMarcus Mohrhttps://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/151Add interpolate() method to CSFVectorClass that supports BoundaryUID2021-05-15T19:01:05+02:00Marcus MohrAdd interpolate() method to CSFVectorClass that supports BoundaryUIDIt would be convenient, if the CSFVectorClass was to provide an ```interpolate()``` method that allows passing a BoundaryUID. This should be straightforward as the **P[12]Function** classes and the underlying **[Vertex|Edge]DoFFunction**...It would be convenient, if the CSFVectorClass was to provide an ```interpolate()``` method that allows passing a BoundaryUID. This should be straightforward as the **P[12]Function** classes and the underlying **[Vertex|Edge]DoFFunction** classes already support that.Marcus MohrMarcus Mohrhttps://i10git.cs.fau.de/hyteg/hyteg/-/issues/146Add parent classes for smoothers2022-08-11T18:09:57+02:00Andreas WagnerAdd parent classes for smoothersIn our meeting today the question was posed, which methods we should include in the abstract Operator interface.
One idea was to include all the smoothing procedures to make it easier to switch between operators and smoothers at runtime...In our meeting today the question was posed, which methods we should include in the abstract Operator interface.
One idea was to include all the smoothing procedures to make it easier to switch between operators and smoothers at runtime. This is currently not possible and since all our operators have different smoothers, changing the operator often requires commenting in large amounts of code to pass compilation, even if the respective code paths are never taken.
This idea was rejected to avoid too much clutter in the base class in the long run and since from an idiomatic point of view a class should be able to implement all the methods of the parent class.
An alternative suggestion was to implement an interface for every smoother type. For the Jacobi, GS and symmetric GS this could look as follows:
```cpp
template < typename Function >
class JacobiSmoothable
{
public:
virtual ~JacobiSmoothable() = default;
virtual void smooth_jac( const Function& dst,
const Function& rhs,
const Function& tmp,
const real_t& relax,
size_t level,
DoFType flag ) const = 0;
};
template < typename Function >
class GSSmoothable
{
public:
virtual ~GSSmoothable() = default;
virtual void smooth_gs( const Function& dst, const Function& rhs, size_t level, DoFType flag ) const = 0;
};
template < typename Function >
class GSBackwardsSmoothable
{
public:
virtual ~GSBackwardsSmoothable() = default;
virtual void smooth_gs_backwards( const Function& dst, const Function& rhs, size_t level, DoFType flag ) const = 0;
};
// ... all the other smoothers ..
```
Operators, which want to support some of the smoothers can then inherit from the respective interfaces.
For a reduced version of the P1ConstantOperator which only support Jacobi and forward GS this could look like this:
```cpp
template < class P1Form, bool Diagonal = false, bool Lumped = false, bool InvertDiagonal = false >
class P1ConstantOperator : public Operator< P1Function< real_t >, P1Function< real_t > >,
public JacobiSmoothable< P1Function< real_t > >,
public GSSmoothable< P1Function< real_t > >
{
public:
void smooth_jac( const P1Function< real_t >& dst,
const P1Function< real_t >& rhs,
const P1Function< real_t >& tmp,
const real_t& relax,
size_t level,
DoFType flag ) const override
{
// implementation ...
}
void smooth_gs( const P1Function< real_t >& dst, const P1Function< real_t >& rhs, size_t level, DoFType flag ) const override
{
// implementation ...
}
};
```
Our solvers can then dynamically at runtime check if the operator supports the given smoother, and throw an exception if this is not the case:
```cpp
template < class OperatorType >
class GaussSeidelSmoother : public Solver< OperatorType >
{
public:
void solve( const OperatorType& A,
const typename OperatorType::srcType& x,
const typename OperatorType::dstType& b,
const size_t level ) override
{
if ( const auto* A_gs = dynamic_cast< const GSSmoothable< typename OperatorType::srcType >* >( &A ) )
{
A_gs->smooth_gs( x, b, level, flag_ );
}
else
{
throw std::runtime_error( "The Gauss-Seidel Operator requires the GSSmoothable interface." );
}
}
// ...
};
template < class OperatorType >
class SymmetricGaussSeidelSmoother : public Solver< OperatorType >
{
public:
void solve( const OperatorType& A,
const typename OperatorType::srcType& x,
const typename OperatorType::dstType& b,
const size_t level ) override
{
if ( const auto* A_gs = dynamic_cast< const GSSmoothable< typename OperatorType::srcType >* >( &A ) )
{
A_gs->smooth_gs( x, b, level, flag_ );
}
else
{
throw std::runtime_error( "The symmetric Gauss-Seidel Operator requires the GSSmoothable interface." );
}
if( const auto* A_gs_bw = dynamic_cast< const GSBackwardsSmoothable< typename OperatorType::srcType >* >( &A ))
{
A_gs_bw->smooth_gs_backwards( x, b, level, flag_ );
}
else
{
throw std::runtime_error( "The symmetric Gauss-Seidel Operator requires the GSBackwardsSmoothable interface." );
}
}
// ...
};
```
Thus in the example above the `GaussSeidelSmoother` would work for the `P1ConstantLaplaceOperator`, while the `SymmetricGaussSeidelSmoother` would throw an exception:
```
P1ConstantLaplaceOperator A( /* ... */ );
GaussSeidelSmoother< P1ConstantLaplaceOperator > GS ( /* ... */ );
SymmetricGaussSeidelSmoother< P1ConstantLaplaceOperator > SGS ( /* ... */ );
// works
GS.solve( A, x, b, level );
std::cout << "GS finished!" << std::endl;
// throws an exception at runtime
SGS.solve( A, x, b, level );
std::cout << "SGS finished!" << std::endl;
```Andreas WagnerAndreas Wagnerhttps://i10git.cs.fau.de/hyteg/hyteg/-/issues/142VectorFunction enhancements2021-04-01T15:54:36+02:00Marcus MohrVectorFunction enhancementsHi,
this issue is based on the comments for !381. That merge enhanced the implementation and functionality of the P1VectorFunction class. My comments were:
- Maybe we can move some of the methods (e.g. add, assign, dotLocal,...) upwar...Hi,
this issue is based on the comments for !381. That merge enhanced the implementation and functionality of the P1VectorFunction class. My comments were:
- Maybe we can move some of the methods (e.g. add, assign, dotLocal,...) upwards to a common base class instead of re-implementing basically the same stuff again in P2VectorFunction. But that is under discussion anyhow. Question is, if we get this to work nicely together with the templatisations.
- In that respect maybe also ```filter``` could be a free function?
- I definitely suggest to rename ```getMaxMagnitude()``` for VectorFunctions, as IMHO the name is very misleading. If I ask for the maximal magnitude of a velocity field e.g. I'd expect to get that back and not the maximal magnitude over all components.
Cheers
MarcusMarcus MohrMarcus Mohrhttps://i10git.cs.fau.de/hyteg/hyteg/-/issues/135Suggested changes to new VectorFunctions2020-09-25T17:30:26+02:00Marcus MohrSuggested changes to new VectorFunctionsHi,
I like the new ```P[12]VectorFunction``` classes. However, I'd like to suggest some changes:
1) I don't think that it is very flexible to have the scalar subfunctions be defined as follows:
```
P2Function< ValueType > u;
P2Function...Hi,
I like the new ```P[12]VectorFunction``` classes. However, I'd like to suggest some changes:
1) I don't think that it is very flexible to have the scalar subfunctions be defined as follows:
```
P2Function< ValueType > u;
P2Function< ValueType > v;
P2Function< ValueType > w;
```
For one the naming is closely tied to the current use with CFD velocity fields. Second you cannot access them in an automated fashion by indexing. Hence I suggest to change this to
```
std::array< P2Function< ValueType >, 3 > component;
P2Function< ValueType >& u = component[0];
P2Function< ValueType >& v = component[1];
P2Function< ValueType >& w = component[2];
```
The constructors then would have change also
```
P2VectorFunction( const std::string& _name,
const std::shared_ptr< PrimitiveStorage >& storage,
size_t minLevel,
size_t maxLevel )
: component( {{
{ _name + "_u", storage, minLevel, maxLevel },
{ _name + "_v", storage, minLevel, maxLevel },
{ storage->hasGlobalCells() ? P2Function< ValueType >( _name + "_w", storage, minLevel, maxLevel ) : P2Function< ValueType >( _name + "_w_dummy", storage ) }}} )
{}
```
2. I'd prefer to place the header files in ```src/hyteg/p[12]functionspace``` instead of ```src/hyteg/composites```. For me this feels more natural. Also the scalar ```P2Function``` is already a _composite_ of a VectorDoF and an EdgeDoFFunction.
3. Is there any reason why the VectorFunctions do not inherit from Function?
Cheers
MarcusAndreas WagnerAndreas Wagnerhttps://i10git.cs.fau.de/hyteg/hyteg/-/issues/125Refactor PETSc matrix assembly API2022-02-04T18:13:00+01:00Nils KohlRefactor PETSc matrix assembly APII suggest improving the `createMatrix()` interface to simplify adding new operators.
Overloading of `createMatrix()` could be one solution - the current templates are
confusing and hard to get right. This especially applies for the impl...I suggest improving the `createMatrix()` interface to simplify adding new operators.
Overloading of `createMatrix()` could be one solution - the current templates are
confusing and hard to get right. This especially applies for the implementation of
more complicated operators or wrappers such as the `UnsteadyDiffusionOperator`.https://i10git.cs.fau.de/hyteg/hyteg/-/issues/124Preconditioner interface2020-04-09T17:20:03+02:00Daniel DrzisgaPreconditioner interfaceCurrently, the preconditioner interface is very similar to the Solver interface. However, when calling `solve(A, x, b)` in a preconditioner, the `b` variable gets changed instead of the `x` variable. This makes it impossible to use a ge...Currently, the preconditioner interface is very similar to the Solver interface. However, when calling `solve(A, x, b)` in a preconditioner, the `b` variable gets changed instead of the `x` variable. This makes it impossible to use a general solver as a preconditioner. For instance, multgrid within a PCG solver.
My proposal: Swap the meaning of `x` and `b` in all preconditioners located in `solvers/preconditioners`. Make this swap consistent in solvers where preconditioners are currently supported, e.g., `CGSolver`, `MinResSolver`, and maybe others?https://i10git.cs.fau.de/hyteg/hyteg/-/issues/123Domain VTK output into single file via MPI2020-12-01T12:44:51+01:00Nils KohlDomain VTK output into single file via MPICurrently each process outputs an individual file when writing the domain partitioning.
This creates unnecessary many VTK files.Currently each process outputs an individual file when writing the domain partitioning.
This creates unnecessary many VTK files.https://i10git.cs.fau.de/hyteg/hyteg/-/issues/114Extend HyTeG to work on level 1 and 02020-02-24T16:04:02+01:00Nils KohlExtend HyTeG to work on level 1 and 0Apparently for some P1 GMG tests (even P1-P1 Stokes) we can plug in level 1 as minimum level and it works.
Level 0 could really be helpful to accelerate our coarse grid solvers.
* [x] indexing
* [x] functions
* [x] operator assembly
* [...Apparently for some P1 GMG tests (even P1-P1 Stokes) we can plug in level 1 as minimum level and it works.
Level 0 could really be helpful to accelerate our coarse grid solvers.
* [x] indexing
* [x] functions
* [x] operator assembly
* [x] single level solvers (CG, PETSc suite, ...)
* [x] grid transfer, multigrid solvers
* [x] fix 2D
* [x] function iterator, PETSc field split indexing
* [x] agglomeration
* [ ] ~~VTK~~ -> #116 Dominik Thoennesdominik.thoennes@fau.deDominik Thoennesdominik.thoennes@fau.dehttps://i10git.cs.fau.de/hyteg/hyteg/-/issues/111Flexible implementation of boundary handling for additive communication2020-02-24T16:01:01+01:00Nils KohlFlexible implementation of boundary handling for additive communicationDominik Thoennesdominik.thoennes@fau.deDominik Thoennesdominik.thoennes@fau.dehttps://i10git.cs.fau.de/hyteg/hyteg/-/issues/101Replace Pointnd with walberla::vector2022-05-30T10:32:16+02:00Dominik Thoennesdominik.thoennes@fau.deReplace Pointnd with walberla::vectorTo reduce code duplication we could possibly replace our own pointnd implementation with walberla::vector.
1. [ ] check whether they have the same behavior
2. [ ] do refactor (could be quite some work)To reduce code duplication we could possibly replace our own pointnd implementation with walberla::vector.
1. [ ] check whether they have the same behavior
2. [ ] do refactor (could be quite some work)