adaptive-cfd / WABBIT

Wavelet Adaptive Block-Based solver for Interactions with Turbulence
https://www.cfd.tu-berlin.de/
GNU General Public License v3.0
56 stars 27 forks source link

reformulated neighboring codes, introducing full wavelet transformation and image denoising #64

Closed Arcadia197 closed 2 months ago

Arcadia197 commented 2 months ago

Reformulated neighboring codes

Previously, the code was only able to handle one neighbor for each ghost layer patch. However, for the full wavelet transformation several neighbors on different levels for the same patch will be needed. A new neighbor formulation enables the coexistence of neighbors on different levels. For all 56 possible patches in 3D we have now 3x56 possible neighbor relations. 2D and 3D codes have been merged to avoid confusion of neighbor codes. Taken from the code, the neighbor formulations are now constructed the following way:

!> neighbor codes: \n
!  ---------------
!>   1- 56 : lvl_diff =  0  (same level)
!>  57-112 : lvl_diff = +1  (coarser neighbor)
!> 113-168 : lvl_diff = -1  (finer   neighbor)
!> For each range, the different 56 entries are:
!> 01-08 : X side (4-,4+)
!> 09-16 : Y-side (4-,4+)
!> 17-24 : Z-side (4-,4+)
!> 25-32 : X-Y edge (2--, 2+-, 2-+, 2++)
!> 33-40 : X-Z edge (2--, 2+-, 2-+, 2++)
!> 41-48 : Y-Z edge (2--, 2+-, 2-+, 2++)
!> 49-56 : corners (---, +--, -+-, ++-, --+, +-+, -++, +++)

As one can see, sides have 4 possible codes for the 4 possible configurations of blocks on finer level, edges have 2 and corners 1 possible configurations. 2D codes onlz have X- and Y-side with 2 configurations and the X-Y edges with 1 configuration each.

With the change in code formulation the concerning code to compute the boundaries has been revisited and greatly simplified. It should hopefully be clear enough now.

Full wavelet transformation - adapt_tree_cvs

The full wavelet transformation on all levels has been implemented. It works by doing the following steps:

  1. Initialize the full grid from the leaf-layer down to level JMin. New blocks are marked as empty with the refinement flag.
  2. Do the full wavelet decomposition. We work downwards from the leaf-layer. Along the way, interior blocks (blocks which are neither leaf-blocks nor root-blocks on level JMin) can only be decomposed if they have correctly synched blocks from all neighbors on the same level. This way we avoid the need of any coarse extension on blocks besides the leaf-layer blocks. It also effectively creates an iterative decomposition algorithm that starts on the leaf-layer and for later steps mostly works level-wise down to JMin due to restrictions of blocks waiting for their medium-lvl neighbors. Luckily, as the grid information is already present, only synching block values and lgt-data is needed for each step.
  3. Adaption of grid is done only once after the full wavelet transformation was completed. With this, mask and norm information is only computed once. All unsignificant blocks can then be deleted
  4. Full wavelet reconstruction is applied, we work level-wise from lowest layer to the leaf-layer and update the SC of all daughters along the way, afterwards all non-leaf blocks are deleted as they are currently not used for further steps

This algorithm has the following benefits:

Image denoising

The iterative algorithm by Azzalini was implemented in order to compute image denoising. This was introduced as a new step after full wavelet decomposition and before the grid adaption. Possible, CVS iterative algorithm will be introduced here as well but is not yet tested. A new post-processing tool (wabbit-post --denoise) deploys the denoising algorithm on wabbit files of uniform grid to be denoised and a test was added to check the results. Tests show that the denoising algorithm works quite well for the bi-orthogonal wavelets.

Miscellaneous changes