REMORA
Regional Modeling of Oceans Refined Adaptively
Loading...
Searching...
No Matches
REMORA_Math.H
Go to the documentation of this file.
1#ifndef REMORA_Math_H
2#define REMORA_Math_H
3#include <AMReX_REAL.H>
4
5/** \brief interpolate 1D arrays
6 *
7 * Interpolates 1D array beta at the location alpha_interp in the array alpha
8 * requiring 1D arrays alpha and beta to be the same size alpha_size.
9 */
10AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
11amrex::Real interpolate_1d(const amrex::Real* alpha, const amrex::Real* beta,
12 const amrex::Real alpha_interp, const int alpha_size)
13{
14
15 amrex::Real beta_interp;
16
17 //in case the interpolation point already exists in the array
18 //just return it
19 for (int i = 0; i < alpha_size; ++i) {
20 if (alpha[i] == alpha_interp) {
21 beta_interp = beta[i];
22 return beta_interp;
23 }
24 }
25
26 // we didn't find an exact match for alpha_interp in alpha,
27 // so we need linear interpolation/extrapolation
28 const amrex::Real* alpha_begin = &alpha[0];
29 const amrex::Real* alpha_end = &alpha[alpha_size];
30 amrex::Real max = *std::max_element(alpha_begin, alpha_end);
31 amrex::Real min = *std::min_element(alpha_begin, alpha_end);
32 if (alpha_interp >= min && alpha_interp <= max) //interpolate
33 {
34 for (int i = 0; i < alpha_size; ++i)
35 {
36 if (alpha_interp >= alpha[i] && alpha_interp <= alpha[i + 1])
37 {
38 //y = y0 + (y1-y0)*(x-x0)/(x1-x0);
39 amrex::Real y0 = beta[i];
40 amrex::Real y1 = beta[i + 1];
41 amrex::Real x = alpha_interp;
42 amrex::Real x0 = alpha[i];
43 amrex::Real x1 = alpha[i + 1];
44 beta_interp = y0 + (y1 - y0)*(x - x0) / (x1 - x0);
45 }
46 }
47 }
48 else //extrapolate
49 {
50 //y = y0 + ((x - x0) / (x1 - x0)) * (y1 - y0)
51 int i;
52 if (alpha_interp >= *alpha_end)
53 {
54 i = alpha_end - alpha_begin - 1;
55 }
56 else
57 {
58 i = 0;
59 }
60 amrex::Real y0 = beta[i];
61 amrex::Real y1 = beta[i + 1];
62 amrex::Real x = alpha_interp;
63 amrex::Real x0 = alpha[i];
64 amrex::Real x1 = alpha[i + 1];
65 beta_interp = y0 + ((x - x0) / (x1 - x0)) * (y1 - y0);
66 }
67
68 return beta_interp;
69
70using namespace amrex;
71}
72#endif
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real interpolate_1d(const amrex::Real *alpha, const amrex::Real *beta, const amrex::Real alpha_interp, const int alpha_size)
interpolate 1D arrays
Definition REMORA_Math.H:11