diff --git a/.github/workflows/pypi-release.yml b/.github/workflows/pypi-release.yml index ad3c789f6e..cb016c524c 100644 --- a/.github/workflows/pypi-release.yml +++ b/.github/workflows/pypi-release.yml @@ -104,7 +104,7 @@ jobs: activate-environment: parcels python-version: "3.10" channels: conda-forge - - run: conda install -c conda-forge c-compiler pip + - run: conda install -c conda-forge pip - run: pip install parcels --no-cache - run: curl https://raw.githubusercontent.com/OceanParcels/parcels/main/docs/examples/example_peninsula.py > example_peninsula.py - run: python example_peninsula.py diff --git a/docs/conf.py b/docs/conf.py index 9af6eb0457..8d60edd01a 100755 --- a/docs/conf.py +++ b/docs/conf.py @@ -372,7 +372,6 @@ def linkcode_resolve(domain, info): nbsphinx_thumbnails = { "examples/tutorial_parcels_structure": "_images/parcels_user_diagram.png", "examples/tutorial_timestamps": "_static/calendar-icon.jpg", - "examples/tutorial_jit_vs_scipy": "_static/clock-icon.png", "examples/documentation_homepage_animation": "_images/homepage.gif", "examples/tutorial_interaction": "_static/pulled_particles_twoatractors_line.gif", "examples/documentation_LargeRunsOutput": "_static/harddrive.png", diff --git a/docs/documentation/index.rst b/docs/documentation/index.rst index 9c849b84fc..c44b532aa3 100644 --- a/docs/documentation/index.rst +++ b/docs/documentation/index.rst @@ -32,7 +32,6 @@ Parcels has several documentation and tutorial Jupyter notebooks and scripts whi :caption: Creating ParticleSets :name: tutorial-particlesets - ../examples/tutorial_jit_vs_scipy.ipynb ../examples/tutorial_delaystart.ipynb diff --git a/docs/installation.rst b/docs/installation.rst index f5775847a6..1f89f951c9 100644 --- a/docs/installation.rst +++ b/docs/installation.rst @@ -36,9 +36,6 @@ The steps below are the installation instructions for Linux, macOS and Windows. python example_peninsula.py --fieldset 100 100 -.. note:: - If you are on macOS and get a compilation error, you may need to accept the Apple xcode license ``xcode-select --install``. If this does not solve the compilation error, you may want to try running ``export CC=gcc``. If the compilation error remains, you may want to check `this solution `_. - *Optionally:* if you want to run all the examples and tutorials, start Jupyter and open the tutorial notebooks: .. code-block:: bash diff --git a/docs/reference.rst b/docs/reference.rst index 67062e6cf9..caecbd39a7 100644 --- a/docs/reference.rst +++ b/docs/reference.rst @@ -11,5 +11,4 @@ Parcels API User-defined Kernels Particle-particle interaction Gridsets and grids - C Code Generation Miscellaneous diff --git a/docs/reference/code_generation.rst b/docs/reference/code_generation.rst deleted file mode 100644 index 94337f149c..0000000000 --- a/docs/reference/code_generation.rst +++ /dev/null @@ -1,16 +0,0 @@ -C Code Generation -================= - -parcels.compilation.codegenerator module ----------------------------------------- - -.. automodule:: parcels.compilation.codegenerator - :members: - :show-inheritance: yes - -parcels.compilation.codecompiler module ---------------------------------------- - -.. automodule:: parcels.compilation.codecompiler - :members: - :show-inheritance: yes diff --git a/parcels/compilation/__init__.py b/parcels/compilation/__init__.py deleted file mode 100644 index e69de29bb2..0000000000 diff --git a/parcels/compilation/codecompiler.py b/parcels/compilation/codecompiler.py deleted file mode 100644 index 4237c87529..0000000000 --- a/parcels/compilation/codecompiler.py +++ /dev/null @@ -1,314 +0,0 @@ -import os -import subprocess -from struct import calcsize - -from parcels._compat import MPI - -_tmp_dir = os.getcwd() - - -class Compiler_parameters: - def __init__(self): - self._compiler = "" - self._cppargs = [] - self._ldargs = [] - self._incdirs = [] - self._libdirs = [] - self._libs = [] - self._dynlib_ext = "" - self._stclib_ext = "" - self._obj_ext = "" - self._exe_ext = "" - - @property - def compiler(self): - return self._compiler - - @property - def cppargs(self): - return self._cppargs - - @property - def ldargs(self): - return self._ldargs - - @property - def incdirs(self): - return self._incdirs - - @property - def libdirs(self): - return self._libdirs - - @property - def libs(self): - return self._libs - - @property - def dynlib_ext(self): - return self._dynlib_ext - - @property - def stclib_ext(self): - return self._stclib_ext - - @property - def obj_ext(self): - return self._obj_ext - - @property - def exe_ext(self): - return self._exe_ext - - -class GNU_parameters(Compiler_parameters): - def __init__(self, cppargs=None, ldargs=None, incdirs=None, libdirs=None, libs=None): - super().__init__() - if cppargs is None: - cppargs = [] - if ldargs is None: - ldargs = [] - if incdirs is None: - incdirs = [] - if libdirs is None: - libdirs = [] - if libs is None: - libs = [] - libs.append("m") - - Iflags = [] - if isinstance(incdirs, list): - for dir in incdirs: - Iflags.append("-I" + dir) - Lflags = [] - if isinstance(libdirs, list): - for dir in libdirs: - Lflags.append("-L" + dir) - lflags = [] - if isinstance(libs, list): - for lib in libs: - lflags.append("-l" + lib) - - cc_env = os.getenv("CC") - mpicc = None - if MPI: - mpicc_env = os.getenv("MPICC") - mpicc = mpicc_env - mpicc = "mpicc" if mpicc is None and os._exists("mpicc") else None - mpicc = "mpiCC" if mpicc is None and os._exists("mpiCC") else None - self._compiler = mpicc if MPI and mpicc is not None else cc_env if cc_env is not None else "gcc" - opt_flags = ["-g", "-O3"] - arch_flag = ["-m64" if calcsize("P") == 8 else "-m32"] - self._cppargs = ["-Wall", "-fPIC", "-std=gnu11"] - self._cppargs += Iflags - self._cppargs += opt_flags + cppargs + arch_flag - self._ldargs = ["-shared"] - self._ldargs += Lflags - self._ldargs += lflags - self._ldargs += ldargs - if len(Lflags) > 0: - self._ldargs += [f"-Wl, -rpath={':'.join(libdirs)}"] - self._ldargs += arch_flag - self._incdirs = incdirs - self._libdirs = libdirs - self._libs = libs - self._dynlib_ext = "so" - self._stclib_ext = "a" - self._obj_ext = "o" - self._exe_ext = "" - - -class Clang_parameters(Compiler_parameters): - def __init__(self, cppargs=None, ldargs=None, incdirs=None, libdirs=None, libs=None): - super().__init__() - if cppargs is None: - cppargs = [] - if ldargs is None: - ldargs = [] - if incdirs is None: - incdirs = [] - if libdirs is None: - libdirs = [] - if libs is None: - libs = [] - self._compiler = "cc" - self._cppargs = cppargs - self._ldargs = ldargs - self._incdirs = incdirs - self._libdirs = libdirs - self._libs = libs - self._dynlib_ext = "dynlib" - self._stclib_ext = "a" - self._obj_ext = "o" - self._exe_ext = "exe" - - -class MinGW_parameters(Compiler_parameters): - def __init__(self, cppargs=None, ldargs=None, incdirs=None, libdirs=None, libs=None): - super().__init__() - if cppargs is None: - cppargs = [] - if ldargs is None: - ldargs = [] - if incdirs is None: - incdirs = [] - if libdirs is None: - libdirs = [] - if libs is None: - libs = [] - self._compiler = "gcc" - self._cppargs = cppargs - self._ldargs = ldargs - self._incdirs = incdirs - self._libdirs = libdirs - self._libs = libs - self._dynlib_ext = "so" - self._stclib_ext = "a" - self._obj_ext = "o" - self._exe_ext = "exe" - - -class VS_parameters(Compiler_parameters): - def __init__(self, cppargs=None, ldargs=None, incdirs=None, libdirs=None, libs=None): - super().__init__() - if cppargs is None: - cppargs = [] - if ldargs is None: - ldargs = [] - if incdirs is None: - incdirs = [] - if libdirs is None: - libdirs = [] - if libs is None: - libs = [] - self._compiler = "cl" - self._cppargs = cppargs - self._ldargs = ldargs - self._incdirs = incdirs - self._libdirs = libdirs - self._libs = libs - self._dynlib_ext = "dll" - self._stclib_ext = "lib" - self._obj_ext = "obj" - self._exe_ext = "exe" - - -class CCompiler: - """A compiler object for creating and loading shared libraries. - - Parameters - ---------- - cc : - C compiler executable (uses environment variable ``CC`` if not provided). - cppargs : - A list of arguments to the C compiler (optional). - ldargs : - A list of arguments to the linker (optional). - """ - - def __init__(self, cc=None, cppargs=None, ldargs=None, incdirs=None, libdirs=None, libs=None, tmp_dir=None): - if tmp_dir is None: - tmp_dir = _tmp_dir - if cppargs is None: - cppargs = [] - if ldargs is None: - ldargs = [] - - self._cc = os.getenv("CC") if cc is None else cc - self._cppargs = cppargs - self._ldargs = ldargs - self._dynlib_ext = "" - self._stclib_ext = "" - self._obj_ext = "" - self._exe_ext = "" - self._tmp_dir = tmp_dir - self._incdirs = incdirs - self._libdirs = libdirs # only possible for already-compiled, external libraries - self._libs = libs # only possible for already-compiled, external libraries - - def compile(self, src, obj, log): - pass - - def _create_compile_process_(self, cmd, src, log): - with open(log, "w") as logfile: - try: - subprocess.check_call(cmd, stdout=logfile, stderr=logfile) - except OSError: - raise RuntimeError(f"OSError during compilation. Please check if compiler exists: {self._cc}") - except subprocess.CalledProcessError: - with open(log) as logfile2: - raise RuntimeError( - f"Error during compilation:\n" - f"Compilation command: {cmd}\n" - f"Source/Destination file: {src}\n" - f"Log file: {logfile.name}\n" - f"Log output: {logfile2.read()}\n" - f"\n" - f"If you are on macOS, it might help to type 'export CC=gcc'" - ) - return True - - -class CCompiler_SS(CCompiler): - """Single-stage C-compiler; used for a SINGLE source file.""" - - def __init__(self, cc=None, cppargs=None, ldargs=None, incdirs=None, libdirs=None, libs=None, tmp_dir=None): - super().__init__( - cc=cc, cppargs=cppargs, ldargs=ldargs, incdirs=incdirs, libdirs=libdirs, libs=libs, tmp_dir=tmp_dir - ) - - def __str__(self): - output = "[CCompiler_SS]: " - output += f"('cc': {self._cc}), " - output += f"('cppargs': {self._cppargs}), " - output += f"('ldargs': {self._ldargs}), " - output += f"('incdirs': {self._incdirs}), " - output += f"('libdirs': {self._libdirs}), " - output += f"('libs': {self._libs}), " - output += f"('tmp_dir': {self._tmp_dir}), " - return output - - def compile(self, src, obj, log): - cc = [self._cc] + self._cppargs + ["-o", obj, src] + self._ldargs - with open(log, "w") as logfile: - logfile.write(f"Compiling: {cc}\n") - self._create_compile_process_(cc, src, log) - - -class GNUCompiler_SS(CCompiler_SS): - """A compiler object for the GNU Linux toolchain. - - Parameters - ---------- - cppargs : - A list of arguments to pass to the C compiler - (optional). - ldargs : - A list of arguments to pass to the linker (optional). - """ - - def __init__(self, cppargs=None, ldargs=None, incdirs=None, libdirs=None, libs=None, tmp_dir=None): - c_params = GNU_parameters(cppargs, ldargs, incdirs, libdirs, libs) - super().__init__( - c_params.compiler, - cppargs=c_params.cppargs, - ldargs=c_params.ldargs, - incdirs=c_params.incdirs, - libdirs=c_params.libdirs, - libs=c_params.libs, - tmp_dir=tmp_dir, - ) - self._dynlib_ext = c_params.dynlib_ext - self._stclib_ext = c_params.stclib_ext - self._obj_ext = c_params.obj_ext - self._exe_ext = c_params.exe_ext - - def compile(self, src, obj, log): - lib_pathfile = os.path.basename(obj) - lib_pathdir = os.path.dirname(obj) - obj = os.path.join(lib_pathdir, lib_pathfile) - - super().compile(src, obj, log) - - -GNUCompiler = GNUCompiler_SS diff --git a/parcels/include/index_search.h b/parcels/include/index_search.h deleted file mode 100644 index 388ba4fa50..0000000000 --- a/parcels/include/index_search.h +++ /dev/null @@ -1,523 +0,0 @@ -#ifndef _INDEX_SEARCH_H -#define _INDEX_SEARCH_H -#ifdef __cplusplus -extern "C" { -#endif - -#include -#include -#include - -#define CHECKSTATUS(res) do {if (res != SUCCESS) return res;} while (0) -#define CHECKSTATUS_KERNELLOOP(res) {if (res == REPEAT) return res;} -#define rtol 1.e-5 -#define atol 1.e-8 - -#ifdef DOUBLE_COORD_VARIABLES -typedef double type_coord; -#else -typedef float type_coord; -#endif - -typedef enum - { - LINEAR=0, NEAREST=1, CGRID_VELOCITY=2, CGRID_TRACER=3, BGRID_VELOCITY=4, BGRID_W_VELOCITY=5, BGRID_TRACER=6, LINEAR_INVDIST_LAND_TRACER=7, PARTIALSLIP=8, FREESLIP=9 - } InterpCode; - -typedef enum - { - NEMO = 0, MITGCM = 1, MOM5 = 2, POP = 3, CROCO = 4 - } GridIndexingType; - -typedef struct -{ - int gtype; - void *grid; -} CGrid; - -typedef struct -{ - int xdim, ydim, zdim, tdim, z4d; - int sphere_mesh, zonal_periodic; - int *chunk_info; - int *load_chunk; - double tfull_min, tfull_max; - int* periods; - float *lonlat_minmax; - float *lon, *lat, *depth; - double *time; -} CStructuredGrid; - - -typedef enum - { - SUCCESS=0, EVALUATE=10, REPEAT=20, DELETE=30, STOPEXECUTION=40, STOPALLEXECUTION=41, ERROR=50, ERRORINTERPOLATION=51, ERROROUTOFBOUNDS=60, ERRORTHROUGHSURFACE=61, ERRORTIMEEXTRAPOLATION=70 - } StatusCode; - -typedef enum - { - RECTILINEAR_Z_GRID=0, RECTILINEAR_S_GRID=1, CURVILINEAR_Z_GRID=2, CURVILINEAR_S_GRID=3 - } GridType; - -// equal/closeness comparison that is equal to numpy (double) -static inline bool is_close_dbl(double a, double b) { - return (fabs(a-b) <= (atol + rtol * fabs(b))); -} - -// customisable equal/closeness comparison (double) -static inline bool is_close_dbl_tol(double a, double b, double tolerance) { - return (fabs(a-b) <= (tolerance + fabs(b))); -} - -// numerically accurate equal/closeness comparison (double) -static inline bool is_equal_dbl(double a, double b) { - return (fabs(a-b) <= (DBL_EPSILON * fabs(b))); -} - -// customisable equal/closeness comparison (float) -static inline bool is_close_flt_tol(float a, float b, float tolerance) { - return (fabs(a-b) <= (tolerance + fabs(b))); -} - -// equal/closeness comparison that is equal to numpy (float) -static inline bool is_close_flt(float a, float b) { - return (fabs(a-b) <= ((float)(atol) + (float)(rtol) * fabs(b))); -} - -// numerically accurate equal/closeness comparison (float) -static inline bool is_equal_flt(float a, float b) { - return (fabs(a-b) <= (FLT_EPSILON * fabs(b))); -} - -static inline bool is_zero_dbl(double a) { - return (fabs(a) <= DBL_EPSILON * fabs(a)); -} - -static inline bool is_zero_flt(float a) { - return (fabs(a) <= FLT_EPSILON * fabs(a)); -} - -static inline StatusCode search_indices_vertical_z(type_coord z, int zdim, float *zvals, int *zi, double *zeta, int gridindexingtype) -{ - if (zvals[zdim-1] > zvals[0]){ - if ((z < zvals[0]) && (gridindexingtype == MOM5) && (z > 2 * zvals[0] - zvals[1])){ - *zi = -1; - *zeta = z / zvals[0]; - return SUCCESS; - } - if ((z > zvals[zdim-1]) && (z < 0) && (gridindexingtype == CROCO)){ - *zi = zdim-2; - *zeta = 1; - return SUCCESS; - } - if (z < zvals[0]) {return ERRORTHROUGHSURFACE;} - if (z > zvals[zdim-1]) {return ERROROUTOFBOUNDS;} - while (*zi < zdim-1 && z > zvals[*zi+1]) ++(*zi); - while (*zi > 0 && z < zvals[*zi]) --(*zi); - } - else{ - if (z > zvals[0]) {return ERRORTHROUGHSURFACE;} - if (z < zvals[zdim-1]) {return ERROROUTOFBOUNDS;} - while (*zi < zdim-1 && z < zvals[*zi+1]) ++(*zi); - while (*zi > 0 && z > zvals[*zi]) --(*zi); - } - if (*zi == zdim-1) {--*zi;} - - *zeta = (z - zvals[*zi]) / (zvals[*zi+1] - zvals[*zi]); - return SUCCESS; -} - -static inline StatusCode search_indices_vertical_s(double time, type_coord z, - int tdim, int zdim, int ydim, int xdim, float *zvals, - int ti, int *zi, int yi, int xi, double *zeta, double eta, double xsi, - int z4d, double t0, double t1, int interp_method) -{ - if (interp_method == BGRID_VELOCITY || interp_method == BGRID_W_VELOCITY || interp_method == BGRID_TRACER){ - xsi = 1; - eta = 1; - } - float zcol[zdim]; - int zii; - if (z4d == 1){ - float (*zvalstab)[zdim][ydim][xdim] = (float (*)[zdim][ydim][xdim]) zvals; - int ti1 = ti; - if (ti < tdim-1) - ti1= ti+1; - double zt0, zt1; - for (zii=0; zii < zdim; zii++){ - zt0 = (1-xsi)*(1-eta) * zvalstab[ti ][zii][yi ][xi ] - + ( xsi)*(1-eta) * zvalstab[ti ][zii][yi ][xi+1] - + ( xsi)*( eta) * zvalstab[ti ][zii][yi+1][xi+1] - + (1-xsi)*( eta) * zvalstab[ti ][zii][yi+1][xi ]; - zt1 = (1-xsi)*(1-eta) * zvalstab[ti1][zii][yi ][xi ] - + ( xsi)*(1-eta) * zvalstab[ti1][zii][yi ][xi+1] - + ( xsi)*( eta) * zvalstab[ti1][zii][yi+1][xi+1] - + (1-xsi)*( eta) * zvalstab[ti1][zii][yi+1][xi ]; - zcol[zii] = zt0 + (zt1 - zt0) * (float)((time - t0) / (t1 - t0)); - } - - } - else{ - float (*zvalstab)[ydim][xdim] = (float (*)[ydim][xdim]) zvals; - for (zii=0; zii < zdim; zii++){ - zcol[zii] = (1-xsi)*(1-eta) * zvalstab[zii][yi ][xi ] - + ( xsi)*(1-eta) * zvalstab[zii][yi ][xi+1] - + ( xsi)*( eta) * zvalstab[zii][yi+1][xi+1] - + (1-xsi)*( eta) * zvalstab[zii][yi+1][xi ]; - } - } - - if (zcol[zdim-1] > zcol[0]){ - if (z < zcol[0]) {return ERRORTHROUGHSURFACE;} - if (z > zcol[zdim-1]) {return ERROROUTOFBOUNDS;} - while (*zi < zdim-1 && z > zcol[*zi+1]) ++(*zi); - while (*zi > 0 && z < zcol[*zi]) --(*zi); - } - else{ - if (z > zcol[0]) {return ERRORTHROUGHSURFACE;} - if (z < zcol[zdim-1]) {return ERROROUTOFBOUNDS;} - while (*zi < zdim-1 && z < zcol[*zi+1]) ++(*zi); - while (*zi > 0 && z > zcol[*zi]) --(*zi); - } - if (*zi == zdim-1) {--*zi;} - - *zeta = (z - zcol[*zi]) / (zcol[*zi+1] - zcol[*zi]); - return SUCCESS; -} - -static inline void reconnect_bnd_indices(int *yi, int *xi, int ydim, int xdim, int onlyX, int sphere_mesh) -{ - if (*xi < 0){ - if (sphere_mesh) - (*xi) = xdim-2; - else - (*xi) = 0; - } - if (*xi > xdim-2){ - if (sphere_mesh) - (*xi) = 0; - else - (*xi) = xdim-2; - } - if (onlyX == 0){ - if (*yi < 0){ - (*yi) = 0; - } - if (*yi > ydim-2){ - (*yi) = ydim-2; - if (sphere_mesh) - (*xi) = xdim - (*xi); - } - } -} - - -static inline StatusCode search_indices_rectilinear(double time, type_coord z, type_coord y, type_coord x, - CStructuredGrid *grid, GridType gtype, - int ti, int *zi, int *yi, int *xi, double *zeta, double *eta, double *xsi, - double t0, double t1, int interp_method, - int gridindexingtype) -{ - int xdim = grid->xdim; - int ydim = grid->ydim; - int zdim = grid->zdim; - int tdim = grid->tdim; - float *xvals = grid->lon; - float *yvals = grid->lat; - float *zvals = grid->depth; - float *xy_minmax = grid->lonlat_minmax; - int sphere_mesh = grid->sphere_mesh; - int zonal_periodic = grid->zonal_periodic; - int z4d = grid->z4d; - - if (zonal_periodic == 0){ - if ((xdim > 1) && ((x < xy_minmax[0]) || (x > xy_minmax[1]))) - return ERROROUTOFBOUNDS; - } - if ((ydim > 1) && ((y < xy_minmax[2]) || (y > xy_minmax[3]))) - return ERROROUTOFBOUNDS; - - if (xdim == 1){ - *xi = 0; - *xsi = 0; - } - else if (sphere_mesh == 0){ - while (*xi < xdim-1 && x > xvals[*xi+1]) ++(*xi); - while (*xi > 0 && x < xvals[*xi]) --(*xi); - *xsi = (x - xvals[*xi]) / (xvals[*xi+1] - xvals[*xi]); - } - else{ - - float xvalsi = xvals[*xi]; - // TODO: this will fail if longitude is e.g. only [-180, 180] (so length 2) - if (xvalsi < x - 225) xvalsi += 360; - if (xvalsi > x + 225) xvalsi -= 360; - float xvalsi1 = xvals[*xi+1]; - if (xvalsi1 < xvalsi - 180) xvalsi1 += 360; - if (xvalsi1 > xvalsi + 180) xvalsi1 -= 360; - - int itMax = 10000; - int it = 0; - while ( (xvalsi > x) || (xvalsi1 < x) ){ - if (xvalsi1 < x) - ++(*xi); - else if (xvalsi > x) - --(*xi); - reconnect_bnd_indices(yi, xi, ydim, xdim, 1, 1); - xvalsi = xvals[*xi]; - if (xvalsi < x - 225) xvalsi += 360; - if (xvalsi > x + 225) xvalsi -= 360; - xvalsi1 = xvals[*xi+1]; - if (xvalsi1 < xvalsi - 180) xvalsi1 += 360; - if (xvalsi1 > xvalsi + 180) xvalsi1 -= 360; - it++; - if (it > itMax){ - return ERROROUTOFBOUNDS; - } - } - - *xsi = (x - xvalsi) / (xvalsi1 - xvalsi); - } - - if (ydim == 1){ - *yi = 0; - *eta = 0; - } - else { - while (*yi < ydim-1 && y > yvals[*yi+1]) ++(*yi); - while (*yi > 0 && y < yvals[*yi]) --(*yi); - *eta = (y - yvals[*yi]) / (yvals[*yi+1] - yvals[*yi]); - } - - StatusCode status; - if (zdim > 1){ - switch(gtype){ - case RECTILINEAR_Z_GRID: - status = search_indices_vertical_z(z, zdim, zvals, zi, zeta, gridindexingtype); - break; - case RECTILINEAR_S_GRID: - status = search_indices_vertical_s(time, z, tdim, zdim, ydim, xdim, zvals, - ti, zi, *yi, *xi, zeta, *eta, *xsi, - z4d, t0, t1, interp_method); - break; - default: - status = ERRORINTERPOLATION; - } - CHECKSTATUS(status); - } - else - *zeta = 0; - - if ( (*xsi < 0) && (is_zero_dbl(*xsi)) ) {*xsi = 0.;} - if ( (*xsi > 1) && (is_close_dbl(*xsi, 1.)) ) {*xsi = 1.;} - if ( (*eta < 0) && (is_zero_dbl(*eta)) ) {*eta = 0.;} - if ( (*eta > 1) && (is_close_dbl(*eta, 1.)) ) {*eta = 1.;} - if ( (*zeta < 0) && (is_zero_dbl(*zeta)) ) {*zeta = 0.;} - if ( (*zeta > 1) && (is_close_dbl(*zeta, 1.)) ) {*zeta = 1.;} - - if ( (*xsi < 0) || (*xsi > 1) ) return ERRORINTERPOLATION; - if ( (*eta < 0) || (*eta > 1) ) return ERRORINTERPOLATION; - if ( (*zeta < 0) || (*zeta > 1) ) return ERRORINTERPOLATION; - - return SUCCESS; -} - - -static inline StatusCode search_indices_curvilinear(double time, type_coord z, type_coord y, type_coord x, - CStructuredGrid *grid, GridType gtype, - int ti, int *zi, int *yi, int *xi, - double *zeta, double *eta, double *xsi, - double t0, double t1, int interp_method, - int gridindexingtype) -{ - int xi_old = *xi; - int yi_old = *yi; - int xdim = grid->xdim; - int ydim = grid->ydim; - int zdim = grid->zdim; - int tdim = grid->tdim; - float *xvals = grid->lon; - float *yvals = grid->lat; - float *zvals = grid->depth; - float *xy_minmax = grid->lonlat_minmax; - int sphere_mesh = grid->sphere_mesh; - int zonal_periodic = grid->zonal_periodic; - int z4d = grid->z4d; - - // NEMO convention - float (* xgrid)[xdim] = (float (*)[xdim]) xvals; - float (* ygrid)[xdim] = (float (*)[xdim]) yvals; - - if (zonal_periodic == 0){ - if ((x < xy_minmax[0]) || (x > xy_minmax[1])){ - if (xgrid[0][0] < xgrid[0][xdim-1]) {return ERROROUTOFBOUNDS;} - else if (x < xgrid[0][0] && x > xgrid[0][xdim-1]) {return ERROROUTOFBOUNDS;} - } - } - if ((y < xy_minmax[2]) || (y > xy_minmax[3])) - return ERROROUTOFBOUNDS; - - double a[4], b[4]; - - *xsi = *eta = -1; - int maxIterSearch = 1e6, it = 0; - double tol = 1e-10; - while ( (*xsi < -tol) || (*xsi > 1+tol) || (*eta < -tol) || (*eta > 1+tol) ){ - double xgrid_loc[4] = {xgrid[*yi][*xi], xgrid[*yi][*xi+1], xgrid[*yi+1][*xi+1], xgrid[*yi+1][*xi]}; - if (sphere_mesh){ //we are on the sphere - int i4; - if (xgrid_loc[0] < x - 225) xgrid_loc[0] += 360; - if (xgrid_loc[0] > x + 225) xgrid_loc[0] -= 360; - for (i4 = 1; i4 < 4; ++i4){ - if (xgrid_loc[i4] < xgrid_loc[0] - 180) xgrid_loc[i4] += 360; - if (xgrid_loc[i4] > xgrid_loc[0] + 180) xgrid_loc[i4] -= 360; - } - } - double ygrid_loc[4] = {ygrid[*yi][*xi], ygrid[*yi][*xi+1], ygrid[*yi+1][*xi+1], ygrid[*yi+1][*xi]}; - - a[0] = xgrid_loc[0]; - a[1] = -xgrid_loc[0] + xgrid_loc[1]; - a[2] = -xgrid_loc[0] + xgrid_loc[3]; - a[3] = xgrid_loc[0] - xgrid_loc[1] + xgrid_loc[2] - xgrid_loc[3]; - b[0] = ygrid_loc[0]; - b[1] = -ygrid_loc[0] + ygrid_loc[1]; - b[2] = -ygrid_loc[0] + ygrid_loc[3]; - b[3] = ygrid_loc[0] - ygrid_loc[1] + ygrid_loc[2] - ygrid_loc[3]; - - double aa = a[3]*b[2] - a[2]*b[3]; - double bb = a[3]*b[0] - a[0]*b[3] + a[1]*b[2] - a[2]*b[1] + x*b[3] - y*a[3]; - double cc = a[1]*b[0] - a[0]*b[1] + x*b[1] - y*a[1]; - if (fabs(aa) < 1e-12) // Rectilinear cell, or quasi - *eta = -cc / bb; - else{ - double det = sqrt(bb*bb-4*aa*cc); - if (det == det) // so, if det is nan we keep the xsi, eta from previous iter - *eta = (-bb+det)/(2*aa); - } - if ( fabs(a[1]+a[3]*(*eta)) < 1e-12 ) // this happens when recti cell rotated of 90deg - *xsi = ( (y-ygrid_loc[0]) / (ygrid_loc[1]-ygrid_loc[0]) + - (y-ygrid_loc[3]) / (ygrid_loc[2]-ygrid_loc[3]) ) * .5; - else - *xsi = (x-a[0]-a[2]* (*eta)) / (a[1]+a[3]* (*eta)); - if ( (*xsi < 0) && (*eta < 0) && (*xi == 0) && (*yi == 0) ) - return ERROROUTOFBOUNDS; - if ( (*xsi > 1) && (*eta > 1) && (*xi == xdim-1) && (*yi == ydim-1) ) - return ERROROUTOFBOUNDS; - if (*xsi < -tol) - (*xi)--; - if (*xsi > 1+tol) - (*xi)++; - if (*eta < -tol) - (*yi)--; - if (*eta > 1+tol) - (*yi)++; - reconnect_bnd_indices(yi, xi, ydim, xdim, 0, sphere_mesh); - it++; - if ( it > maxIterSearch){ - printf("Correct cell not found for (lat, lon) = (%f, %f) after %d iterations\n", y, x, maxIterSearch); - printf("Debug info: old particle indices: (yi, xi) %d %d\n", yi_old, xi_old); - printf(" new particle indices: (yi, xi) %d %d\n", *yi, *xi); - printf(" Mesh 2d shape: %d %d\n", ydim, xdim); - printf(" Relative particle position: (eta, xsi) %1.16e %1.16e\n", *eta, *xsi); - return ERROROUTOFBOUNDS; - } - } - if ( (*xsi != *xsi) || (*eta != *eta) ){ // check if nan - printf("Correct cell not found for (lat, lon) = (%f, %f))\n", y, x); - printf("Debug info: old particle indices: (yi, xi) %d %d\n", yi_old, xi_old); - printf(" new particle indices: (yi, xi) %d %d\n", *yi, *xi); - printf(" Mesh 2d shape: %d %d\n", ydim, xdim); - printf(" Relative particle position: (eta, xsi) %1.16e %1.16e\n", *eta, *xsi); - return ERROROUTOFBOUNDS; - } - if (*xsi < 0) *xsi = 0; - if (*xsi > 1) *xsi = 1; - if (*eta < 0) *eta = 0; - if (*eta > 1) *eta = 1; - - StatusCode status; - if (zdim > 1){ - switch(gtype){ - case CURVILINEAR_Z_GRID: - status = search_indices_vertical_z(z, zdim, zvals, zi, zeta, gridindexingtype); - break; - case CURVILINEAR_S_GRID: - status = search_indices_vertical_s(time, z, tdim, ydim, xdim, zdim, zvals, - ti, zi, *yi, *xi, zeta, *eta, *xsi, - z4d, t0, t1, interp_method); - break; - default: - status = ERRORINTERPOLATION; - } - CHECKSTATUS(status); - } - else - *zeta = 0; - - if ( (*xsi < 0) || (*xsi > 1) ) return ERRORINTERPOLATION; - if ( (*eta < 0) || (*eta > 1) ) return ERRORINTERPOLATION; - if ( (*zeta < 0) || (*zeta > 1) ) return ERRORINTERPOLATION; - - return SUCCESS; -} - -/* Local linear search to update grid index - * params ti, sizeT, time. t0, t1 are only used for 4D S grids - * */ -static inline StatusCode search_indices(double time, type_coord z, type_coord y, type_coord x, - CStructuredGrid *grid, - int ti, int *zi, int *yi, int *xi, - double *zeta, double *eta, double *xsi, - GridType gtype, double t0, double t1, int interp_method, - int gridindexingtype) -{ - switch(gtype){ - case RECTILINEAR_Z_GRID: - case RECTILINEAR_S_GRID: - return search_indices_rectilinear(time, z, y, x, grid, gtype, ti, zi, yi, xi, zeta, eta, xsi, - t0, t1, interp_method, gridindexingtype); - break; - case CURVILINEAR_Z_GRID: - case CURVILINEAR_S_GRID: - return search_indices_curvilinear(time, z, y, x, grid, gtype, ti, zi, yi, xi, zeta, eta, xsi, - t0, t1, interp_method, gridindexingtype); - break; - default: - printf("Only RECTILINEAR_Z_GRID, RECTILINEAR_S_GRID, CURVILINEAR_Z_GRID and CURVILINEAR_S_GRID grids are currently implemented\n"); - return ERROR; - } -} - -/* Local linear search to update time index */ -static inline StatusCode search_time_index(double *t, int size, double *tvals, int *ti, int time_periodic, double tfull_min, double tfull_max, int *periods) -{ - if (*ti < 0) - *ti = 0; - if (time_periodic == 1){ - if (*t < tvals[0]){ - *ti = size-1; - *periods = (int) floor( (*t-tfull_min)/(tfull_max-tfull_min)); - *t -= *periods * (tfull_max-tfull_min); - if (*t < tvals[0]){ // e.g. t=5, tfull_min=0, t_full_max=5 -> periods=1 but we want periods = 0 - *periods -= 1; - *t -= *periods * (tfull_max-tfull_min); - } - search_time_index(t, size, tvals, ti, time_periodic, tfull_min, tfull_max, periods); - } - else if (*t > tvals[size-1]){ - *ti = 0; - *periods = (int) floor( (*t-tfull_min)/(tfull_max-tfull_min)); - *t -= *periods * (tfull_max-tfull_min); - search_time_index(t, size, tvals, ti, time_periodic, tfull_min, tfull_max, periods); - } - } - while (*ti < size-1 && *t > tvals[*ti+1]) ++(*ti); - while (*ti > 0 && *t < tvals[*ti]) --(*ti); - return SUCCESS; -} - - -#ifdef __cplusplus -} -#endif -#endif diff --git a/parcels/include/interpolation_utils.h b/parcels/include/interpolation_utils.h deleted file mode 100644 index b92fd33273..0000000000 --- a/parcels/include/interpolation_utils.h +++ /dev/null @@ -1,139 +0,0 @@ -#include -#include -#include - -typedef enum - { - ZONAL=0, MERIDIONAL=1, VERTICAL=2, - } Orientation; - - -static inline void phi2D_lin(double eta, double xsi, double *phi) -{ - phi[0] = (1-xsi) * (1-eta); - phi[1] = xsi * (1-eta); - phi[2] = xsi * eta ; - phi[3] = (1-xsi) * eta ; -} - - -static inline void phi1D_quad(double xsi, double *phi) -{ - phi[0] = 2*xsi*xsi-3*xsi+1; - phi[1] = -4*xsi*xsi+4*xsi; - phi[2] = 2*xsi*xsi-xsi; -} - - -static inline void dphidxsi3D_lin(double zeta, double eta, double xsi, double *dphidzeta, double *dphideta, double *dphidxsi) -{ - dphidxsi[0] = - (1-eta) * (1-zeta); - dphidxsi[1] = (1-eta) * (1-zeta); - dphidxsi[2] = ( eta) * (1-zeta); - dphidxsi[3] = - ( eta) * (1-zeta); - dphidxsi[4] = - (1-eta) * ( zeta); - dphidxsi[5] = (1-eta) * ( zeta); - dphidxsi[6] = ( eta) * ( zeta); - dphidxsi[7] = - ( eta) * ( zeta); - - dphideta[0] = - (1-xsi) * (1-zeta); - dphideta[1] = - ( xsi) * (1-zeta); - dphideta[2] = ( xsi) * (1-zeta); - dphideta[3] = (1-xsi) * (1-zeta); - dphideta[4] = - (1-xsi) * ( zeta); - dphideta[5] = - ( xsi) * ( zeta); - dphideta[6] = ( xsi) * ( zeta); - dphideta[7] = (1-xsi) * ( zeta); - - dphidzeta[0] = - (1-xsi) * (1-eta); - dphidzeta[1] = - ( xsi) * (1-eta); - dphidzeta[2] = - ( xsi) * ( eta); - dphidzeta[3] = - (1-xsi) * ( eta); - dphidzeta[4] = (1-xsi) * (1-eta); - dphidzeta[5] = ( xsi) * (1-eta); - dphidzeta[6] = ( xsi) * ( eta); - dphidzeta[7] = (1-xsi) * ( eta); -} - -static inline void dxdxsi3D_lin(double *pz, double *py, double *px, double zeta, double eta, double xsi, double *jacM, int sphere_mesh) -{ - double dphidxsi[8], dphideta[8], dphidzeta[8]; - dphidxsi3D_lin(zeta, eta, xsi, dphidzeta, dphideta, dphidxsi); - - int i; - for(i=0; i<9; ++i) - jacM[i] = 0; - - double deg2m = 1852 * 60.; - double rad = M_PI / 180.; - double lat = (1-xsi) * (1-eta) * py[0]+ - xsi * (1-eta) * py[1]+ - xsi * eta * py[2]+ - (1-xsi) * eta * py[3]; - double jac_lon = (sphere_mesh == 1) ? (deg2m * cos(rad * lat) ) : 1; - double jac_lat = (sphere_mesh == 1) ? deg2m : 1; - - for(i=0; i<8; ++i){ - jacM[3*0+0] += px[i] * dphidxsi[i] * jac_lon; // dxdxsi - jacM[3*0+1] += px[i] * dphideta[i] * jac_lon; // dxdeta - jacM[3*0+2] += px[i] * dphidzeta[i] * jac_lon; // dxdzeta - jacM[3*1+0] += py[i] * dphidxsi[i] * jac_lat; // dydxsi - jacM[3*1+1] += py[i] * dphideta[i] * jac_lat; // dydeta - jacM[3*1+2] += py[i] * dphidzeta[i] * jac_lat; // dydzeta - jacM[3*2+0] += pz[i] * dphidxsi[i]; // dzdxsi - jacM[3*2+1] += pz[i] * dphideta[i]; // dzdeta - jacM[3*2+2] += pz[i] * dphidzeta[i]; // dzdzeta - } -} - -static inline double jacobian3D_lin_face(double *pz, double *py, double *px, - double zeta, double eta, double xsi, - Orientation orientation, int sphere_mesh) -{ - double jacM[9]; - dxdxsi3D_lin(pz, py, px, zeta, eta, xsi, jacM, sphere_mesh); - - double j[3]; - - if (orientation == ZONAL){ - j[0] = jacM[3*1+1]*jacM[3*2+2]-jacM[3*1+2]*jacM[3*2+1]; - j[1] =-jacM[3*0+1]*jacM[3*2+2]+jacM[3*0+2]*jacM[3*2+1]; - j[2] = jacM[3*0+1]*jacM[3*1+2]-jacM[3*0+2]*jacM[3*1+1]; - } - else if (orientation == MERIDIONAL){ - j[0] = jacM[3*1+0]*jacM[3*2+2]-jacM[3*1+2]*jacM[3*2+0]; - j[1] =-jacM[3*0+0]*jacM[3*2+2]+jacM[3*0+2]*jacM[3*2+0]; - j[2] = jacM[3*0+0]*jacM[3*1+2]-jacM[3*0+2]*jacM[3*1+0]; - } - else if (orientation == VERTICAL){ - j[0] = jacM[3*1+0]*jacM[3*2+1]-jacM[3*1+1]*jacM[3*2+0]; - j[1] =-jacM[3*0+0]*jacM[3*2+1]+jacM[3*0+1]*jacM[3*2+0]; - j[2] = jacM[3*0+0]*jacM[3*1+1]-jacM[3*0+1]*jacM[3*1+0]; - } - - return sqrt(j[0]*j[0]+j[1]*j[1]+j[2]*j[2]); -} - -static inline double jacobian3D_lin(double *pz, double *py, double *px, - double zeta, double eta, double xsi, - int sphere_mesh) -{ - double jacM[9]; - dxdxsi3D_lin(pz, py, px, zeta, eta, xsi, jacM, sphere_mesh); - - double jac = jacM[3*0+0] * (jacM[3*1+1]*jacM[3*2+2] - jacM[3*2+1]*jacM[3*1+2]) - - jacM[3*0+1] * (jacM[3*1+0]*jacM[3*2+2] - jacM[3*2+0]*jacM[3*1+2]) - + jacM[3*0+2] * (jacM[3*1+0]*jacM[3*2+1] - jacM[3*2+0]*jacM[3*1+1]); - - - return jac; -} - -static inline double dot_prod(double *a, double *b, size_t n) -{ - double val = 0; - int i = 0; - for(i=0; i -#include -#include -#include -#include -#include "random.h" -#include "index_search.h" -#include "interpolation_utils.h" - -#define min(X, Y) (((X) < (Y)) ? (X) : (Y)) -#define max(X, Y) (((X) > (Y)) ? (X) : (Y)) - -typedef struct -{ - int xdim, ydim, zdim, tdim, igrid, allow_time_extrapolation, time_periodic; - float ****data_chunks; - CGrid *grid; -} CField; - -/* Bilinear interpolation routine for 2D grid */ -static inline StatusCode spatial_interpolation_bilinear(double eta, double xsi, - float data[2][2], float *value) -{ - *value = (1-xsi)*(1-eta) * data[0][0] - + xsi *(1-eta) * data[0][1] - + xsi * eta * data[1][1] - + (1-xsi)* eta * data[1][0]; - return SUCCESS; -} - -/* Bilinear interpolation routine for 2D grid for tracers with squared inverse distance weighting near land*/ -static inline StatusCode spatial_interpolation_bilinear_invdist_land(double eta, double xsi, - float data[2][2], float *value) -{ - int i, j, k, l, nb_land = 0, land[2][2] = {{0}}; - float w_sum = 0.; - // count the number of surrounding land points (assume land is where the value is close to zero) - for (i = 0; i < 2; i++) { - for (j = 0; j < 2; j++) { - if (is_zero_flt(data[i][j])) { - land[i][j] = 1; - nb_land++; - } - else { - // record the coordinates of the last non-land point - // (for the case where this is the only location with valid data) - k = i; - l = j; - } - } - } - switch (nb_land) { - case 0: // no land, use usual routine - return spatial_interpolation_bilinear(eta, xsi, data, value); - case 3: // single non-land point - *value = data[k][l]; - return SUCCESS; - case 4: // only land - *value = 0.; - return SUCCESS; - default: - break; - } - // interpolate with 1 or 2 land points - *value = 0.; - for (i = 0; i < 2; i++) { - for (j = 0; j < 2; j++) { - float distance = pow((xsi - j), 2) + pow((eta - i), 2); - if (is_zero_flt(distance)) { - if (land[i][j] == 1) { // index search led us directly onto land - *value = 0.; - return SUCCESS; - } - else { - *value = data[i][j]; - return SUCCESS; - } - } - else if (land[i][j] == 0) { - *value += data[i][j] / distance; - w_sum += 1 / distance; - } - } - } - *value /= w_sum; - return SUCCESS; -} - -/* Trilinear interpolation routine for 3D grid */ -static inline StatusCode spatial_interpolation_trilinear(double zeta, double eta, double xsi, - float data[2][2][2], float *value) -{ - float f0, f1; - f0 = (1-xsi)*(1-eta) * data[0][0][0] - + xsi *(1-eta) * data[0][0][1] - + xsi * eta * data[0][1][1] - + (1-xsi)* eta * data[0][1][0]; - f1 = (1-xsi)*(1-eta) * data[1][0][0] - + xsi *(1-eta) * data[1][0][1] - + xsi * eta * data[1][1][1] - + (1-xsi)* eta * data[1][1][0]; - *value = (1-zeta) * f0 + zeta * f1; - return SUCCESS; -} - -/* Trilinear interpolation routine for MOM surface 3D grid */ -static inline StatusCode spatial_interpolation_trilinear_surface(double zeta, double eta, double xsi, - float data[2][2][2], float *value) -{ - float f1; - f1 = (1-xsi)*(1-eta) * data[0][0][0] - + xsi *(1-eta) * data[0][0][1] - + xsi * eta * data[0][1][1] - + (1-xsi)* eta * data[0][1][0]; - *value = zeta * f1; - return SUCCESS; -} - -static inline StatusCode spatial_interpolation_trilinear_bottom(double zeta, double eta, double xsi, - float data[2][2][2], float *value) -{ - float f1; - f1 = (1-xsi)*(1-eta) * data[1][0][0] - + xsi *(1-eta) * data[1][0][1] - + xsi * eta * data[1][1][1] - + (1-xsi)* eta * data[1][1][0]; - *value = (1 - zeta) * f1; - return SUCCESS; -} - -/* Trilinear interpolation routine for 3D grid for tracers with squared inverse distance weighting near land*/ -static inline StatusCode spatial_interpolation_trilinear_invdist_land(double zeta, double eta, double xsi, - float data[2][2][2], float *value) -{ - int i, j, k, l, m, n, nb_land = 0, land[2][2][2] = {{{0}}}; - float w_sum = 0.; - // count the number of surrounding land points (assume land is where the value is close to zero) - for (i = 0; i < 2; i++) { - for (j = 0; j < 2; j++) { - for (k = 0; k < 2; k++) { - if(is_zero_flt(data[i][j][k])) { - land[i][j][k] = 1; - nb_land++; - } - else { - // record the coordinates of the last non-land point - // (for the case where this is the only location with valid data) - l = i; - m = j; - n = k; - } - } - } - } - switch (nb_land) { - case 0: // no land, use usual routine - return spatial_interpolation_trilinear(zeta, eta, xsi, data, value); - case 7: // single non-land point - *value = data[l][m][n]; - return SUCCESS; - case 8: // only land - *value = 0.; - return SUCCESS; - default: - break; - } - // interpolate with 1 to 6 land points - *value = 0.; - for (i = 0; i < 2; i++) { - for (j = 0; j < 2; j++) { - for (k = 0; k < 2; k++) { - float distance = pow((zeta - i), 2) + pow((eta - j), 2) + pow((xsi - k), 2); - if (is_zero_flt(distance)) { - if (land[i][j][k] == 1) { - // index search led us directly onto land - *value = 0.; - return SUCCESS; - } else { - *value = data[i][j][k]; - return SUCCESS; - } - } - else if (land[i][j][k] == 0) { - *value += data[i][j][k] / distance; - w_sum += 1 / distance; - } - } - } - } - *value /= w_sum; - return SUCCESS; -} - -/* Nearest neighbor interpolation routine for 2D grid */ -static inline StatusCode spatial_interpolation_nearest2D(double eta, double xsi, - float data[2][2], float *value) -{ - /* Cast data array into data[lat][lon] as per NEMO convention */ - int i, j; - if (xsi < .5) {i = 0;} else {i = 1;} - if (eta < .5) {j = 0;} else {j = 1;} - *value = data[j][i]; - return SUCCESS; -} - -/* C grid interpolation routine for tracers on 2D grid */ -static inline StatusCode spatial_interpolation_tracer_bc_grid_2D(double _eta, double _xsi, - float data[2][2], float *value) -{ - *value = data[1][1]; - return SUCCESS; -} - -/* C grid interpolation routine for tracers on 3D grid */ -static inline StatusCode spatial_interpolation_tracer_bc_grid_3D(double _zeta, double _eta, double _xsi, - float data[2][2][2], float *value) -{ - *value = data[0][1][1]; - return SUCCESS; -} - -static inline StatusCode spatial_interpolation_tracer_bc_grid_bottom(double _zeta, double _eta, double _xsi, - float data[2][2][2], float *value) -{ - *value = data[1][1][1]; - return SUCCESS; -} - -/* Nearest neighbor interpolation routine for 3D grid */ -static inline StatusCode spatial_interpolation_nearest3D(double zeta, double eta, double xsi, - float data[2][2][2], float *value) -{ - int i, j, k; - if (xsi < .5) {i = 0;} else {i = 1;} - if (eta < .5) {j = 0;} else {j = 1;} - if (zeta < .5) {k = 0;} else {k = 1;} - *value = data[k][j][i]; - return SUCCESS; -} - -static inline int getBlock2D(int *chunk_info, int yi, int xi, int *block, int *index_local) -{ - int ndim = chunk_info[0]; - if (ndim != 2) - exit(-1); - int i, j; - - int shape[ndim]; - int index[2] = {yi, xi}; - for(i=0; igrid->grid; - int *chunk_info = grid->chunk_info; - int ndim = chunk_info[0]; - int block[ndim]; - int ilocal[ndim]; - - int tii, yii, xii; - - int blockid = getBlock2D(chunk_info, yi, xi, block, ilocal); - if (grid->load_chunk[blockid] < 2){ - grid->load_chunk[blockid] = 1; - return REPEAT; - } - grid->load_chunk[blockid] = 2; - int zdim = 1; - int ydim = chunk_info[1+ndim+block[0]]; - int yshift = chunk_info[1]; - int xdim = chunk_info[1+ndim+yshift+block[1]]; - - if (((ilocal[0] == ydim-1) && (ydim > 1)) || ((ilocal[1] == xdim-1) && (xdim > 1))) - { - // Cell is on multiple chunks - for (tii=0; tii<2; ++tii){ - for (yii=0; yii<2; ++yii){ - for (xii=0; xii<2; ++xii){ - blockid = getBlock2D(chunk_info, yi+yii, xi+xii, block, ilocal); - if (grid->load_chunk[blockid] < 2){ - grid->load_chunk[blockid] = 1; - return REPEAT; - } - grid->load_chunk[blockid] = 2; - zdim = 1; - ydim = chunk_info[1+ndim+block[0]]; - yshift = chunk_info[1]; - xdim = chunk_info[1+ndim+yshift+block[1]]; - float (*data_block)[zdim][ydim][xdim] = (float (*)[zdim][ydim][xdim]) f->data_chunks[blockid]; - float (*data)[xdim] = (float (*)[xdim]) (data_block[ti+tii]); - cell_data[tii][yii][xii] = data[ilocal[0]][ilocal[1]]; - } - } - if (first_tstep_only == 1) - break; - } - } - else - { - float (*data_block)[zdim][ydim][xdim] = (float (*)[zdim][ydim][xdim]) f->data_chunks[blockid]; - for (tii=0; tii<2; ++tii){ - float (*data)[xdim] = (float (*)[xdim]) (data_block[ti+tii]); - int xiid = ((xdim==1) ? 0 : 1); - int yiid = ((ydim==1) ? 0 : 1); - for (yii=0; yii<2; yii++) - for (xii=0; xii<2; xii++) - cell_data[tii][yii][xii] = data[ilocal[0]+(yii*yiid)][ilocal[1]+(xii*xiid)]; - if (first_tstep_only == 1) - break; - } - } - return SUCCESS; -} - -static inline int getBlock3D(int *chunk_info, int zi, int yi, int xi, int *block, int *index_local) -{ - int ndim = chunk_info[0]; - if (ndim != 3) - exit(-1); - int i, j; - - int shape[ndim]; - int index[3] = {zi, yi, xi}; - for(i=0; igrid->grid; - int *chunk_info = grid->chunk_info; - int ndim = chunk_info[0]; - int block[ndim]; - int ilocal[ndim]; - - int tii, zii, yii, xii; - - int blockid = getBlock3D(chunk_info, zi, yi, xi, block, ilocal); - if (grid->load_chunk[blockid] < 2){ - grid->load_chunk[blockid] = 1; - return REPEAT; - } - grid->load_chunk[blockid] = 2; - int zdim = chunk_info[1+ndim+block[0]]; - int zshift = chunk_info[1]; - int ydim = chunk_info[1+ndim+zshift+block[1]]; - int yshift = chunk_info[1+1]; - int xdim = chunk_info[1+ndim+zshift+yshift+block[2]]; - - if (((ilocal[0] == zdim-1) && zdim > 1) || ((ilocal[1] == ydim-1) && ydim > 1) || ((ilocal[2] == xdim-1) && xdim >1)) - { - // Cell is on multiple chunks\n - for (tii=0; tii<2; ++tii){ - for (zii=0; zii<2; ++zii){ - for (yii=0; yii<2; ++yii){ - for (xii=0; xii<2; ++xii){ - blockid = getBlock3D(chunk_info, zi+zii, yi+yii, xi+xii, block, ilocal); - if (grid->load_chunk[blockid] < 2){ - grid->load_chunk[blockid] = 1; - return REPEAT; - } - grid->load_chunk[blockid] = 2; - zdim = chunk_info[1+ndim+block[0]]; - zshift = chunk_info[1]; - ydim = chunk_info[1+ndim+zshift+block[1]]; - yshift = chunk_info[1+1]; - xdim = chunk_info[1+ndim+zshift+yshift+block[2]]; - float (*data_block)[zdim][ydim][xdim] = (float (*)[zdim][ydim][xdim]) f->data_chunks[blockid]; - float (*data)[ydim][xdim] = (float (*)[ydim][xdim]) (data_block[ti+tii]); - cell_data[tii][zii][yii][xii] = data[ilocal[0]][ilocal[1]][ilocal[2]]; - } - } - } - if (first_tstep_only == 1) - break; - } - } - else - { - float (*data_block)[zdim][ydim][xdim] = (float (*)[zdim][ydim][xdim]) f->data_chunks[blockid]; - for (tii=0; tii<2; ++tii){ - float (*data)[ydim][xdim] = (float (*)[ydim][xdim]) (data_block[ti+tii]); - int xiid = ((xdim==1) ? 0 : 1); - int yiid = ((ydim==1) ? 0 : 1); - int ziid = ((zdim==1) ? 0 : 1); - for (zii=0; zii<2; zii++) - for (yii=0; yii<2; yii++) - for (xii=0; xii<2; xii++) - cell_data[tii][zii][yii][xii] = data[ilocal[0]+(zii*ziid)][ilocal[1]+(yii*yiid)][ilocal[2]+(xii*xiid)]; - if (first_tstep_only == 1) - break; - } - } - return SUCCESS; -} - - -/* Linear interpolation along the time axis */ -static inline StatusCode temporal_interpolation_structured_grid(double time, type_coord z, type_coord y, type_coord x, - CField *f, - GridType gtype, int *ti, int *zi, int *yi, int *xi, - float *value, int interp_method, int gridindexingtype) -{ - StatusCode status; - CStructuredGrid *grid = f->grid->grid; - int igrid = f->igrid; - - /* Find time index for temporal interpolation */ - if (f->time_periodic == 0 && f->allow_time_extrapolation == 0 && (time < grid->time[0] || time > grid->time[grid->tdim-1])){ - return ERRORTIMEEXTRAPOLATION; - } - status = search_time_index(&time, grid->tdim, grid->time, &ti[igrid], f->time_periodic, grid->tfull_min, grid->tfull_max, grid->periods); CHECKSTATUS(status); - - double xsi, eta, zeta; - - float data2D[2][2][2]; - float data3D[2][2][2][2]; - - // if we're in between time indices, and not at the end of the timeseries, - // we'll make sure to interpolate data between the two time values - // otherwise, we'll only use the data at the current time index - int tii = (ti[igrid] < grid->tdim-1 && time > grid->time[ti[igrid]]) ? 2 : 1; - - float val[2] = {0.0f, 0.0f}; - double t0 = grid->time[ti[igrid]]; - // we set our second time bound and search time depending on the - // index critereon above - double t1 = (tii == 2) ? grid->time[ti[igrid]+1] : t0+1; - double tsrch = (tii == 2) ? time : t0; - - status = search_indices(tsrch, z, y, x, grid, - ti[igrid], &zi[igrid], &yi[igrid], &xi[igrid], - &zeta, &eta, &xsi, gtype, - t0, t1, interp_method, gridindexingtype); - CHECKSTATUS(status); - - if (grid->zdim == 1) { - // last param is a flag, which denotes that we only want the first timestep - // (rather than both) - status = getCell2D(f, ti[igrid], yi[igrid], xi[igrid], data2D, tii == 1); CHECKSTATUS(status); - } else { - if ((gridindexingtype == MOM5) && (zi[igrid] == -1)) { - status = getCell3D(f, ti[igrid], 0, yi[igrid], xi[igrid], data3D, tii == 1); CHECKSTATUS(status); - } else if ((gridindexingtype == POP) && (zi[igrid] == grid->zdim-2)) { - status = getCell3D(f, ti[igrid], zi[igrid]-1, yi[igrid], xi[igrid], data3D, tii == 1); CHECKSTATUS(status); - } else { - status = getCell3D(f, ti[igrid], zi[igrid], yi[igrid], xi[igrid], data3D, tii == 1); CHECKSTATUS(status); - } - } - - // define a helper macro that will select the appropriate interpolation method - // depending on whether we need 2D or 3D -#define INTERP(fn_2d, fn_3d) \ - do { \ - if (grid->zdim == 1) { \ - for (int i = 0; i < tii; i++) { \ - status = fn_2d(eta, xsi, data2D[i], &val[i]); \ - CHECKSTATUS(status); \ - } \ - } else { \ - for (int i = 0; i < tii; i++) { \ - status = fn_3d(zeta, eta, xsi, data3D[i], &val[i]); \ - CHECKSTATUS(status); \ - } \ - } \ - } while (0) - - if ((interp_method == LINEAR) || (interp_method == CGRID_VELOCITY) || - (interp_method == BGRID_VELOCITY) || (interp_method == BGRID_W_VELOCITY)) { - // adjust the normalised coordinate for flux-based interpolation methods - if ((interp_method == CGRID_VELOCITY) || (interp_method == BGRID_W_VELOCITY)) { - if ((gridindexingtype == NEMO) || (gridindexingtype == MOM5) || (gridindexingtype == POP)) { - // velocity is on the northeast of a tracer cell - xsi = 1; - eta = 1; - } else if ((gridindexingtype == MITGCM) || (gridindexingtype == CROCO)) { - // velocity is on the southwest of a tracer cell - xsi = 0; - eta = 0; - } - } else if (interp_method == BGRID_VELOCITY) { - if (gridindexingtype == MOM5) { - zeta = 1; - } else { - zeta = 0; - } - } - if ((gridindexingtype == MOM5) && (zi[igrid] == -1)) { - INTERP(spatial_interpolation_bilinear, spatial_interpolation_trilinear_surface); - } else if ((gridindexingtype == POP) && (zi[igrid] == grid->zdim-2)) { - INTERP(spatial_interpolation_bilinear, spatial_interpolation_trilinear_bottom); - } else { - INTERP(spatial_interpolation_bilinear, spatial_interpolation_trilinear); - } - } else if (interp_method == NEAREST) { - INTERP(spatial_interpolation_nearest2D, spatial_interpolation_nearest3D); - } else if ((interp_method == CGRID_TRACER) || (interp_method == BGRID_TRACER)) { - if ((gridindexingtype == POP) && (zi[igrid] == grid->zdim-2)) { - INTERP(spatial_interpolation_tracer_bc_grid_2D, spatial_interpolation_tracer_bc_grid_bottom); - } else { - INTERP(spatial_interpolation_tracer_bc_grid_2D, spatial_interpolation_tracer_bc_grid_3D); - } - } else if (interp_method == LINEAR_INVDIST_LAND_TRACER) { - INTERP(spatial_interpolation_bilinear_invdist_land, spatial_interpolation_trilinear_invdist_land); - } else { - return ERROR; - } - - // tsrch = t0 in the case where val[1] isn't populated, so this - // gives the right interpolation in either case - *value = val[0] + (val[1] - val[0]) * (float)((tsrch - t0) / (t1 - t0)); - - return SUCCESS; -#undef INTERP -} - -static double dist(double lat1, double lat2, double lon1, double lon2, int sphere_mesh, double lat) -{ - if (sphere_mesh == 1){ - double rad = M_PI / 180.; - double deg2m = 1852 * 60.; - return sqrt((lon2-lon1)*(lon2-lon1) * deg2m * deg2m * cos(rad * lat) * cos(rad * lat) + (lat2-lat1)*(lat2-lat1) * deg2m * deg2m); - } - else{ - return sqrt((lon2-lon1)*(lon2-lon1) + (lat2-lat1)*(lat2-lat1)); - } -} - -/* Linear interpolation routine for 2D C grid */ -static inline StatusCode spatial_interpolation_UV_c_grid(double eta, double xsi, - int yi, int xi, - CStructuredGrid *grid, GridType gtype, - float dataU[2][2], float dataV[2][2], - float *u, float *v) -{ - /* Cast data array into data[lat][lon] as per NEMO convention */ - int xdim = grid->xdim; - - double xgrid_loc[4]; - double ygrid_loc[4]; - int iN; - if( (gtype == RECTILINEAR_Z_GRID) || (gtype == RECTILINEAR_S_GRID) ){ - float *xgrid = grid->lon; - float *ygrid = grid->lat; - for (iN=0; iN < 4; ++iN){ - xgrid_loc[iN] = xgrid[xi+min(1, (iN%3))]; - ygrid_loc[iN] = ygrid[yi+iN/2]; - } - } - else{ - float (* xgrid)[xdim] = (float (*)[xdim]) grid->lon; - float (* ygrid)[xdim] = (float (*)[xdim]) grid->lat; - for (iN=0; iN < 4; ++iN){ - xgrid_loc[iN] = xgrid[yi+iN/2][xi+min(1, (iN%3))]; - ygrid_loc[iN] = ygrid[yi+iN/2][xi+min(1, (iN%3))]; - } - } - int i4; - for (i4 = 1; i4 < 4; ++i4){ - if (xgrid_loc[i4] < xgrid_loc[0] - 180) xgrid_loc[i4] += 360; - if (xgrid_loc[i4] > xgrid_loc[0] + 180) xgrid_loc[i4] -= 360; - } - - - double phi[4]; - phi2D_lin(eta, 0.,phi); - double U0 = dataU[1][0] * dist(ygrid_loc[3], ygrid_loc[0], xgrid_loc[3], xgrid_loc[0], grid->sphere_mesh, dot_prod(phi, ygrid_loc, 4)); - phi2D_lin(eta, 1., phi); - double U1 = dataU[1][1] * dist(ygrid_loc[1], ygrid_loc[2], xgrid_loc[1], xgrid_loc[2], grid->sphere_mesh, dot_prod(phi, ygrid_loc, 4)); - phi2D_lin(0., xsi, phi); - double V0 = dataV[0][1] * dist(ygrid_loc[0], ygrid_loc[1], xgrid_loc[0], xgrid_loc[1], grid->sphere_mesh, dot_prod(phi, ygrid_loc, 4)); - phi2D_lin(1., xsi, phi); - double V1 = dataV[1][1] * dist(ygrid_loc[2], ygrid_loc[3], xgrid_loc[2], xgrid_loc[3], grid->sphere_mesh, dot_prod(phi, ygrid_loc, 4)); - double U = (1-xsi) * U0 + xsi * U1; - double V = (1-eta) * V0 + eta * V1; - - double dphidxsi[4] = {eta-1, 1-eta, eta, -eta}; - double dphideta[4] = {xsi-1, -xsi, xsi, 1-xsi}; - double dxdxsi = 0; double dxdeta = 0; - double dydxsi = 0; double dydeta = 0; - int i; - for(i=0; i<4; ++i){ - dxdxsi += xgrid_loc[i] *dphidxsi[i]; - dxdeta += xgrid_loc[i] *dphideta[i]; - dydxsi += ygrid_loc[i] *dphidxsi[i]; - dydeta += ygrid_loc[i] *dphideta[i]; - } - double meshJac = 1; - if (grid->sphere_mesh == 1){ - double deg2m = 1852 * 60.; - double rad = M_PI / 180.; - phi2D_lin(eta, xsi, phi); - double lat = dot_prod(phi, ygrid_loc, 4); - meshJac = deg2m * deg2m * cos(rad * lat); - } - double jac = (dxdxsi*dydeta - dxdeta * dydxsi) * meshJac; - - *u = ( (-(1-eta) * U - (1-xsi) * V ) * xgrid_loc[0] + - ( (1-eta) * U - xsi * V ) * xgrid_loc[1] + - ( eta * U + xsi * V ) * xgrid_loc[2] + - ( -eta * U + (1-xsi) * V ) * xgrid_loc[3] ) / jac; - *v = ( (-(1-eta) * U - (1-xsi) * V ) * ygrid_loc[0] + - ( (1-eta) * U - xsi * V ) * ygrid_loc[1] + - ( eta * U + xsi * V ) * ygrid_loc[2] + - ( -eta * U + (1-xsi) * V ) * ygrid_loc[3] ) / jac; - - return SUCCESS; -} - - - -static inline StatusCode temporal_interpolationUV_c_grid(double time, type_coord z, type_coord y, type_coord x, - CField *U, CField *V, - GridType gtype, int *ti, int *zi, int *yi, int *xi, - float *u, float *v, int gridindexingtype) -{ - StatusCode status; - CStructuredGrid *grid = U->grid->grid; - int igrid = U->igrid; - - /* Find time index for temporal interpolation */ - if (U->time_periodic == 0 && U->allow_time_extrapolation == 0 && (time < grid->time[0] || time > grid->time[grid->tdim-1])){ - return ERRORTIMEEXTRAPOLATION; - } - status = search_time_index(&time, grid->tdim, grid->time, &ti[igrid], U->time_periodic, grid->tfull_min, grid->tfull_max, grid->periods); CHECKSTATUS(status); - - double xsi, eta, zeta; - - - if (ti[igrid] < grid->tdim-1 && time > grid->time[ti[igrid]]) { - float u0, u1, v0, v1; - double t0 = grid->time[ti[igrid]]; double t1 = grid->time[ti[igrid]+1]; - /* Identify grid cell to sample through local linear search */ - status = search_indices(time, z, y, x, grid, ti[igrid], &zi[igrid], &yi[igrid], &xi[igrid], &zeta, &eta, &xsi, gtype, t0, t1, CGRID_VELOCITY, gridindexingtype); CHECKSTATUS(status); - if (grid->zdim==1){ - float data2D_U[2][2][2], data2D_V[2][2][2]; - if (gridindexingtype == NEMO) { - status = getCell2D(U, ti[igrid], yi[igrid], xi[igrid], data2D_U, 0); CHECKSTATUS(status); - status = getCell2D(V, ti[igrid], yi[igrid], xi[igrid], data2D_V, 0); CHECKSTATUS(status); - } - else if ((gridindexingtype == MITGCM) || (gridindexingtype == CROCO)) { - status = getCell2D(U, ti[igrid], yi[igrid]-1, xi[igrid], data2D_U, 0); CHECKSTATUS(status); - status = getCell2D(V, ti[igrid], yi[igrid], xi[igrid]-1, data2D_V, 0); CHECKSTATUS(status); - } - status = spatial_interpolation_UV_c_grid(eta, xsi, yi[igrid], xi[igrid], grid, gtype, data2D_U[0], data2D_V[0], &u0, &v0); CHECKSTATUS(status); - status = spatial_interpolation_UV_c_grid(eta, xsi, yi[igrid], xi[igrid], grid, gtype, data2D_U[1], data2D_V[1], &u1, &v1); CHECKSTATUS(status); - - } else { - float data3D_U[2][2][2][2], data3D_V[2][2][2][2]; - if (gridindexingtype == NEMO) { - status = getCell3D(U, ti[igrid], zi[igrid], yi[igrid], xi[igrid], data3D_U, 0); CHECKSTATUS(status); - status = getCell3D(V, ti[igrid], zi[igrid], yi[igrid], xi[igrid], data3D_V, 0); CHECKSTATUS(status); - } - else if ((gridindexingtype == MITGCM) || (gridindexingtype == CROCO)) { - status = getCell3D(U, ti[igrid], zi[igrid], yi[igrid]-1, xi[igrid], data3D_U, 0); CHECKSTATUS(status); - status = getCell3D(V, ti[igrid], zi[igrid], yi[igrid], xi[igrid]-1, data3D_V, 0); CHECKSTATUS(status); - } - status = spatial_interpolation_UV_c_grid(eta, xsi, yi[igrid], xi[igrid], grid, gtype, data3D_U[0][0], data3D_V[0][0], &u0, &v0); CHECKSTATUS(status); - status = spatial_interpolation_UV_c_grid(eta, xsi, yi[igrid], xi[igrid], grid, gtype, data3D_U[1][0], data3D_V[1][0], &u1, &v1); CHECKSTATUS(status); - } - *u = u0 + (u1 - u0) * (float)((time - t0) / (t1 - t0)); - *v = v0 + (v1 - v0) * (float)((time - t0) / (t1 - t0)); - return SUCCESS; - } else { - double t0 = grid->time[ti[igrid]]; - status = search_indices(t0, z, y, x, grid, ti[igrid], &zi[igrid], &yi[igrid], &xi[igrid], &zeta, &eta, &xsi, gtype, t0, t0+1, CGRID_VELOCITY, gridindexingtype); CHECKSTATUS(status); - if (grid->zdim==1){ - float data2D_U[2][2][2], data2D_V[2][2][2]; - if (gridindexingtype == NEMO) { - status = getCell2D(U, ti[igrid], yi[igrid], xi[igrid], data2D_U, 1); CHECKSTATUS(status); - status = getCell2D(V, ti[igrid], yi[igrid], xi[igrid], data2D_V, 1); CHECKSTATUS(status); - } - else if ((gridindexingtype == MITGCM) || (gridindexingtype == CROCO)) { - status = getCell2D(U, ti[igrid], yi[igrid]-1, xi[igrid], data2D_U, 1); CHECKSTATUS(status); - status = getCell2D(V, ti[igrid], yi[igrid], xi[igrid]-1, data2D_V, 1); CHECKSTATUS(status); - } - status = spatial_interpolation_UV_c_grid(eta, xsi, yi[igrid], xi[igrid], grid, gtype, data2D_U[0], data2D_V[0], u, v); CHECKSTATUS(status); - } - else{ - float data3D_U[2][2][2][2], data3D_V[2][2][2][2]; - if (gridindexingtype == NEMO) { - status = getCell3D(U, ti[igrid], zi[igrid], yi[igrid], xi[igrid], data3D_U, 1); CHECKSTATUS(status); - status = getCell3D(V, ti[igrid], zi[igrid], yi[igrid], xi[igrid], data3D_V, 1); CHECKSTATUS(status); - } - else if ((gridindexingtype == MITGCM) || (gridindexingtype == CROCO)) { - status = getCell3D(U, ti[igrid], zi[igrid], yi[igrid]-1, xi[igrid], data3D_U, 1); CHECKSTATUS(status); - status = getCell3D(V, ti[igrid], zi[igrid], yi[igrid], xi[igrid]-1, data3D_V, 1); CHECKSTATUS(status); - } - status = spatial_interpolation_UV_c_grid(eta, xsi, yi[igrid], xi[igrid], grid, gtype, data3D_U[0][0], data3D_V[0][0], u, v); CHECKSTATUS(status); - } - return SUCCESS; - } -} - -/* Quadratic interpolation routine for 3D C grid */ -static inline StatusCode spatial_interpolation_UVW_c_grid(double zeta, double eta, double xsi, - int ti, int zi, int yi, int xi, - CStructuredGrid *grid, GridType gtype, - float dataU[2][2][2], float dataV[2][2][2], float dataW[2][2][2], - float *u, float *v, float *w, int gridindexingtype) -{ - /* Cast data array into data[lat][lon] as per NEMO convention */ - int xdim = grid->xdim; - int ydim = grid->ydim; - int zdim = grid->zdim; - - float xgrid_loc[4]; - float ygrid_loc[4]; - int iN; - if( gtype == RECTILINEAR_S_GRID ){ - float *xgrid = grid->lon; - float *ygrid = grid->lat; - for (iN=0; iN < 4; ++iN){ - xgrid_loc[iN] = xgrid[xi+min(1, (iN%3))]; - ygrid_loc[iN] = ygrid[yi+iN/2]; - } - } - else{ - float (* xgrid)[xdim] = (float (*)[xdim]) grid->lon; - float (* ygrid)[xdim] = (float (*)[xdim]) grid->lat; - for (iN=0; iN < 4; ++iN){ - xgrid_loc[iN] = xgrid[yi+iN/2][xi+min(1, (iN%3))]; - ygrid_loc[iN] = ygrid[yi+iN/2][xi+min(1, (iN%3))]; - } - } - int i4; - for (i4 = 1; i4 < 4; ++i4){ - if (xgrid_loc[i4] < xgrid_loc[0] - 180) xgrid_loc[i4] += 360; - if (xgrid_loc[i4] > xgrid_loc[0] + 180) xgrid_loc[i4] -= 360; - } - - float u0 = dataU[0][1][0]; - float u1 = dataU[0][1][1]; - float v0 = dataV[0][0][1]; - float v1 = dataV[0][1][1]; - float w0 = dataW[0][1][1]; - float w1 = dataW[1][1][1]; - - double px[8] = {xgrid_loc[0], xgrid_loc[1], xgrid_loc[2], xgrid_loc[3], - xgrid_loc[0], xgrid_loc[1], xgrid_loc[2], xgrid_loc[3]}; - double py[8] = {ygrid_loc[0], ygrid_loc[1], ygrid_loc[2], ygrid_loc[3], - ygrid_loc[0], ygrid_loc[1], ygrid_loc[2], ygrid_loc[3]}; - double pz[8]; - if (grid->z4d == 1){ - float (*zvals)[zdim][ydim][xdim] = (float (*)[zdim][ydim][xdim]) grid->depth; - for (iN=0; iN < 4; ++iN){ - pz[iN] = zvals[ti][zi][yi+iN/2][xi+min(1, (iN%3))]; - pz[iN+4] = zvals[ti][zi+1][yi+iN/2][xi+min(1, (iN%3))]; - } - } - else{ - float (*zvals)[ydim][xdim] = (float (*)[ydim][xdim]) grid->depth; - for (iN=0; iN < 4; ++iN){ - pz[iN] = zvals[zi][yi+iN/2][xi+min(1, (iN%3))]; - pz[iN+4] = zvals[zi+1][yi+iN/2][xi+min(1, (iN%3))]; - } - } - - double U0 = u0 * jacobian3D_lin_face(pz, py, px, zeta, eta, 0, ZONAL, grid->sphere_mesh); - double U1 = u1 * jacobian3D_lin_face(pz, py, px, zeta, eta, 1, ZONAL, grid->sphere_mesh); - double V0 = v0 * jacobian3D_lin_face(pz, py, px, zeta, 0, xsi, MERIDIONAL, grid->sphere_mesh); - double V1 = v1 * jacobian3D_lin_face(pz, py, px, zeta, 1, xsi, MERIDIONAL, grid->sphere_mesh); - double W0 = w0 * jacobian3D_lin_face(pz, py, px, 0, eta, xsi, VERTICAL, grid->sphere_mesh); - double W1 = w1 * jacobian3D_lin_face(pz, py, px, 1, eta, xsi, VERTICAL, grid->sphere_mesh); - - // Computing fluxes in half left hexahedron -> flux_u05 - double xxu[8] = {px[0], (px[0]+px[1])/2, (px[2]+px[3])/2, px[3], px[4], (px[4]+px[5])/2, (px[6]+px[7])/2, px[7]}; - double yyu[8] = {py[0], (py[0]+py[1])/2, (py[2]+py[3])/2, py[3], py[4], (py[4]+py[5])/2, (py[6]+py[7])/2, py[7]}; - double zzu[8] = {pz[0], (pz[0]+pz[1])/2, (pz[2]+pz[3])/2, pz[3], pz[4], (pz[4]+pz[5])/2, (pz[6]+pz[7])/2, pz[7]}; - double flux_u0 = u0 * jacobian3D_lin_face(zzu, yyu, xxu, .5, .5, 0, ZONAL, grid->sphere_mesh); - double flux_v0_halfx = v0 * jacobian3D_lin_face(zzu, yyu, xxu, .5, 0, .5, MERIDIONAL, grid->sphere_mesh); - double flux_v1_halfx = v1 * jacobian3D_lin_face(zzu, yyu, xxu, .5, 1, .5, MERIDIONAL, grid->sphere_mesh); - double flux_w0_halfx = w0 * jacobian3D_lin_face(zzu, yyu, xxu, 0, .5, .5, VERTICAL, grid->sphere_mesh); - double flux_w1_halfx = w1 * jacobian3D_lin_face(zzu, yyu, xxu, 1, .5, .5, VERTICAL, grid->sphere_mesh); - double flux_u05 = flux_u0 + flux_v0_halfx - flux_v1_halfx + flux_w0_halfx - flux_w1_halfx; - - // Computing fluxes in half front hexahedron -> flux_v05 - double xxv[8] = {px[0], px[1], (px[1]+px[2])/2, (px[0]+px[3])/2, px[4], px[5], (px[5]+px[6])/2, (px[4]+px[7])/2}; - double yyv[8] = {py[0], py[1], (py[1]+py[2])/2, (py[0]+py[3])/2, py[4], py[5], (py[5]+py[6])/2, (py[4]+py[7])/2}; - double zzv[8] = {pz[0], pz[1], (pz[1]+pz[2])/2, (pz[0]+pz[3])/2, pz[4], pz[5], (pz[5]+pz[6])/2, (pz[4]+pz[7])/2}; - double flux_u0_halfy = u0 * jacobian3D_lin_face(zzv, yyv, xxv, .5, .5, 0, ZONAL, grid->sphere_mesh); - double flux_u1_halfy = u1 * jacobian3D_lin_face(zzv, yyv, xxv, .5, .5, 1, ZONAL, grid->sphere_mesh); - double flux_v0 = v0 * jacobian3D_lin_face(zzv, yyv, xxv, .5, 0, .5, MERIDIONAL, grid->sphere_mesh); - double flux_w0_halfy = w0 * jacobian3D_lin_face(zzv, yyv, xxv, 0, .5, .5, VERTICAL, grid->sphere_mesh); - double flux_w1_halfy = w1 * jacobian3D_lin_face(zzv, yyv, xxv, 1, .5, .5, VERTICAL, grid->sphere_mesh); - double flux_v05 = flux_u0_halfy - flux_u1_halfy + flux_v0 + flux_w0_halfy - flux_w1_halfy; - - // Computing fluxes in half lower hexahedron -> flux_w05 - double xx[8] = {px[0], px[1], px[2], px[3], (px[0]+px[4])/2, (px[1]+px[5])/2, (px[2]+px[6])/2, (px[3]+px[7])/2}; - double yy[8] = {py[0], py[1], py[2], py[3], (py[0]+py[4])/2, (py[1]+py[5])/2, (py[2]+py[6])/2, (py[3]+py[7])/2}; - double zz[8] = {pz[0], pz[1], pz[2], pz[3], (pz[0]+pz[4])/2, (pz[1]+pz[5])/2, (pz[2]+pz[6])/2, (pz[3]+pz[7])/2}; - double flux_u0_halfz = u0 * jacobian3D_lin_face(zz, yy, xx, .5, .5, 0, ZONAL, grid->sphere_mesh); - double flux_u1_halfz = u1 * jacobian3D_lin_face(zz, yy, xx, .5, .5, 1, ZONAL, grid->sphere_mesh); - double flux_v0_halfz = v0 * jacobian3D_lin_face(zz, yy, xx, .5, 0, .5, MERIDIONAL, grid->sphere_mesh); - double flux_v1_halfz = v1 * jacobian3D_lin_face(zz, yy, xx, .5, 1, .5, MERIDIONAL, grid->sphere_mesh); - double flux_w0 = w0 * jacobian3D_lin_face(zz, yy, xx, 0, .5, .5, VERTICAL, grid->sphere_mesh); - double flux_w05 = flux_u0_halfz - flux_u1_halfz + flux_v0_halfz - flux_v1_halfz + flux_w0; - - double surf_u05 = jacobian3D_lin_face(pz, py, px, .5, .5, .5, ZONAL, grid->sphere_mesh); - double jac_u05 = jacobian3D_lin_face(pz, py, px, zeta, eta, .5, ZONAL, grid->sphere_mesh); - double U05 = flux_u05 / surf_u05 * jac_u05; - - double surf_v05 = jacobian3D_lin_face(pz, py, px, .5, .5, .5, MERIDIONAL, grid->sphere_mesh); - double jac_v05 = jacobian3D_lin_face(pz, py, px, zeta, .5, xsi, MERIDIONAL, grid->sphere_mesh); - double V05 = flux_v05 / surf_v05 * jac_v05; - - double surf_w05 = jacobian3D_lin_face(pz, py, px, .5, .5, .5, VERTICAL, grid->sphere_mesh); - double jac_w05 = jacobian3D_lin_face(pz, py, px, .5, eta, xsi, VERTICAL, grid->sphere_mesh); - double W05 = flux_w05 / surf_w05 * jac_w05; - - double jac = jacobian3D_lin(pz, py, px, zeta, eta, xsi, grid->sphere_mesh); - - double phi[3]; - phi1D_quad(xsi, phi); - double uvec[3] = {U0, U05, U1}; - double dxsidt = dot_prod(phi, uvec, 3) / jac; - phi1D_quad(eta, phi); - double vvec[3] = {V0, V05, V1}; - double detadt = dot_prod(phi, vvec, 3) / jac; - phi1D_quad(zeta, phi); - double wvec[3] = {W0, W05, W1}; - double dzetdt = dot_prod(phi, wvec, 3) / jac; - - double dphidxsi[8], dphideta[8], dphidzeta[8]; - dphidxsi3D_lin(zeta, eta, xsi, dphidzeta, dphideta, dphidxsi); - - *u = dot_prod(dphidxsi, px, 8) * dxsidt + dot_prod(dphideta, px, 8) * detadt + dot_prod(dphidzeta, px, 8) * dzetdt; - *v = dot_prod(dphidxsi, py, 8) * dxsidt + dot_prod(dphideta, py, 8) * detadt + dot_prod(dphidzeta, py, 8) * dzetdt; - *w = dot_prod(dphidxsi, pz, 8) * dxsidt + dot_prod(dphideta, pz, 8) * detadt + dot_prod(dphidzeta, pz, 8) * dzetdt; - - return SUCCESS; -} - -static inline StatusCode temporal_interpolationUVW_c_grid(double time, type_coord z, type_coord y, type_coord x, - CField *U, CField *V, CField *W, - GridType gtype, int *ti, int *zi, int *yi, int *xi, - float *u, float *v, float *w, int gridindexingtype) -{ - StatusCode status; - CStructuredGrid *grid = U->grid->grid; - int igrid = U->igrid; - - /* Find time index for temporal interpolation */ - if (U->time_periodic == 0 && U->allow_time_extrapolation == 0 && (time < grid->time[0] || time > grid->time[grid->tdim-1])){ - return ERRORTIMEEXTRAPOLATION; - } - status = search_time_index(&time, grid->tdim, grid->time, &ti[igrid], U->time_periodic, grid->tfull_min, grid->tfull_max, grid->periods); CHECKSTATUS(status); - - double xsi, eta, zeta; - float data3D_U[2][2][2][2]; - float data3D_V[2][2][2][2]; - float data3D_W[2][2][2][2]; - - - if (ti[igrid] < grid->tdim-1 && time > grid->time[ti[igrid]]) { - float u0, u1, v0, v1, w0, w1; - double t0 = grid->time[ti[igrid]]; double t1 = grid->time[ti[igrid]+1]; - /* Identify grid cell to sample through local linear search */ - status = search_indices(time, z, y, x, grid, ti[igrid], &zi[igrid], &yi[igrid], &xi[igrid], &zeta, &eta, &xsi, gtype, t0, t1, CGRID_VELOCITY, gridindexingtype); CHECKSTATUS(status); - status = getCell3D(U, ti[igrid], zi[igrid], yi[igrid], xi[igrid], data3D_U, 0); CHECKSTATUS(status); - status = getCell3D(V, ti[igrid], zi[igrid], yi[igrid], xi[igrid], data3D_V, 0); CHECKSTATUS(status); - status = getCell3D(W, ti[igrid], zi[igrid], yi[igrid], xi[igrid], data3D_W, 0); CHECKSTATUS(status); - if (grid->zdim==1){ - return ERROR; - } else { - status = spatial_interpolation_UVW_c_grid(zeta, eta, xsi, ti[igrid], zi[igrid], yi[igrid], xi[igrid], grid, gtype, data3D_U[0], data3D_V[0], data3D_W[0], &u0, &v0, &w0, gridindexingtype); CHECKSTATUS(status); - status = spatial_interpolation_UVW_c_grid(zeta, eta, xsi, ti[igrid]+1, zi[igrid], yi[igrid], xi[igrid], grid, gtype, data3D_U[1], data3D_V[1], data3D_W[1], &u1, &v1, &w1, gridindexingtype); CHECKSTATUS(status); - } - *u = u0 + (u1 - u0) * (float)((time - t0) / (t1 - t0)); - *v = v0 + (v1 - v0) * (float)((time - t0) / (t1 - t0)); - *w = w0 + (w1 - w0) * (float)((time - t0) / (t1 - t0)); - return SUCCESS; - } else { - double t0 = grid->time[ti[igrid]]; - status = search_indices(t0, z, y, x, grid, ti[igrid], &zi[igrid], &yi[igrid], &xi[igrid], &zeta, &eta, &xsi, gtype, t0, t0+1, CGRID_VELOCITY, gridindexingtype); CHECKSTATUS(status); - status = getCell3D(U, ti[igrid], zi[igrid], yi[igrid], xi[igrid], data3D_U, 1); CHECKSTATUS(status); - status = getCell3D(V, ti[igrid], zi[igrid], yi[igrid], xi[igrid], data3D_V, 1); CHECKSTATUS(status); - status = getCell3D(W, ti[igrid], zi[igrid], yi[igrid], xi[igrid], data3D_W, 1); CHECKSTATUS(status); - if (grid->zdim==1){ - return ERROR; - } - else{ - status = spatial_interpolation_UVW_c_grid(zeta, eta, xsi, ti[igrid], zi[igrid], yi[igrid], xi[igrid], grid, gtype, data3D_U[0], data3D_V[0], data3D_W[0], u, v, w, gridindexingtype); CHECKSTATUS(status); - } - return SUCCESS; - } -} - -static inline StatusCode calculate_slip_conditions_2D(double eta, double xsi, - float dataU[2][2], float dataV[2][2], float dataW[2][2], - float *u, float *v, float *w, int interp_method, int withW) -{ - float f_u = 1, f_v = 1, f_w = 1; - if ((is_zero_flt(dataU[0][0])) && (is_zero_flt(dataU[0][1])) && - (is_zero_flt(dataV[0][0])) && (is_zero_flt(dataV[0][1])) && eta > 0.){ - if (interp_method == PARTIALSLIP) { - f_u = f_u * (.5 + .5 * eta) / eta; - if (withW) { - f_w = f_w * (.5 + .5 * eta) / eta; - } - } else if (interp_method == FREESLIP) { - f_u = f_u / eta; - if (withW) { - f_w = f_w / eta; - } - } - } - if ((is_zero_flt(dataU[1][0])) && (is_zero_flt(dataU[1][1])) && - (is_zero_flt(dataV[1][0])) && (is_zero_flt(dataV[1][1])) && eta < 1.){ - if (interp_method == PARTIALSLIP) { - f_u = f_u * (1 - .5 * eta) / (1 - eta); - if (withW) { - f_w = f_w * (1 - .5 * eta) / (1 - eta); - } - } else if (interp_method == FREESLIP) { - f_u = f_u / (1 - eta); - if (withW) { - f_w = f_w / (1 - eta); - } - } - } - if ((is_zero_flt(dataU[0][0])) && (is_zero_flt(dataU[1][0])) && - (is_zero_flt(dataV[0][0])) && (is_zero_flt(dataV[1][0])) && xsi > 0.){ - if (interp_method == PARTIALSLIP) { - f_v = f_v * (.5 + .5 * xsi) / xsi; - if (withW) { - f_w = f_w * (.5 + .5 * xsi) / xsi; - } - } else if (interp_method == FREESLIP) { - f_v = f_v / xsi; - if (withW) { - f_w = f_w / xsi; - } - } - } - if ((is_zero_flt(dataU[0][1])) && (is_zero_flt(dataU[1][1])) && - (is_zero_flt(dataV[0][1])) && (is_zero_flt(dataV[1][1])) && xsi < 1.){ - if (interp_method == PARTIALSLIP) { - f_v = f_v * (1 - .5 * xsi) / (1 - xsi); - if (withW) { - f_w = f_w * (1 - .5 * xsi) / (1 - xsi); - } - } else if (interp_method == FREESLIP) { - f_v = f_v / (1 - xsi); - if (withW) { - f_w = f_w / (1 - xsi); - } - } - } - *u *= f_u; - *v *= f_v; - if (withW) { - *w *= f_w; - } - - return SUCCESS; -} - -static inline StatusCode calculate_slip_conditions_3D(double zeta, double eta, double xsi, - float dataU[2][2][2], float dataV[2][2][2], float dataW[2][2][2], - float *u, float *v, float *w, int interp_method, int withW) -{ - float f_u = 1, f_v = 1, f_w = 1; - if ((is_zero_flt(dataU[0][0][0])) && (is_zero_flt(dataU[0][0][1])) && (is_zero_flt(dataU[1][0][0])) && (is_zero_flt(dataU[1][0][1])) && - (is_zero_flt(dataV[0][0][0])) && (is_zero_flt(dataV[0][0][1])) && (is_zero_flt(dataV[1][0][0])) && (is_zero_flt(dataV[1][0][1])) && - eta > 0.){ - if (interp_method == PARTIALSLIP) { - f_u = f_u * (.5 + .5 * eta) / eta; - if (withW) { - f_w = f_w * (.5 + .5 * eta) / eta; - } - } else if (interp_method == FREESLIP) { - f_u = f_u / eta; - if (withW) { - f_w = f_w / eta; - } - } - } - if ((is_zero_flt(dataU[0][1][0])) && (is_zero_flt(dataU[0][1][1])) && (is_zero_flt(dataU[1][1][0])) && (is_zero_flt(dataU[1][1][1])) && - (is_zero_flt(dataV[0][1][0])) && (is_zero_flt(dataV[0][1][1])) && (is_zero_flt(dataV[1][1][0])) && (is_zero_flt(dataV[1][1][1])) && - eta < 1.){ - if (interp_method == PARTIALSLIP) { - f_u = f_u * (1 - .5 * eta) / (1 - eta); - if (withW) { - f_w = f_w * (1 - .5 * eta) / (1 - eta); - } - } else if (interp_method == FREESLIP) { - f_u = f_u / (1 - eta); - if (withW) { - f_w = f_w / (1 - eta); - } - } - } - if ((is_zero_flt(dataU[0][0][0])) && (is_zero_flt(dataU[0][1][0])) && (is_zero_flt(dataU[1][0][0])) && (is_zero_flt(dataU[1][1][0])) && - (is_zero_flt(dataV[0][0][0])) && (is_zero_flt(dataV[0][1][0])) && (is_zero_flt(dataV[1][0][0])) && (is_zero_flt(dataV[1][1][0])) && - xsi > 0.){ - if (interp_method == PARTIALSLIP) { - f_v = f_v * (.5 + .5 * xsi) / xsi; - if (withW) { - f_w = f_w * (.5 + .5 * xsi) / xsi; - } - } else if (interp_method == FREESLIP) { - f_v = f_v / xsi; - if (withW) { - f_w = f_w / xsi; - } - } - } - if ((is_zero_flt(dataU[0][0][1])) && (is_zero_flt(dataU[0][1][1])) && (is_zero_flt(dataU[1][0][1])) && (is_zero_flt(dataU[1][1][1])) && - (is_zero_flt(dataV[0][0][1])) && (is_zero_flt(dataV[0][1][1])) && (is_zero_flt(dataV[1][0][1])) && (is_zero_flt(dataV[1][1][1])) && - xsi < 1.){ - if (interp_method == PARTIALSLIP) { - f_v = f_v * (1 - .5 * xsi) / (1 - xsi); - if (withW) { - f_w = f_w * (1 - .5 * xsi) / (1 - xsi); - } - } else if (interp_method == FREESLIP) { - f_v = f_v / (1 - xsi); - if (withW) { - f_w = f_w / (1 - xsi); - } - } - } - if ((is_zero_flt(dataU[0][0][0])) && (is_zero_flt(dataU[0][0][1])) && (is_zero_flt(dataU[0][1][0])) && (is_zero_flt(dataU[0][1][1])) && - (is_zero_flt(dataV[0][0][0])) && (is_zero_flt(dataV[0][0][1])) && (is_zero_flt(dataV[0][1][0])) && (is_zero_flt(dataV[0][1][1])) && - zeta > 0.){ - if (interp_method == PARTIALSLIP) { - f_u = f_u * (.5 + .5 * zeta) / zeta; - f_v = f_v * (.5 + .5 * zeta) / zeta; - } else if (interp_method == FREESLIP) { - f_u = f_u / zeta; - f_v = f_v / zeta; - } - } - if ((is_zero_flt(dataU[1][0][0])) && (is_zero_flt(dataU[1][0][1])) && (is_zero_flt(dataU[1][1][0])) && (is_zero_flt(dataU[1][1][1])) && - (is_zero_flt(dataV[1][0][0])) && (is_zero_flt(dataV[1][0][1])) && (is_zero_flt(dataV[1][1][0])) && (is_zero_flt(dataV[1][1][1])) && - zeta < 1.){ - if (interp_method == PARTIALSLIP) { - f_u = f_u * (1 - .5 * zeta) / (1 - zeta); - f_v = f_v * (1 - .5 * zeta) / (1 - zeta); - } else if (interp_method == FREESLIP) { - f_u = f_u / (1 - zeta); - f_v = f_v / (1 - zeta); - } - } - *u *= f_u; - *v *= f_v; - if (withW) { - *w *= f_w; - } - - return SUCCESS; -} - -static inline StatusCode temporal_interpolation_slip(double time, type_coord z, type_coord y, type_coord x, - CField *U, CField *V, CField *W, - GridType gtype, int *ti, int *zi, int *yi, int *xi, - float *u, float *v, float *w, int interp_method, int gridindexingtype, int withW) -{ - StatusCode status; - CStructuredGrid *grid = U->grid->grid; - int igrid = U->igrid; - - /* Find time index for temporal interpolation */ - if (U->time_periodic == 0 && U->allow_time_extrapolation == 0 && (time < grid->time[0] || time > grid->time[grid->tdim-1])){ - return ERRORTIMEEXTRAPOLATION; - } - status = search_time_index(&time, grid->tdim, grid->time, &ti[igrid], U->time_periodic, grid->tfull_min, grid->tfull_max, grid->periods); CHECKSTATUS(status); - - double xsi, eta, zeta; - - if (ti[igrid] < grid->tdim-1 && time > grid->time[ti[igrid]]) { - float u0, u1, v0, v1, w0, w1; - double t0 = grid->time[ti[igrid]]; double t1 = grid->time[ti[igrid]+1]; - /* Identify grid cell to sample through local linear search */ - status = search_indices(time, z, y, x, grid, ti[igrid], &zi[igrid], &yi[igrid], &xi[igrid], &zeta, &eta, &xsi, gtype, t0, t1, interp_method, gridindexingtype); CHECKSTATUS(status); - if (grid->zdim==1){ - float data2D_U[2][2][2], data2D_V[2][2][2], data2D_W[2][2][2]; - status = getCell2D(U, ti[igrid], yi[igrid], xi[igrid], data2D_U, 0); CHECKSTATUS(status); - status = getCell2D(V, ti[igrid], yi[igrid], xi[igrid], data2D_V, 0); CHECKSTATUS(status); - if (withW){ - status = getCell2D(W, ti[igrid], yi[igrid], xi[igrid], data2D_W, 0); CHECKSTATUS(status); - status = spatial_interpolation_bilinear(eta, xsi, data2D_W[0], &w0); CHECKSTATUS(status); - status = spatial_interpolation_bilinear(eta, xsi, data2D_W[1], &w1); CHECKSTATUS(status); - } - - status = spatial_interpolation_bilinear(eta, xsi, data2D_U[0], &u0); CHECKSTATUS(status); - status = spatial_interpolation_bilinear(eta, xsi, data2D_V[0], &v0); CHECKSTATUS(status); - status = calculate_slip_conditions_2D(eta, xsi, data2D_U[0], data2D_V[0], data2D_W[0], &u0, &v0, &w0, interp_method, withW); CHECKSTATUS(status); - - status = spatial_interpolation_bilinear(eta, xsi, data2D_U[1], &u1); CHECKSTATUS(status); - status = spatial_interpolation_bilinear(eta, xsi, data2D_V[1], &v1); CHECKSTATUS(status); - status = calculate_slip_conditions_2D(eta, xsi, data2D_U[1], data2D_V[1], data2D_W[1], &u1, &v1, &w1, interp_method, withW); CHECKSTATUS(status); - } else { - float data3D_U[2][2][2][2], data3D_V[2][2][2][2], data3D_W[2][2][2][2]; - status = getCell3D(U, ti[igrid], zi[igrid], yi[igrid], xi[igrid], data3D_U, 0); CHECKSTATUS(status); - status = getCell3D(V, ti[igrid], zi[igrid], yi[igrid], xi[igrid], data3D_V, 0); CHECKSTATUS(status); - if (withW){ - status = getCell3D(W, ti[igrid], zi[igrid], yi[igrid], xi[igrid], data3D_W, 0); CHECKSTATUS(status); - status = spatial_interpolation_trilinear(zeta, eta, xsi, data3D_W[0], &w0); CHECKSTATUS(status); - status = spatial_interpolation_trilinear(zeta, eta, xsi, data3D_W[1], &w1); CHECKSTATUS(status); - } - status = spatial_interpolation_trilinear(zeta, eta, xsi, data3D_U[0], &u0); CHECKSTATUS(status); - status = spatial_interpolation_trilinear(zeta, eta, xsi, data3D_V[0], &v0); CHECKSTATUS(status); - status = calculate_slip_conditions_3D(zeta, eta, xsi, data3D_U[0], data3D_V[0], data3D_W[0], &u0, &v0, &w0, interp_method, withW); CHECKSTATUS(status); - - status = spatial_interpolation_trilinear(zeta, eta, xsi, data3D_U[1], &u1); CHECKSTATUS(status); - status = spatial_interpolation_trilinear(zeta, eta, xsi, data3D_V[1], &v1); CHECKSTATUS(status); - status = calculate_slip_conditions_3D(zeta, eta, xsi, data3D_U[1], data3D_V[1], data3D_W[1], &u1, &v1, &w1, interp_method, withW); CHECKSTATUS(status); - } - *u = u0 + (u1 - u0) * (float)((time - t0) / (t1 - t0)); - *v = v0 + (v1 - v0) * (float)((time - t0) / (t1 - t0)); - if (withW){ - *w = w0 + (w1 - w0) * (float)((time - t0) / (t1 - t0)); - } - - } else { - double t0 = grid->time[ti[igrid]]; - status = search_indices(t0, z, y, x, grid, ti[igrid], &zi[igrid], &yi[igrid], &xi[igrid], &zeta, &eta, &xsi, gtype, t0, t0+1, interp_method, gridindexingtype); CHECKSTATUS(status); - if (grid->zdim==1){ - float data2D_U[2][2][2], data2D_V[2][2][2], data2D_W[2][2][2]; - status = getCell2D(U, ti[igrid], yi[igrid], xi[igrid], data2D_U, 1); CHECKSTATUS(status); - status = getCell2D(V, ti[igrid], yi[igrid], xi[igrid], data2D_V, 1); CHECKSTATUS(status); - if (withW){ - status = getCell2D(W, ti[igrid], yi[igrid], xi[igrid], data2D_W, 1); CHECKSTATUS(status); - status = spatial_interpolation_bilinear(eta, xsi, data2D_W[0], w); CHECKSTATUS(status); - } - - status = spatial_interpolation_bilinear(eta, xsi, data2D_U[0], u); CHECKSTATUS(status); - status = spatial_interpolation_bilinear(eta, xsi, data2D_V[0], v); CHECKSTATUS(status); - - status = calculate_slip_conditions_2D(eta, xsi, data2D_U[0], data2D_V[0], data2D_W[0], u, v, w, interp_method, withW); CHECKSTATUS(status); - } else { - float data3D_U[2][2][2][2], data3D_V[2][2][2][2], data3D_W[2][2][2][2]; - status = getCell3D(U, ti[igrid], zi[igrid], yi[igrid], xi[igrid], data3D_U, 1); CHECKSTATUS(status); - status = getCell3D(V, ti[igrid], zi[igrid], yi[igrid], xi[igrid], data3D_V, 1); CHECKSTATUS(status); - if (withW){ - status = getCell3D(W, ti[igrid], zi[igrid], yi[igrid], xi[igrid], data3D_W, 1); CHECKSTATUS(status); - status = spatial_interpolation_trilinear(zeta, eta, xsi, data3D_W[0], w); CHECKSTATUS(status); - } - status = spatial_interpolation_trilinear(zeta, eta, xsi, data3D_U[0], u); CHECKSTATUS(status); - status = spatial_interpolation_trilinear(zeta, eta, xsi, data3D_V[0], v); CHECKSTATUS(status); - status = calculate_slip_conditions_3D(zeta, eta, xsi, data3D_U[0], data3D_V[0], data3D_W[0], u, v, w, interp_method, withW); CHECKSTATUS(status); - } - } - return SUCCESS; -} - -static inline StatusCode temporal_interpolation(double time, type_coord z, type_coord y, type_coord x, - CField *f, - int *ti, int *zi, int *yi, int *xi, - float *value, int interp_method, int gridindexingtype) -{ - CGrid *_grid = f->grid; - GridType gtype = _grid->gtype; - - if (gtype == RECTILINEAR_Z_GRID || gtype == RECTILINEAR_S_GRID || gtype == CURVILINEAR_Z_GRID || gtype == CURVILINEAR_S_GRID){ - return temporal_interpolation_structured_grid(time, z, y, x, f, gtype, ti, zi, yi, xi, value, interp_method, gridindexingtype); - } - else{ - printf("Only RECTILINEAR_Z_GRID, RECTILINEAR_S_GRID, CURVILINEAR_Z_GRID and CURVILINEAR_S_GRID grids are currently implemented\n"); - return ERROR; - } -} - -static inline StatusCode temporal_interpolationUV(double time, type_coord z, type_coord y, type_coord x, - CField *U, CField *V, - int *ti, int *zi, int *yi, int *xi, - float *valueU, float *valueV, int interp_method, int gridindexingtype) -{ - StatusCode status; - if (interp_method == CGRID_VELOCITY){ - CGrid *_grid = U->grid; - GridType gtype = _grid->gtype; - status = temporal_interpolationUV_c_grid(time, z, y, x, U, V, gtype, ti, zi, yi, xi, valueU, valueV, gridindexingtype); CHECKSTATUS(status); - return SUCCESS; - } else if ((interp_method == PARTIALSLIP) || (interp_method == FREESLIP)){ - CGrid *_grid = U->grid; - CField *W = U; - GridType gtype = _grid->gtype; - int withW = 0; - status = temporal_interpolation_slip(time, z, y, x, U, V, W, gtype, ti, zi, yi, xi, valueU, valueV, 0, interp_method, gridindexingtype, withW); CHECKSTATUS(status); - return SUCCESS; - } else { - status = temporal_interpolation(time, z, y, x, U, ti, zi, yi, xi, valueU, interp_method, gridindexingtype); CHECKSTATUS(status); - status = temporal_interpolation(time, z, y, x, V, ti, zi, yi, xi, valueV, interp_method, gridindexingtype); CHECKSTATUS(status); - return SUCCESS; - } -} - -static inline StatusCode temporal_interpolationUVW(double time, type_coord z, type_coord y, type_coord x, - CField *U, CField *V, CField *W, - int *ti, int *zi, int *yi, int *xi, - float *valueU, float *valueV, float *valueW, int interp_method, int gridindexingtype) -{ - StatusCode status; - if (interp_method == CGRID_VELOCITY){ - CGrid *_grid = U->grid; - GridType gtype = _grid->gtype; - if (gtype == RECTILINEAR_S_GRID || gtype == CURVILINEAR_S_GRID){ - status = temporal_interpolationUVW_c_grid(time, z, y, x, U, V, W, gtype, ti, zi, yi, xi, valueU, valueV, valueW, gridindexingtype); CHECKSTATUS(status); - return SUCCESS; - } - } else if ((interp_method == PARTIALSLIP) || (interp_method == FREESLIP)){ - CGrid *_grid = U->grid; - GridType gtype = _grid->gtype; - int withW = 1; - status = temporal_interpolation_slip(time, z, y, x, U, V, W, gtype, ti, zi, yi, xi, valueU, valueV, valueW, interp_method, gridindexingtype, withW); CHECKSTATUS(status); - return SUCCESS; - } - status = temporal_interpolationUV(time, z, y, x, U, V, ti, zi, yi, xi, valueU, valueV, interp_method, gridindexingtype); CHECKSTATUS(status); - if (interp_method == BGRID_VELOCITY) - interp_method = BGRID_W_VELOCITY; - if (gridindexingtype == CROCO) // Linear vertical interpolation for CROCO - interp_method = LINEAR; - status = temporal_interpolation(time, z, y, x, W, ti, zi, yi, xi, valueW, interp_method, gridindexingtype); CHECKSTATUS(status); - return SUCCESS; -} - - -static inline double croco_from_z_to_sigma(double time, type_coord z, type_coord y, type_coord x, - CField *U, CField *H, CField *Zeta, - int *ti, int *zi, int *yi, int *xi, double hc, float *cs_w) -{ - float local_h, local_zeta, z0; - int status, zii; - CStructuredGrid *grid = U->grid->grid; - float *sigma_levels = grid->depth; - int zdim = grid->zdim; - float zvec[zdim]; - status = temporal_interpolation(time, 0, y, x, H, ti, zi, yi, xi, &local_h, LINEAR, CROCO); CHECKSTATUS(status); - status = temporal_interpolation(time, 0, y, x, Zeta, ti, zi, yi, xi, &local_zeta, LINEAR, CROCO); CHECKSTATUS(status); - for (zii = 0; zii < zdim; zii++) { - z0 = hc*sigma_levels[zii] + (local_h - hc) *cs_w[zii]; - zvec[zii] = z0 + local_zeta * (1 + z0 / local_h); - } - if (z >= zvec[zdim-1]) - zii = zdim - 2; - else - for (zii = 0; zii < zdim-1; zii++) - if ((z >= zvec[zii]) && (z < zvec[zii+1])) - break; - - return sigma_levels[zii] + (z - zvec[zii]) * (sigma_levels[zii + 1] - sigma_levels[zii]) / (zvec[zii + 1] - zvec[zii]); -} - -#ifdef __cplusplus -} -#endif -#endif diff --git a/parcels/include/random.h b/parcels/include/random.h deleted file mode 100644 index 07a461cef5..0000000000 --- a/parcels/include/random.h +++ /dev/null @@ -1,111 +0,0 @@ -#ifndef _PARCELS_RANDOM_H -#define _PARCELS_RANDOM_H -#ifdef __cplusplus -extern "C" { -#endif - -/**************************************************/ - - -/**************************************************/ -/* Random number generation (RNG) functions */ -/**************************************************/ - -static inline void parcels_seed(int seed) -{ - srand(seed); -} - -static inline float parcels_random() -{ - return (float)rand()/(float)(RAND_MAX); -} - -static inline float parcels_uniform(float low, float high) -{ - return (float)rand()/(float)((float)(RAND_MAX) / (high-low)) + low; -} - -static inline int parcels_randint(int low, int high) -{ - return (rand() % (high-low)) + low; -} - -static inline float parcels_normalvariate(float loc, float scale) -/* Function to create a Gaussian random variable with mean loc and standard deviation scale */ -/* Uses Box-Muller transform, adapted from ftp://ftp.taygeta.com/pub/c/boxmuller.c */ -/* (c) Copyright 1994, Everett F. Carter Jr. Permission is granted by the author to use */ -/* this software for any application provided this copyright notice is preserved. */ -{ - float x1, x2, w, y1; - - do { - x1 = 2.0 * (float)rand()/(float)(RAND_MAX) - 1.0; - x2 = 2.0 * (float)rand()/(float)(RAND_MAX) - 1.0; - w = x1 * x1 + x2 * x2; - } while ( w >= 1.0 ); - - w = sqrt( (-2.0 * log( w ) ) / w ); - y1 = x1 * w; - return( loc + y1 * scale ); -} - -static inline float parcels_expovariate(float lamb) -//Function to create an exponentially distributed random variable -{ - float u; - u = (float)rand()/((float)(RAND_MAX) + 1.0); - return (-log(1.0-u)/lamb); -} - -static inline float parcels_vonmisesvariate(float mu, float kappa) -/* Circular data distribution. */ -/* Returns a float between 0 and 2*pi */ -/* mu is the mean angle, expressed in radians between 0 and 2*pi, and */ -/* kappa is the concentration parameter, which must be greater than or */ -/* equal to zero. If kappa is equal to zero, this distribution reduces */ -/* to a uniform random angle over the range 0 to 2*pi. */ -/* Based upon an algorithm published in: Fisher, N.I., */ -/* Statistical Analysis of Circular Data", Cambridge University Press, 1993.*/ -{ - float u1, u2, u3, r, s, z, d, f, q, theta; - - if (kappa <= 1e-6){ - return (2.0 * M_PI * (float)rand()/(float)(RAND_MAX)); - } - - s = 0.5 / kappa; - if (fabs(s) <= FLT_EPSILON * fabs(s)){ - return mu; - } - r = s + sqrt(1.0 + s * s); - - do { - u1 = (float)rand()/(float)(RAND_MAX); - z = cos(M_PI * u1); - - d = z / (r + z); - u2 = (float)rand()/(float)(RAND_MAX); - } while ( ( u2 >= (1.0 - d * d) ) && ( u2 > (1.0 - d) * exp(d) ) ); - - q = 1.0 / r; - f = (q + z) / (1.0 + q * z); - u3 = (float)rand()/(float)(RAND_MAX); - - if (u3 > 0.5){ - theta = fmod(mu + acos(f), 2.0*M_PI); - } - else { - theta = fmod(mu - acos(f), 2.0*M_PI); - } - if (theta < 0){ - theta = 2.0*M_PI+theta; - } - - return theta; -} - -#ifdef __cplusplus -} -#endif -#endif diff --git a/pyproject.toml b/pyproject.toml index 70d202689a..991a6691cf 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -135,7 +135,6 @@ known-first-party = ["parcels"] [tool.mypy] files = [ - "parcels/compilation/codegenerator.py", "parcels/_typing.py", "parcels/tools/*.py", "parcels/grid.py", diff --git a/tests/test_fieldset_sampling.py b/tests/test_fieldset_sampling.py index 2318c291f0..5cc6bda7a1 100644 --- a/tests/test_fieldset_sampling.py +++ b/tests/test_fieldset_sampling.py @@ -451,8 +451,6 @@ def test_fieldset_sample_geographic_polar(fieldset_geometric_polar): pset = ParticleSet(fieldset, pclass=pclass(), lat=lat, lon=np.zeros(npart) - 45.0) pset.execute(pset.Kernel(SampleUV), endtime=1.0, dt=1.0) - # Note: 1.e-2 is a very low rtol, so there seems to be a rather - # large sampling error for the JIT correction. assert np.allclose(pset.u, lat, rtol=1e-2) diff --git a/tests/test_kernel_execution.py b/tests/test_kernel_execution.py index 883991b979..51d0b7e8cd 100644 --- a/tests/test_kernel_execution.py +++ b/tests/test_kernel_execution.py @@ -308,17 +308,3 @@ def MoveWest(particle, fieldset, time): # pragma: no cover pset.execute(pset.Kernel(MoveEast), endtime=1.0, dt=1.0) pset.execute(pset.Kernel(MoveWest), endtime=3.0, dt=1.0) assert np.allclose(pset.lon, 0.3, rtol=1e-5) # should be 0.5 + 0.1 - 0.3 = 0.3 - - -def test_compilers(): - from parcels.compilation.codecompiler import ( - CCompiler_SS, - Clang_parameters, - MinGW_parameters, - VS_parameters, - ) - - for param_class in [Clang_parameters, MinGW_parameters, VS_parameters]: - params = param_class() # noqa - - print(CCompiler_SS())