Skip to content
Snippets Groups Projects

Compare revisions

Changes are shown as if the source revision was being merged into the target revision. Learn more about comparing revisions.

Source

Select target project
No results found

Target

Select target project
  • icon-libraries/libiconmath
1 result
Show changes
Commits on Source (8)
Hide whitespace changes
Inline Side-by-side
Showing
with 685 additions and 76 deletions
.gitignore
View file @ e219c7cc
......@@ -22,4 +22,5 @@ iconmath_Tests
# Test stage files:
/**/Testing/*
iconmath_Tests
run_tests.f90
CMakeLists.txt
View file @ e219c7cc
......@@ -14,7 +14,9 @@ cmake_minimum_required(VERSION 3.18)
project(
iconmath
VERSION 1.0.0
LANGUAGES Fortran)
LANGUAGES Fortran CXX)
set(CMAKE_CXX_STANDARD 17)
option(BUILD_SHARED_LIBS "Build shared libraries" ON)
option(BUILD_TESTING "Build tests" ON)
......@@ -23,6 +25,7 @@ option(BUILD_ICONMATH_HORIZONTAL "Build horizontal library" ON)
option(IM_ENABLE_MIXED_PRECISION "Enable mixed precision" OFF)
option(IM_ENABLE_LOOP_EXCHANGE "Enable loop exchange" OFF)
option(IM_USE_CPP_BINDINGS "Use C++ bindings" OFF)
option(IM_ENABLE_DIM_SWAP "Enable dimension swap" OFF)
option(IM_ENABLE_OPENACC "Enable OpenACC support" OFF)
option(IM_ENABLE_OPENMP "Enable OpenMP support" OFF)
......
src/support/CMakeLists.txt
View file @ e219c7cc
......@@ -13,9 +13,11 @@ add_library(
iconmath-support
mo_gridman_constants.f90
mo_lib_grid_geometry_info.f90
mo_lib_loopindices.f90
mo_lib_loopindices.cpp
mo_lib_loopindices.F90
mo_math_constants.f90
mo_math_types.f90
mo_math_utilities.cpp
mo_math_utilities.F90
mo_random_number_generators.F90)
......@@ -40,6 +42,10 @@ if(IM_ENABLE_DIM_SWAP)
target_compile_definitions(iconmath-support PRIVATE __SWAPDIM)
endif()
if(IM_USE_CPP_BINDINGS)
target_compile_definitions(iconmath-support PRIVATE __USE_CPP_BINDINGS)
endif()
if(IM_ENABLE_OPENACC)
# If _OPENACC is defined, assume that the required compiler flags are already
# provided, e.g. in CMAKE_Fortran_FLAGS:
......@@ -57,7 +63,16 @@ target_include_directories(
# Path to the Fortran modules:
$<BUILD_INTERFACE:$<$<COMPILE_LANGUAGE:Fortran>:${Fortran_MODULE_DIRECTORY}>>
$<INSTALL_INTERFACE:$<$<COMPILE_LANGUAGE:Fortran>:$<INSTALL_PREFIX>/${CMAKE_INSTALL_INCLUDEDIR}>>
)
INTERFACE
# Path to the internal C/C++ headers (for testing): Requires CMake 3.15+ for
# multiple compile languages
# https://cmake.org/cmake/help/latest/manual/cmake-generator-expressions.7.html
$<BUILD_INTERFACE:$<$<COMPILE_LANGUAGE:C,CXX>:${CMAKE_CURRENT_SOURCE_DIR}>>
PRIVATE
# Path to config.h (for C and C++ only): Requires CMake 3.15+ for multiple
# compile languages
# https://cmake.org/cmake/help/latest/manual/cmake-generator-expressions.7.html
$<BUILD_INTERFACE:$<$<COMPILE_LANGUAGE:C,CXX>:${CMAKE_CURRENT_BINARY_DIR}>>)
target_link_libraries(iconmath-support PUBLIC fortran-support::fortran-support)
......
src/support/mo_lib_loopindices.f90 → src/support/mo_lib_loopindices.F90
View file @ e219c7cc
......@@ -16,12 +16,18 @@
MODULE mo_lib_loopindices
#ifdef __USE_CPP_BINDINGS
USE, INTRINSIC :: ISO_C_BINDING
#endif
IMPLICIT NONE
PRIVATE
PUBLIC :: get_indices_c_lib, get_indices_e_lib, get_indices_v_lib
#ifndef __USE_CPP_BINDINGS
CONTAINS
!-------------------------------------------------------------------------
......@@ -121,5 +127,31 @@ CONTAINS
END SUBROUTINE get_indices_v_lib
END MODULE mo_lib_loopindices
#else
INTERFACE
SUBROUTINE get_indices_c_lib(i_startidx_in, i_endidx_in, nproma, i_blk, i_startblk, i_endblk, &
i_startidx_out, i_endidx_out) BIND(C, NAME="get_indices_c_lib")
IMPORT :: C_INT
INTEGER(C_INT), VALUE :: i_startidx_in, i_endidx_in, nproma, i_blk, i_startblk, i_endblk
INTEGER(C_INT) :: i_startidx_out, i_endidx_out
END SUBROUTINE get_indices_c_lib
SUBROUTINE get_indices_e_lib(i_startidx_in, i_endidx_in, nproma, i_blk, i_startblk, i_endblk, &
i_startidx_out, i_endidx_out) BIND(C, NAME="get_indices_e_lib")
IMPORT :: C_INT
INTEGER(C_INT), VALUE :: i_startidx_in, i_endidx_in, nproma, i_blk, i_startblk, i_endblk
INTEGER(C_INT) :: i_startidx_out, i_endidx_out
END SUBROUTINE get_indices_e_lib
SUBROUTINE get_indices_v_lib(i_startidx_in, i_endidx_in, nproma, i_blk, i_startblk, i_endblk, &
i_startidx_out, i_endidx_out) BIND(C, NAME="get_indices_v_lib")
IMPORT :: C_INT
INTEGER(C_INT), VALUE :: i_startidx_in, i_endidx_in, nproma, i_blk, i_startblk, i_endblk
INTEGER(C_INT) :: i_startidx_out, i_endidx_out
END SUBROUTINE get_indices_v_lib
END INTERFACE
#endif
END MODULE mo_lib_loopindices
src/support/mo_lib_loopindices.cpp 0 → 100644
View file @ e219c7cc
#include <algorithm> // For std::max
extern "C" {
// get_indices_c_lib function
void get_indices_c_lib(int i_startidx_in, int i_endidx_in, int nproma, int i_blk, int i_startblk, int i_endblk,
int &i_startidx_out, int &i_endidx_out) {
if (i_blk == i_startblk) {
i_startidx_out = std::max(1, i_startidx_in);
i_endidx_out = nproma;
if (i_blk == i_endblk) {
i_endidx_out = i_endidx_in;
}
} else if (i_blk == i_endblk) {
i_startidx_out = 1;
i_endidx_out = i_endidx_in;
} else {
i_startidx_out = 1;
i_endidx_out = nproma;
}
}
// get_indices_e_lib function
void get_indices_e_lib(int i_startidx_in, int i_endidx_in, int nproma, int i_blk, int i_startblk, int i_endblk,
int &i_startidx_out, int &i_endidx_out) {
i_startidx_out = (i_blk != i_startblk) ? 1 : std::max(1, i_startidx_in);
i_endidx_out = (i_blk != i_endblk) ? nproma : i_endidx_in;
}
// get_indices_v_lib function
void get_indices_v_lib(int i_startidx_in, int i_endidx_in, int nproma, int i_blk, int i_startblk, int i_endblk,
int &i_startidx_out, int &i_endidx_out) {
if (i_blk == i_startblk) {
i_startidx_out = i_startidx_in;
i_endidx_out = nproma;
if (i_blk == i_endblk) {
i_endidx_out = i_endidx_in;
}
} else if (i_blk == i_endblk) {
i_startidx_out = 1;
i_endidx_out = i_endidx_in;
} else {
i_startidx_out = 1;
i_endidx_out = nproma;
}
}
}
src/support/mo_math_utilities.F90
View file @ e219c7cc
......@@ -22,6 +22,7 @@
! #endif
MODULE mo_math_utilities
USE, INTRINSIC :: ISO_C_BINDING
USE mo_iconlib_kind, ONLY: wp, dp, sp
USE mo_math_constants, ONLY: pi, pi_2, dbl_eps
USE mo_gridman_constants, ONLY: SUCCESS, TORUS_MAX_LAT
......@@ -160,7 +161,98 @@ MODULE mo_math_utilities
CHARACTER(LEN=*), PARAMETER :: modname = 'mo_math_utilities'
!-------------------------------------------------------------------------
!>
!! TDMA tridiagonal matrix solver for a_i*x_(i-1) + b_i*x_i + c_i*x_(i+1) = d_i
!!
!! a - sub-diagonal (means it is the diagonal below the main diagonal)
!! b - the main diagonal
!! c - sup-diagonal (means it is the diagonal above the main diagonal)
!! d - right part
!! varout - the answer (identical to x in description above)
!! slev - start level (top)
!! elev - end level (bottom)
! Preprocessor directive to conditionally include the tdma_solver_vec implementation
#ifndef __USE_CPP_BINDINGS
CONTAINS
SUBROUTINE tdma_solver_vec(a, b, c, d, slev, elev, startidx, endidx, varout, opt_acc_queue)
INTEGER, INTENT(IN) :: slev, elev
INTEGER, INTENT(IN) :: startidx, endidx
REAL(wp), INTENT(IN) :: a(:, :), b(:, :), c(:, :), d(:, :)
REAL(wp), INTENT(OUT) :: varout(:, :)
INTEGER, OPTIONAL, INTENT(IN) :: opt_acc_queue
!
! local
REAL(wp):: m, c_p(SIZE(a, 1), SIZE(a, 2)), d_p(SIZE(a, 1), SIZE(a, 2))
INTEGER :: i
INTEGER :: jc
INTEGER :: acc_queue
IF (PRESENT(opt_acc_queue)) THEN
acc_queue = opt_acc_queue
ELSE
acc_queue = 1
END IF
! initialize c-prime and d-prime
!$ACC PARALLEL DEFAULT(PRESENT) CREATE(c_p, d_p) ASYNC(acc_queue)
!$ACC LOOP GANG(STATIC: 1) VECTOR
DO jc = startidx, endidx
c_p(jc, slev) = c(jc, slev)/b(jc, slev)
d_p(jc, slev) = d(jc, slev)/b(jc, slev)
END DO
! solve for vectors c-prime and d-prime
!$ACC LOOP SEQ
!NEC$ outerloop_unroll(4)
DO i = slev + 1, elev
!$ACC LOOP GANG(STATIC: 1) VECTOR PRIVATE(m)
DO jc = startidx, endidx
m = 1._wp/(b(jc, i) - c_p(jc, i - 1)*a(jc, i))
c_p(jc, i) = c(jc, i)*m
d_p(jc, i) = (d(jc, i) - d_p(jc, i - 1)*a(jc, i))*m
END DO
END DO
! initialize varout
!$ACC LOOP GANG(STATIC: 1) VECTOR
DO jc = startidx, endidx
varout(jc, elev) = d_p(jc, elev)
END DO
! solve for varout from the vectors c-prime and d-prime
!$ACC LOOP SEQ
!NEC$ outerloop_unroll(4)
DO i = elev - 1, slev, -1
!$ACC LOOP GANG(STATIC: 1) VECTOR
DO jc = startidx, endidx
varout(jc, i) = d_p(jc, i) - c_p(jc, i)*varout(jc, i + 1)
END DO
END DO
!$ACC END PARALLEL
IF (.NOT. PRESENT(opt_acc_queue)) THEN
!$ACC WAIT(acc_queue)
END IF
END SUBROUTINE tdma_solver_vec
#else
! C++ binding for tdma_solver_vec
INTERFACE
SUBROUTINE tdma_solver_vec(a, b, c, d, slev, elev, startidx, endidx, varout, opt_acc_queue) BIND(C, NAME="tdma_solver_vec")
IMPORT :: C_DOUBLE, C_INT
REAL(C_DOUBLE), INTENT(IN) :: a(*), b(*), c(*), d(*)
INTEGER(C_INT), VALUE :: slev, elev, startidx, endidx
REAL(C_DOUBLE), INTENT(OUT) :: varout(*)
INTEGER(C_INT), OPTIONAL :: opt_acc_queue
END SUBROUTINE tdma_solver_vec
END INTERFACE
CONTAINS
#endif
!-------------------------------------------------------------------------
! Variant for double-precision (or working-precision=dp) lon+lat in ICON
......@@ -2019,78 +2111,6 @@ CONTAINS
END SUBROUTINE tdma_solver
!-------------------------------------------------------------------------
!>
!! TDMA tridiagonal matrix solver for a_i*x_(i-1) + b_i*x_i + c_i*x_(i+1) = d_i
!!
!! a - sub-diagonal (means it is the diagonal below the main diagonal)
!! b - the main diagonal
!! c - sup-diagonal (means it is the diagonal above the main diagonal)
!! d - right part
!! varout - the answer (identical to x in description above)
!! slev - start level (top)
!! elev - end level (bottom)
SUBROUTINE tdma_solver_vec(a, b, c, d, slev, elev, startidx, endidx, varout, opt_acc_queue)
INTEGER, INTENT(IN) :: slev, elev
INTEGER, INTENT(IN) :: startidx, endidx
REAL(wp), INTENT(IN) :: a(:, :), b(:, :), c(:, :), d(:, :)
REAL(wp), INTENT(OUT) :: varout(:, :)
INTEGER, OPTIONAL, INTENT(IN) :: opt_acc_queue
!
! local
REAL(wp):: m, c_p(SIZE(a, 1), SIZE(a, 2)), d_p(SIZE(a, 1), SIZE(a, 2))
INTEGER :: i
INTEGER :: jc
INTEGER :: acc_queue
IF (PRESENT(opt_acc_queue)) THEN
acc_queue = opt_acc_queue
ELSE
acc_queue = 1
END IF
! initialize c-prime and d-prime
!$ACC PARALLEL DEFAULT(PRESENT) CREATE(c_p, d_p) ASYNC(acc_queue)
!$ACC LOOP GANG(STATIC: 1) VECTOR
DO jc = startidx, endidx
c_p(jc, slev) = c(jc, slev)/b(jc, slev)
d_p(jc, slev) = d(jc, slev)/b(jc, slev)
END DO
! solve for vectors c-prime and d-prime
!$ACC LOOP SEQ
!NEC$ outerloop_unroll(4)
DO i = slev + 1, elev
!$ACC LOOP GANG(STATIC: 1) VECTOR PRIVATE(m)
DO jc = startidx, endidx
m = 1._wp/(b(jc, i) - c_p(jc, i - 1)*a(jc, i))
c_p(jc, i) = c(jc, i)*m
d_p(jc, i) = (d(jc, i) - d_p(jc, i - 1)*a(jc, i))*m
END DO
END DO
! initialize varout
!$ACC LOOP GANG(STATIC: 1) VECTOR
DO jc = startidx, endidx
varout(jc, elev) = d_p(jc, elev)
END DO
! solve for varout from the vectors c-prime and d-prime
!$ACC LOOP SEQ
!NEC$ outerloop_unroll(4)
DO i = elev - 1, slev, -1
!$ACC LOOP GANG(STATIC: 1) VECTOR
DO jc = startidx, endidx
varout(jc, i) = d_p(jc, i) - c_p(jc, i)*varout(jc, i + 1)
END DO
END DO
!$ACC END PARALLEL
IF (.NOT. PRESENT(opt_acc_queue)) THEN
!$ACC WAIT(acc_queue)
END IF
END SUBROUTINE tdma_solver_vec
!-------------------------------------------------------------------------
!-------------------------------------------------------------------------
!
!> Helper functions for computing the vertical layer structure
......
src/support/mo_math_utilities.cpp 0 → 100644
View file @ e219c7cc
#include <vector>
#include <iostream>
#include <chrono> // For timing
extern "C" {
void tdma_solver_vec(double *a, double *b, double *c, double *d,
int slev, int elev, int startidx, int endidx,
double* varout, int opt_acc_queue = -1) {
int acc_queue = (opt_acc_queue == -1) ? 1 : opt_acc_queue; // Use 1 as the default if opt_acc_queue is not provided
// Determine array sizes based on startidx and endidx
int nrows = endidx - startidx;
int ncols = elev - slev;
// Temporary arrays for c-prime and d-prime
std::vector<double> cp(nrows * ncols, 0.0);
std::vector<double> dp(nrows * ncols, 0.0);
// Helper function to access 2D arrays stored as 1D
auto idx = [&](int row, int col) { return col * nrows + row; };
// Start timing
auto start_time = std::chrono::high_resolution_clock::now();
// OpenACC Parallel Region
#pragma acc parallel default(present) create(cp[:nrows*ncols], dp[:nrows*ncols]) async(acc_queue)
{
// Initialize c-prime and d-prime
#pragma acc loop gang(static: 1) vector
for (int jc = startidx; jc < endidx; ++jc) {
cp[idx(jc, slev)] = c[idx(jc, slev)] / b[idx(jc, slev)];
dp[idx(jc, slev)] = d[idx(jc, slev)] / b[idx(jc, slev)];
}
// Solve for vectors c-prime and d-prime
#pragma acc loop seq
for (int i = slev + 1; i < elev; ++i) {
#pragma acc loop gang(static: 1) vector
for (int jc = startidx; jc < endidx; ++jc) {
double m = 1.0 / (b[idx(jc, i)] - cp[idx(jc, i - 1)] * a[idx(jc, i)]);
cp[idx(jc, i)] = c[idx(jc, i)] * m;
dp[idx(jc, i)] = (d[idx(jc, i)] - dp[idx(jc, i - 1)] * a[idx(jc, i)]) * m;
}
}
// Initialize varout
#pragma acc loop gang(static: 1) vector
for (int jc = startidx; jc < endidx; ++jc) {
varout[idx(jc, elev-1)] = dp[idx(jc, elev-1)];
}
// Solve for varout from the vectors c-prime and d-prime
#pragma acc loop seq
for (int i = elev - 2; i >= slev; --i) {
#pragma acc loop gang(static: 1) vector
for (int jc = startidx; jc < endidx; ++jc) {
varout[idx(jc, i)] = dp[idx(jc, i)] - cp[idx(jc, i)] * varout[idx(jc, i + 1)];
}
}
}
printf("tdma_solver_vec: completed using C++\n");
// Wait for OpenACC asynchronous operations to complete if acc_queue is not optional
if (opt_acc_queue == -1) {
#pragma acc wait(acc_queue)
}
// End timing
auto end_time = std::chrono::high_resolution_clock::now();
std::chrono::duration<double> elapsed_time = end_time - start_time;
std::cout << "Elapsed time for tdma_solver_vec (C++): " << elapsed_time.count() << " seconds" << std::endl;
}
}
test/fortran/test_loopindices.f90 0 → 100644
View file @ e219c7cc
! ICON
!
! ---------------------------------------------------------------
! Copyright (C) 2004-2024, DWD, MPI-M, DKRZ, KIT, ETH, MeteoSwiss
! Contact information: icon-model.org
!
! See AUTHORS.TXT for a list of authors
! See LICENSES/ for license information
! SPDX-License-Identifier: BSD-3-Clause
! ---------------------------------------------------------------
MODULE TEST_mo_lib_loopindices
USE FORTUTF
IMPLICIT NONE
PRIVATE
PUBLIC :: TEST_get_indices_lib
CONTAINS
SUBROUTINE TEST_get_indices_lib
USE mo_lib_loopindices, ONLY: get_indices_c_lib, get_indices_e_lib, get_indices_v_lib
INTEGER :: i_startidx_in ! Start index as input
INTEGER :: i_endidx_in ! End index as input
INTEGER :: nproma ! inner loop length/vector length
INTEGER :: i_blk ! Current block (variable jb in do loops)
INTEGER :: i_startblk ! Start block of do loop
INTEGER :: i_endblk ! End block of do loop
INTEGER :: i_startidx_out, i_endidx_out ! Start and end indices (jc loop), as output
i_startidx_in = 2
i_endidx_in = 15
nproma = 32
i_startblk = 1
i_endblk = 40
! CASE: I -> i_blk == i_startblk
i_blk = 1
CALL get_indices_c_lib(i_startidx_in, i_endidx_in, nproma, i_blk, i_startblk, i_endblk, &
i_startidx_out, i_endidx_out)
CALL TAG_TEST("TEST_get_indices_c_start_1")
CALL ASSERT_EQUAL(i_startidx_out, MAX(1, i_startidx_in))
CALL TAG_TEST("TEST_get_indices_c_end_1")
CALL ASSERT_EQUAL(i_endidx_out, nproma)
i_startidx_out = 0
i_endidx_out = 0
CALL get_indices_e_lib(i_startidx_in, i_endidx_in, nproma, i_blk, i_startblk, i_endblk, &
i_startidx_out, i_endidx_out)
CALL TAG_TEST("TEST_get_indices_e_start_1")
CALL ASSERT_EQUAL(i_startidx_out, MAX(1, i_startidx_in))
CALL TAG_TEST("TEST_get_indices_e_end_1")
CALL ASSERT_EQUAL(i_endidx_out, nproma)
i_startidx_out = 0
i_endidx_out = 0
CALL get_indices_v_lib(i_startidx_in, i_endidx_in, nproma, i_blk, i_startblk, i_endblk, &
i_startidx_out, i_endidx_out)
CALL TAG_TEST("TEST_get_indices_v_start_1")
CALL ASSERT_EQUAL(i_startidx_out, MAX(1, i_startidx_in))
CALL TAG_TEST("TEST_get_indices_v_end_1")
CALL ASSERT_EQUAL(i_endidx_out, nproma)
i_startidx_out = 0
i_endidx_out = 0
! CASE: II -> i_blk == i_endblk
i_blk = 40
CALL get_indices_c_lib(i_startidx_in, i_endidx_in, nproma, i_blk, i_startblk, i_endblk, &
i_startidx_out, i_endidx_out)
CALL TAG_TEST("TEST_get_indices_c_start_2")
CALL ASSERT_EQUAL(i_startidx_out, 1)
CALL TAG_TEST("TEST_get_indices_c_end_2")
CALL ASSERT_EQUAL(i_endidx_out, i_endidx_in)
i_startidx_out = 0
i_endidx_out = 0
CALL get_indices_e_lib(i_startidx_in, i_endidx_in, nproma, i_blk, i_startblk, i_endblk, &
i_startidx_out, i_endidx_out)
CALL TAG_TEST("TEST_get_indices_e_start_2")
CALL ASSERT_EQUAL(i_startidx_out, 1)
CALL TAG_TEST("TEST_get_indices_e_end_2")
CALL ASSERT_EQUAL(i_endidx_out, i_endidx_in)
i_startidx_out = 0
i_endidx_out = 0
CALL get_indices_v_lib(i_startidx_in, i_endidx_in, nproma, i_blk, i_startblk, i_endblk, &
i_startidx_out, i_endidx_out)
CALL TAG_TEST("TEST_get_indices_v_start_2")
CALL ASSERT_EQUAL(i_startidx_out, 1)
CALL TAG_TEST("TEST_get_indices_v_end_2")
CALL ASSERT_EQUAL(i_endidx_out, i_endidx_in)
i_startidx_out = 0
i_endidx_out = 0
! CASE: III -> Every other cases
i_blk = 20
CALL get_indices_c_lib(i_startidx_in, i_endidx_in, nproma, i_blk, i_startblk, i_endblk, &
i_startidx_out, i_endidx_out)
CALL TAG_TEST("TEST_get_indices_c_start_3")
CALL ASSERT_EQUAL(i_startidx_out, 1)
CALL TAG_TEST("TEST_get_indices_c_end_3")
CALL ASSERT_EQUAL(i_endidx_out, nproma)
i_startidx_out = 0
i_endidx_out = 0
CALL get_indices_e_lib(i_startidx_in, i_endidx_in, nproma, i_blk, i_startblk, i_endblk, &
i_startidx_out, i_endidx_out)
CALL TAG_TEST("TEST_get_indices_e_start_3")
CALL ASSERT_EQUAL(i_startidx_out, 1)
CALL TAG_TEST("TEST_get_indices_e_end_3")
CALL ASSERT_EQUAL(i_endidx_out, nproma)
i_startidx_out = 0
i_endidx_out = 0
CALL get_indices_v_lib(i_startidx_in, i_endidx_in, nproma, i_blk, i_startblk, i_endblk, &
i_startidx_out, i_endidx_out)
CALL TAG_TEST("TEST_get_indices_v_start_3")
CALL ASSERT_EQUAL(i_startidx_out, 1)
CALL TAG_TEST("TEST_get_indices_v_end_3")
CALL ASSERT_EQUAL(i_endidx_out, nproma)
END SUBROUTINE TEST_get_indices_lib
END MODULE TEST_mo_lib_loopindices
test/fortran/test_math_utilities.f90 0 → 100644
View file @ e219c7cc
! ICON
!
! ---------------------------------------------------------------
! Copyright (C) 2004-2024, DWD, MPI-M, DKRZ, KIT, ETH, MeteoSwiss
! Contact information: icon-model.org
!
! See AUTHORS.TXT for a list of authors
! See LICENSES/ for license information
! SPDX-License-Identifier: BSD-3-Clause
! ---------------------------------------------------------------
MODULE TEST_mo_math_utilities
USE FORTUTF
USE mo_math_types, ONLY: t_cartesian_coordinates, t_geographical_coordinates, &
& t_line, t_tangent_vectors
USE mo_math_constants, ONLY: pi, pi_2, pi_4
USE mo_lib_grid_geometry_info
! USE mo_physical_constants, ONLY: earth_radius
USE, INTRINSIC :: ISO_FORTRAN_ENV, ONLY: wp => real64
IMPLICIT NONE
PRIVATE
PUBLIC :: TEST_cc2gc, TEST_gc2cc, TEST_cc2tv, TEST_gvec2cvec, TEST_cvec2gvec, TEST_tdma_solver_vec
REAL(wp), PARAMETER :: earth_radius = 6.371229e6_wp !! [m] average radius
CONTAINS
SUBROUTINE TEST_cc2gc
USE mo_math_utilities, ONLY: cc2gc
TYPE(t_cartesian_coordinates) :: coord
TYPE(t_grid_geometry_info) :: geometry_info
TYPE(t_geographical_coordinates) :: pos
REAL(wp) :: lon_ref, lat_ref
coord%x(1) = 10.0_wp
coord%x(2) = 10.0_wp
coord%x(3) = 5.0_wp
geometry_info%geometry_type = sphere_geometry
pos = cc2gc(coord, geometry_info)
lon_ref = 0.78539816339744828_wp
lat_ref = 0.33983690945412193_wp
CALL TAG_TEST("TEST_cc2gc_sphere_lon")
CALL ASSERT_EQUAL(pos%lon, lon_ref)
CALL TAG_TEST("TEST_cc2gc_sphere_lat")
CALL ASSERT_EQUAL(pos%lat, lat_ref)
geometry_info%geometry_type = planar_torus_geometry
geometry_info%domain_length = 2.0_wp*pi*earth_radius
geometry_info%domain_height = 2.0_wp*pi*earth_radius
coord%x(1) = 1000000.0_wp
coord%x(2) = 100.0_wp
pos = cc2gc(coord, geometry_info)
lon_ref = 0.15695558894524117_wp
lat_ref = -0.17453205322393880_wp
CALL TAG_TEST("TEST_cc2gc_torus_lon")
CALL ASSERT_EQUAL(pos%lon, lon_ref)
CALL TAG_TEST("TEST_cc2gc_torus_lat")
CALL ASSERT_EQUAL(pos%lat, lat_ref)
END SUBROUTINE TEST_cc2gc
SUBROUTINE TEST_gc2cc
USE mo_math_utilities, ONLY: gc2cc
TYPE(t_cartesian_coordinates) :: coord, coord_ref
TYPE(t_grid_geometry_info) :: geometry_info
TYPE(t_geographical_coordinates) :: pos
REAL(wp) :: lon_ref, lat_ref, tol
INTEGER :: i
CHARACTER(LEN=32) :: tag
tol = 1d-15
pos%lon = 0.78_wp
pos%lat = 0.34_wp
geometry_info%geometry_type = sphere_geometry
coord = gc2cc(pos, geometry_info)
coord_ref%x(1) = 0.6702170547483377_wp
coord_ref%x(2) = 0.6630199536212522_wp
coord_ref%x(3) = 0.3334870921408144_wp
DO i = 1, SIZE(coord%x)
! write(*,"(i4,a,f24.16)") i, ' coord%x(i): ', coord%x(i)
WRITE (tag, '(a,i1)') "TEST_gc2cc_sphere_", i
CALL TAG_TEST(tag)
CALL ASSERT_ALMOST_EQUAL(coord%x(i), coord_ref%x(i), tol)
END DO
geometry_info%geometry_type = planar_torus_geometry
geometry_info%domain_length = 2.0_wp*pi*earth_radius
geometry_info%domain_height = 2.0_wp*pi*earth_radius
pos%lon = 0.15_wp
pos%lat = -0.17_wp
coord = gc2cc(pos, geometry_info)
coord_ref%x(1) = 955684.3499999998603016_wp
coord_ref%x(2) = 519845.4807382422150113_wp
coord_ref%x(3) = 0.0_wp
DO i = 1, SIZE(coord%x)
! write(*,"(i4,a,f24.16)") i, ' coord%x(i): ', coord%x(i)
WRITE (tag, '(a,i1)') "TEST_gc2cc_torus_", i
CALL TAG_TEST(tag)
CALL ASSERT_EQUAL(coord%x(i), coord_ref%x(i))
END DO
END SUBROUTINE TEST_gc2cc
SUBROUTINE TEST_cc2tv
USE mo_math_utilities, ONLY: cc2tv
TYPE(t_cartesian_coordinates) :: coord
TYPE(t_geographical_coordinates) :: pos
TYPE(t_tangent_vectors):: tt
REAL(wp) :: v1_ref, v2_ref
pos%lon = pi_2
pos%lat = pi_4
coord%x(1) = 10.0_wp
coord%x(2) = 15.0_wp
coord%x(3) = 5.0_wp
tt = cc2tv(coord, pos)
v1_ref = -9.9999999999999982_wp
v2_ref = -7.0710678118654737_wp
CALL TAG_TEST("TEST_cc2tv_v1")
CALL ASSERT_EQUAL(tt%v1, v1_ref)
CALL TAG_TEST("TEST_cc2tv_v2")
CALL ASSERT_EQUAL(tt%v2, v2_ref)
END SUBROUTINE TEST_cc2tv
SUBROUTINE TEST_gvec2cvec
USE mo_math_utilities, ONLY: gvec2cvec
REAL(wp) :: p_gu, p_gv ! zonal and meridional vec. component
REAL(wp) :: p_long, p_lat ! geo. coord. of data point
TYPE(t_grid_geometry_info) :: geometry_info
REAL(wp) :: p_cu, p_cv, p_cw ! Cart. vector
REAL(wp) :: p_cu_ref, p_cv_ref, p_cw_ref ! Cart. vector ref
geometry_info%geometry_type = sphere_geometry
p_gu = 10.0_wp
p_gv = 5.0_wp
p_long = pi_2
p_lat = pi_4
CALL gvec2cvec(p_gu, p_gv, p_long, p_lat, p_cu, p_cv, p_cw, geometry_info)
p_cu_ref = -10.0_wp
p_cv_ref = -3.5355339059327369_wp
p_cw_ref = 3.5355339059327378_wp
CALL TAG_TEST("TEST_gvec2cvec_sphere_cu")
CALL ASSERT_EQUAL(p_cu, p_cu_ref)
CALL TAG_TEST("TEST_gvec2cvec_sphere_cv")
CALL ASSERT_EQUAL(p_cv, p_cv_ref)
CALL TAG_TEST("TEST_gvec2cvec_sphere_cw")
CALL ASSERT_EQUAL(p_cw, p_cw_ref)
geometry_info%geometry_type = planar_torus_geometry
CALL gvec2cvec(p_gu, p_gv, p_long, p_lat, p_cu, p_cv, p_cw, geometry_info)
CALL TAG_TEST("TEST_gvec2cvec_torus_cu")
CALL ASSERT_EQUAL(p_cu, p_gu)
CALL TAG_TEST("TEST_gvec2cvec_torus_cv")
CALL ASSERT_EQUAL(p_cv, p_gv)
CALL TAG_TEST("TEST_gvec2cvec_torus_cw")
CALL ASSERT_EQUAL(p_cw, 0.0_wp)
END SUBROUTINE TEST_gvec2cvec
SUBROUTINE TEST_cvec2gvec
USE mo_math_utilities, ONLY: cvec2gvec
REAL(wp) :: p_cu, p_cv, p_cw ! Cart. vector
REAL(wp) :: p_long, p_lat ! geo. coord. of data point
TYPE(t_grid_geometry_info) :: geometry_info
REAL(wp) :: p_gu, p_gv ! zonal and meridional vec. comp.
REAL(wp) :: p_gu_ref, p_gv_ref
geometry_info%geometry_type = sphere_geometry
p_cu = -10.0_wp
p_cv = -3.0_wp
p_cw = 3.0_wp
p_long = pi_2
p_lat = pi_4
CALL cvec2gvec(p_cu, p_cv, p_cw, p_long, p_lat, p_gu, p_gv, geometry_info)
! write(*,"(a,f24.16)") ' p_gu: ', p_gu
! write(*,"(a,f24.16)") ' p_gv: ', p_gv
p_gu_ref = 10.0_wp
p_gv_ref = 4.2426406871192857_wp
CALL TAG_TEST("TEST_cvec2gvec_sphere_gu")
CALL ASSERT_EQUAL(p_gu, p_gu_ref)
CALL TAG_TEST("TEST_cvec2gvec_sphere_gv")
CALL ASSERT_EQUAL(p_gv, p_gv_ref)
geometry_info%geometry_type = planar_torus_geometry
CALL cvec2gvec(p_cu, p_cv, p_cw, p_long, p_lat, p_gu, p_gv, geometry_info)
CALL TAG_TEST("TEST_cvec2gvec_torus_gu")
CALL ASSERT_EQUAL(p_gu, p_cu)
CALL TAG_TEST("TEST_cvec2gvec_torus_gv")
CALL ASSERT_EQUAL(p_gv, p_cv)
END SUBROUTINE TEST_cvec2gvec
SUBROUTINE TEST_tdma_solver_vec
USE mo_math_utilities, ONLY: tdma_solver_vec
INTEGER, PARAMETER :: n = 10
REAL(wp) :: a(n, n), b(n, n), c(n, n), d(n, n), x(n, n)
INTEGER :: i, j
REAL(wp) :: sum, sum_ref
DO i = 1, n
DO j = 1, n
a(i, j) = 1.0_wp
b(i, j) = 2.0_wp
c(i, j) = 1.0_wp
d(i, j) = 1.0_wp
END DO
END DO
CALL tdma_solver_vec(a, b, c, d, 1, n, 1, n, x)
sum = 0.0_wp
DO i = 1, n
sum = sum + x(i, 1)
END DO
sum_ref = 4.5454545454545467_wp
CALL TAG_TEST("TEST_tdma_solver_vec")
CALL ASSERT_EQUAL(sum, sum_ref)
END SUBROUTINE TEST_tdma_solver_vec
END MODULE TEST_mo_math_utilities

Imprint | Privacy Policy