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/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/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/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/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 Ilangohttps://i10git.cs.fau.de/hyteg/hyteg/-/issues/231ADIOS2 output for unresolved particles.2023-10-19T12:10:36+02:00Nils KohlADIOS2 output for unresolved particles.(See also discussion in https://i10git.cs.fau.de/hyteg/hyteg/-/merge_requests/653#note_25667.)
ADIOS2 output is currently not supported for the unresolved particle feature. Either this has to be implemented as part of MESA-PD or through...(See also discussion in https://i10git.cs.fau.de/hyteg/hyteg/-/merge_requests/653#note_25667.)
ADIOS2 output is currently not supported for the unresolved particle feature. Either this has to be implemented as part of MESA-PD or through the coupling on the HyTeG side. I suppose the latter is simpler, but not sure yet.https://i10git.cs.fau.de/hyteg/hyteg/-/issues/229VTK output is not possible for functions in single precision2024-01-16T13:40:33+01:00Michael ZikeliVTK output is not possible for functions in single precisionWhile writing VTK output, functions that have another type than double precision are simply skipped and only the mesh structure and the functions of double or integer type are printed. There is no user notification, that this is done, th...While writing VTK output, functions that have another type than double precision are simply skipped and only the mesh structure and the functions of double or integer type are printed. There is no user notification, that this is done, the function values are just missing in the output file.
A possible workaround is to use the `.copyfrom()` member function to create a printable copied function object of value type double, however this works only for P1 functions.
Allowing VTK output for any precision (e.g. single, half, ...) would be a nice to have.Michael ZikeliMichael Zikelihttps://i10git.cs.fau.de/hyteg/hyteg/-/issues/212Incorporation of weak Dirichlet boundary conditions in test cases2024-02-20T09:55:24+01:00Fabian BöhmIncorporation of weak Dirichlet boundary conditions in test casesCurrently, the application of weak boundary conditions from the user view in a test case is complicated:
- in HyTeG solvers, the DoFtype has to be set to 'All' via solver.setDoFType( All ); because the DoFs on the boundary are now par...Currently, the application of weak boundary conditions from the user view in a test case is complicated:
- in HyTeG solvers, the DoFtype has to be set to 'All' via solver.setDoFType( All ); because the DoFs on the boundary are now part of the solution
- the standard strong boundary application has to be switched off in PETSc solvers by solver.disableApplicationBC( true );
- the user can/has to set the penalty depending on his application and operator (ensure coercivity of the operator: e.g. larger constant required in 3D)
Two ways to avoid giving the responsibility to do these small yet crucial operations to the user are:
- Inspecting the given operator within the solver and set DoFs and boundary handling according to the operator type (if it applies its BCs weakly)
(con: it is not the solver's responsibility to check how BCs should be applied but the operator's, the solver should 'just solve')
- introducing another boundary flag WeakDirichletBC and handling the case similarly to the existing free slip BC.
(con: additional functions for the new boundary flag would have to be implemented, and the boundary has to be set with the new flag at the start of the test case
pro: IMO intuitive way to handle weak BCs, no modifications in the solvers)
Are there other/better ideas to handle this?https://i10git.cs.fau.de/hyteg/hyteg/-/issues/171DG implementation2022-04-11T16:48:50+02:00Nils KohlDG implementation# Planned features
* Discontinuous Galerkin implementation with arbitrary choice of (polynomial) shape functions
* arbitrary order
* p-refinement over macro-elements
* matrix-free operator evaluation
* weak form generation via HFG
* eva...# Planned features
* Discontinuous Galerkin implementation with arbitrary choice of (polynomial) shape functions
* arbitrary order
* p-refinement over macro-elements
* matrix-free operator evaluation
* weak form generation via HFG
* evaluation of functions and linear functionals for right-hand side assembly generated via HFG
* direct volume to volume communication (interface primitives not used)
* SIP for Poisson
Probably also(?)
* hp-multigrid (maybe also projection into continuous element space?)
* weakly imposed boundary conditions
* upwind operator
# Roadmap
- [x] Volume DoF data structure + DGFunction class w/ basic indexing functionality
- [x] DGOperator class, apply, toMatrix()
- [x] mass form, interpolation
- [x] linear form for RHS
- [x] Laplace form (SIP, interface integrals)
- [x] communication (PackInfo etc.), macro-interface integrals
- [x] inhom. Dirichlet BCs.
- [ ] higher order shape functions
- [ ] HFG refactoring and fixed DGForm API
- [x] 3D extensions
- [ ] p-adaptivity
# Possible optimizations / open refactoring / dev notes
- [ ] optimize evaluation by precomputing and storing trafo matrices on each macro-face / macro-cell (it's recomputed every time in the evaluation function)
- [ ] make all methods in DGBasisInfo const
- [ ] give DGForms the following virtual methods
* `isConstant()`: true if the form is constant over the macro - so optimization can be applied
* `isInlined()`: true if entire kernel is available
* `onlyInnerIntegrals()`: true if only volume integrals are non-zero (think mass matrix) - this allows to implement efficient(?) local inversion of the block-diagonal system without communication
- [ ] probably we should not double the information about the local polynomial degree in each macro in the DGFunction AND VolumeDoFFunction
- [x] optimize VTK evaluation by new DGFunction::evaluate( faceID, microIdx, faceType, ... ) method that skips looking for the micro-element - we know already which one we are targeting (_this is not only an optimization but also fixes slightly incorrect visualization_)
- [ ] what to do with the DoFType 'flags'?
- [x] implement affine element-local evaluation directly in HFG
- [ ] Optimize HFG trafo matrix computations by letting it know which element vertices belong to the interface.
This should allow for some cancellations since some of the input points are actually equal.
- [ ] There is an issue if the boundary macros are marked as 'Inner' as done for Stokes usually - since then the DGOperator 'thinks' that there must be another volume macro on the other side. Boom.
- [ ] Instead of resizing the element matrices, the HFG should rather assert that the size is already correct.Nils KohlNils Kohlhttps://i10git.cs.fau.de/hyteg/hyteg/-/issues/134FunctionIterator< P2Function > not working correctly2023-05-24T15:56:57+02:00Marcus MohrFunctionIterator< P2Function > not working correctlyHi,
there seems (in my understanding) to be a problem with the implementation of the ````FunctionIterator```` class, when used on a ````P2Function````. When I run this test:
````
P2Function< real_t > p2Func( "P2 test func", storage, lev...Hi,
there seems (in my understanding) to be a problem with the implementation of the ````FunctionIterator```` class, when used on a ````P2Function````. When I run this test:
````
P2Function< real_t > p2Func( "P2 test func", storage, level, level );
for ( auto val : FunctionIterator< P2Function< real_t > >( p2Func, level ) )
{
if ( !val.isEdgeDoF() && !val.isVertexDoF() )
{
WALBERLA_ABORT( "P2 DoF must be either on micro edge or micro vertex!" );
}
}
````
it aborts, see [affa19070].
Cheers
MarcusNils KohlNils Kohlhttps://i10git.cs.fau.de/hyteg/hyteg/-/issues/129Blending PSPG inv diag operator2021-12-07T14:38:49+01:00Benjamin MannBlending PSPG inv diag operatorIn `P2P1ElementwiseBlendingStokesOperator` the constant `P1PSPGInvDiagOperator` is used.
Corresponding blending operator required?In `P2P1ElementwiseBlendingStokesOperator` the constant `P1PSPGInvDiagOperator` is used.
Corresponding blending operator required?Marcus MohrMarcus Mohrhttps://i10git.cs.fau.de/hyteg/hyteg/-/issues/121Dynamic function allocation on arbitrary levels and primitives2020-07-31T11:31:48+02:00Nils KohlDynamic function allocation on arbitrary levels and primitivesFor adaptive simulations we require a more flexible and dynamic mechanism to allocate function memory.
This especially involves arbitrary levels and primitives.
Potential issues:
- ~~min / maxLevel currently fixed~~
- internal functions...For adaptive simulations we require a more flexible and dynamic mechanism to allocate function memory.
This especially involves arbitrary levels and primitives.
Potential issues:
- ~~min / maxLevel currently fixed~~
- internal functions (e.g. in solvers) must be adapted as well (this depends strongly on the application)
- operators and communicators need to be adaptedAndreas WagnerAndreas Wagnerhttps://i10git.cs.fau.de/hyteg/hyteg/-/issues/96Black-Box Matrix Output2020-10-07T14:04:44+02:00Nils KohlBlack-Box Matrix OutputFor debugging purposes it would be nice to have some sort of black-box matrix output (to PETSc / Matlab) for arbitrary matrices - especially if they are assembled on-the-fly or are built with inverse matrices (e.g. Schur-complement).
Th...For debugging purposes it would be nice to have some sort of black-box matrix output (to PETSc / Matlab) for arbitrary matrices - especially if they are assembled on-the-fly or are built with inverse matrices (e.g. Schur-complement).
Therefore we could repeatedly apply the matrix to unit vectors to emulate a multiplication with the identity.
The resulting vectors can then be collected and written into multiple files or (probably better) be assembled into
a PETSc data structure (and then written using the available PETSc -> Matlab functions).Daniel DrzisgaDaniel Drzisgahttps://i10git.cs.fau.de/hyteg/hyteg/-/issues/88Fix/test SchurCG2020-07-31T11:23:44+02:00Nils KohlFix/test SchurCGMarcus MohrMarcus Mohr