Use interp::linear instead of hand rolled interpolation

[moprak-nrel]

Summary

There are a few existing places where we have hand rolled interpolation code, this PR aims to identify and change those cases to use interp::linear and interp::bilinear where appropriate.

Pull request type

Please check the type of change introduced:

• Bugfix
• Feature
• Code style update (formatting, renaming)
• Refactoring (no functional changes, no api changes)
• Build related changes
• Documentation content changes
• Other (please describe):

Checklist

The following is included:

• new unit-test(s)
• new regression test(s)
• documentation for new capability

This PR was tested by running:

• the unit tests
  • on GPU
  • on CPU
• the regression tests
  • on GPU
  • on CPU

Additional background

Issue Number:

[moprak-nrel]

abl_godunov_mpl_amr is failing with large diffs on CPU, but not abl_godunov_mpl. Diffs:

            variable name            absolute error            relative error
                                        (||A - B||)         (||A - B||/||A||)
 ----------------------------------------------------------------------------
 level = 0
 density                                          0                         0
 gpx                                0.0009849055158             0.02256242477
 gpy                                 0.000992176435             0.01955978549
 gpz                                 0.004357075842               0.196405001
 mu_turb                             0.001029929218            0.001484382638
 p                                    0.06302397209             0.02391226761
 temperature                        0.0008945087092           2.883477209e-06
 temperature_mueff                   0.003090096663            0.001484372451
 velocityx                           0.003600159854           0.0003721130846
 velocityy                           0.003598752766           0.0003739956677
 velocityz                           0.006753192143            0.003581550337
 velocity_mueff                      0.001029929218            0.001484361244
 level = 1
 density                                          0                         0
 gpx                                 0.003049076587             0.04314988174
 gpy                                 0.003049700883             0.05081915065
 gpz                                  0.01039557679              0.2624832425
 mu_turb                            0.0006245218267            0.002057266101
 p                                    0.07482930219             0.03434654357
 temperature                         0.001680129268           5.564554645e-06
 temperature_mueff                   0.001873752855            0.002057233833
 velocityx                           0.005453110754           0.0005738780405
 velocityy                           0.005452998487           0.0005772748247
 velocityz                            0.01251170997            0.006587137699
 velocity_mueff                     0.0006245218267            0.002057198333

[moprak-nrel]

I'm fairly certain this diff is due to a bug in the current interpolation, this PR should fix it.

Here's a sample output where the old and new interpolations differ
ht: 276.0416667 theights[il]: 281.25 theights[ir]: 302.0833333 old_interp: 300.0001731 new_interp: 300.0005171

For the older interpolation to work correctly, theights[il] < ht <= theights[ir] needs to hold, which is clearly not the case as seen in the output. I suspect this is due to how the indicies are computed in the presence of refinement here:

amr-wind/amr-wind/equation_systems/icns/source_terms/ABLMeanBoussinesq.cpp

Lines 90 to 96 in baefdc9

	const amrex::Real ht = problo[idir] + (iv[idir] + 0.5) * dx[idir];
	const int il = amrex::min(k / lp1, nh_max);
	const int ir = il + 1;
	temp = tvals[il] +
	((tvals[ir] - tvals[il]) / (theights[ir] - theights[il])) *
	(ht - theights[il]);

[moprak-nrel]

TODOs, other places to check:

• equation_systems/icns/source_terms/HurricaneForcing.cpp
• equation_systems/icns/source_terms/BodyForce.cpp
• equation_systems/temperature/source_terms/BodyForce.cpp
• equation_systems/temperature/source_terms/HurricaneTempForcing.cpp

[moprak-nrel]

I'm fairly certain this diff is due to a bug in the current interpolation, this PR should fix it.

Here's a sample output where the old and new interpolations differ ht: 276.0416667 theights[il]: 281.25 theights[ir]: 302.0833333 old_interp: 300.0001731 new_interp: 300.0005171

For the older interpolation to work correctly, theights[il] < ht <= theights[ir] needs to hold, which is clearly not the case as seen in the output. I suspect this is due to how the indicies are computed in the presence of refinement here:

amr-wind/amr-wind/equation_systems/icns/source_terms/ABLMeanBoussinesq.cpp

Lines 90 to 96 in baefdc9

	const amrex::Real ht = problo[idir] + (iv[idir] + 0.5) * dx[idir];
	const int il = amrex::min(k / lp1, nh_max);
	const int ir = il + 1;
	temp = tvals[il] +
	((tvals[ir] - tvals[il]) / (theights[ir] - theights[il])) *
	(ht - theights[il]);

Further evidence that the new interpolation is doing the right thing:
ht: 234.375 old interp: [ 239.5833333, 260.4166667 ] new interp: [ 218.75, 239.5833333 ]

[marchdf]

Fantastic work finding and fixing these bugs!

[marchdf]

How do we know we caught all the hand-rolled interps?

[moprak-nrel]

How do we know we caught all the hand-rolled interps?

All of these seemed to have a similar code pattern of looking at nh_max = size() - 2 and k/lp1 where lp1 = lev +1. I grepped for variants of any of these, but can't say for sure that I haven't missed cases that name things differently.