feat: Add Taper Layers in wave solvers

[acitrain]

This PR add a first version of Taper layers for wave solvers
These layers are set to improve the efficiency of absoption of boundaries.
This is how the version version works:

• We compute a taper profile which is equal to one when you are outside the taper layers and to :
  $d_x = \displaystyle \frac{-3V_{max}}{2L}log(R)(\frac{x}{L})^2$ where $V_{max}$ is the maximum P-wavespeed, L the length of the taper layer and R the reflectivitity coeff (which is between $10^{-3}$ and $10^{-6}$
• Then we multiply by the array at time n and n+1 by the taperCoeff

For this first version only second order solver contains the taper functionality and the computation of the profile as been simplified in particular in corner cases.

Those two functionalities will be add in a next PR.

[acitrain]

Hi @rrsettgast

Could you take a look at it and approve if it seems good to you?
It is ready except integrated test but I want to do it after a code owner review

Thank you !

PS: I don't know why the docs build does not pass, it is only saying "We encountered a problem with a command while building your project. To resolve this error, double check your project configuration and installed dependencies are correct and have not changed recently."

[acitrain]

@sframba @rrsettgast @CusiniM can you review this please ?

Thanks

[Victor-M-Gomes]

This PR brings an important addition to the SEM wavepropagation module of GEOS. With this new absorbing boundary I've checked that one can achieve significantly better modelling results.

I suggest some small modifications to avoid some problems encountered during testing.

[Victor-M-Gomes]

Hello,

The current definition is not optimal.

For example: When m_timestepStabilityLimit==1 a correct dt (m_timeStep/dtOut) value is calculated, however, by setting m_timestepStabilityLimit==0 obliges all other calls to initializePostInitialConditionsPreSubGroups() to use the user-input dt (either from xml or pygeosx). This scenario is often present when interfacing GEOS through python (pygeosx).

One option would be to never set m_timestepStabilityLimit==0, however that implicates a dt calculation at each new call to initializePostInitialConditionsPreSubGroups(). This may be useful for certain applications, but one often just wants to calculate the dt once. Besides that, at the current state, if initializePostInitialConditionsPreSubGroups() is called a second time to calculate the dt (meaning that m_timestepStabilityLimit==1 is always true), the dt calculation does not work, giving nan values.

After exploring a bit, the nan values are due to a de-synchronization between GPU and CPU copies in the computeTimeStep function. By calling initializePostInitialConditionsPreSubGroups(), and therefore computeTimeStep, a second time, the ux_n, uy_n, and uz_n arrays are not updated by the loop in lines 553-558, thus the division operations in lines 569-574 will give nan.

A simple workaround is to avoid:

for( localIndex a = 0; a < sizeNode; ++a )
    {
      ux_n[a] = (real64)rand()/(real64) RAND_MAX;
      uy_n[a] = (real64)rand()/(real64) RAND_MAX;
      uz_n[a] = (real64)rand()/(real64) RAND_MAX;
}

simply writing:

for( localIndex a = 0; a < sizeNode; ++a )
    {
      ux_n[a] = 1.0;
      uy_n[a] = 1.0;
      uz_n[a] = 1.0;
}

and add the following loop to ensure the GPU values are not zero:

forAll< EXEC_POLICY >( sizeNode, [=] GEOS_HOST_DEVICE ( localIndex const a )
    {
      ux_n[a] = 1.0;
      uy_n[a] = 1.0;
      uz_n[a] = 1.0;
} );

Finally, it is not clear to me why we would need: return m_timeStep * m_cflFactor; in line 678.

Notice that these comments are valid to the acoustic case also.

Coming back to the origin of the problem, it may be wise to add a third state m_timestepStabilityLimit==2 (for example). When m_timestepStabilityLimit==1 the dt is always calculated, if m_timestepStabilityLimit==2 the dt is calculated once and then set to state 3 (one that does not exist, implying that no further modification of m_timeStep will take place), and if m_timestepStabilityLimit==0 the user-input dt is used (xml or pygeosx).

[Victor-M-Gomes]

Great job with the taper, choosing vMin over vMax was a really good choice since it usually leads to better results. This needs to be modified in the description of the PR though:
$d_x = \displaystyle \frac{-3V_{max}}{2L}log(R)(\frac{x}{L})^2$ becomes $d_x = \exp\left( \displaystyle \frac{3V_{min}}{2L}log(R)(\frac{L-x}{L})^2 dt \right)$ (the dt also appears in the definition).

After much testing the border reflections have been successfully attenuated and modelling results are significantly improved.

For the future (in another PR maybe) I would suggest substituting $V_{min} \times dt$ by some $\lambda_{min} = V_{min}/f_{max}$, where $f_{max}$ is the maximum frequency that one wishes to consider in the modelling. Right now $f_{max} == f_{sampling}$, the maximum frequency that exists in the data. By allowing the user to input $\lambda_{min}$ or $f_{max}$, the equation would have more physical meaning and the use of $dt$ could be avoided. Another option would be to set $f_{max} \simeq f_{Ricker} \times 3$ (or a similar factor) to avoid defining new inputs. The goal of the suggestion is to avoid using dt in the definition of the taper, while also inserting additional physical meaning.

Finally, it may be important to mention that the choice of L and R are somewhat empirical. For a given R, the greater the L, the better the attenuation, but that also comes with the cost of increasing the model to fit the taper. Thus, in general a compromise between L and R need to be found. For the bounds of R, I would (based on some tests) place them in the interval [ $10^{-3}$ , 1).