From 9199211ab072e41d78b7c9d1f2e89cf39c8f4c81 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Einar=20=C3=96rn=20=C3=93lason?= Date: Mon, 3 Jun 2024 12:54:32 +0200 Subject: [PATCH 01/33] Polynya test case Initial commit of a python script to create the setup and a config file to run the test. TODO: - Change the wind direction - Test BBM - Hourly output doesn't work --- run/config_polynya.cfg | 32 ++++++++++++++++++++++++++++++ run/make_init_polynya.py | 43 ++++++++++++++++++++++++++++++++++++++++ 2 files changed, 75 insertions(+) create mode 100644 run/config_polynya.cfg create mode 100644 run/make_init_polynya.py diff --git a/run/config_polynya.cfg b/run/config_polynya.cfg new file mode 100644 index 000000000..85e4af0be --- /dev/null +++ b/run/config_polynya.cfg @@ -0,0 +1,32 @@ +[model] +init_file = init_polynya.nc +start = 2010-01-01T00:00:00Z +stop = 2010-01-02T00:00:00Z +time_step = P0-0T00:01:00 + +[Modules] +DiagnosticOutputModule = Nextsim::ConfigOutput +DynamicsModule = Nextsim::MEVPDynamics +IceAlbedoModule = Nextsim::WintonAlbedo +AtmosphereBoundaryModule = Nextsim::FluxConfiguredAtmosphere +OceanBoundaryModule = Nextsim::FluxConfiguredOcean +IceThermodynamicsModule = Nextsim::ThermoWinton + +[ConfigOutput] +field_names = hsnow,hice,tice +period = P0-0T01:00:00 + +[FluxConfiguredOcean] +qio = 2.30 +sss = 28 +sst = -1.54 +mld = 10. +current_u = 0. +current_v = 0. + +[FluxConfiguredAtmosphere] +Q_ia = 2 +Q_ow = 200 +dQia_dT = 0 +wind_u = 20 +wind_v = 0 diff --git a/run/make_init_polynya.py b/run/make_init_polynya.py new file mode 100644 index 000000000..8317ac353 --- /dev/null +++ b/run/make_init_polynya.py @@ -0,0 +1,43 @@ +from make_init_base import initMaker + +# Creates initial conditions for the Bjornsson et al. (2001) polynya case + +# Domain size [km] +x = 100 +y = 50 +res = 2 + +nfirst = int(y / res) +nsecond = int(x / res) +nLayers = 3 + +fname = f"init_polynya.nc" + +# The model expects everything in metres +initializer = initMaker(fname, nfirst, nsecond, nLayers, res*1e3) + +# Ice everywhere and all boundaries closed, except the x = 100 km end +initializer.mask[:, :] = 1. +initializer.mask[0, :] = 0. +initializer.mask[-1, :] = 0. +initializer.mask[:, 0] = 0. + +# Uniform concentration of 90% +initializer.cice[:, :] = 0.9 + +# Uniform thickness of 20 cm +initializer.hice[:, :] = 0.2 + +# Ice and ocean temperature and salinity at the freezing point +ice_salinity = 5 # should match Ice::s in constants.hpp +mu: float = -0.055 # should match Water::mu in constants.hpp +ocean_temperature = -1.54 +ocean_salinity = ocean_temperature / mu + +initializer.sss[:, :] = ocean_salinity +initializer.sst[:, :] = ocean_temperature +initializer.tice[:, :, :] = ice_salinity * mu + +# All other variables are zero or not needed + +# The file is written when initializer goes out of scope From db69db8243cb3d4f90543bfa1cbc3b4845ce739f Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Einar=20=C3=96rn=20=C3=93lason?= Date: Fri, 14 Jun 2024 10:28:39 +0200 Subject: [PATCH 02/33] Minor updates to the config file Most importantly, the wind is blowing at an angle now. --- run/config_polynya.cfg | 17 +++++++++-------- 1 file changed, 9 insertions(+), 8 deletions(-) diff --git a/run/config_polynya.cfg b/run/config_polynya.cfg index 85e4af0be..a65c97292 100644 --- a/run/config_polynya.cfg +++ b/run/config_polynya.cfg @@ -1,20 +1,21 @@ [model] init_file = init_polynya.nc -start = 2010-01-01T00:00:00Z -stop = 2010-01-02T00:00:00Z -time_step = P0-0T00:01:00 +start = 2023-01-01T00:00:00Z +stop = 2023-01-03T00:00:00Z +time_step = P0-0T00:02:00 +missing_value = 1e20 [Modules] DiagnosticOutputModule = Nextsim::ConfigOutput DynamicsModule = Nextsim::MEVPDynamics -IceAlbedoModule = Nextsim::WintonAlbedo AtmosphereBoundaryModule = Nextsim::FluxConfiguredAtmosphere OceanBoundaryModule = Nextsim::FluxConfiguredOcean IceThermodynamicsModule = Nextsim::ThermoWinton [ConfigOutput] -field_names = hsnow,hice,tice -period = P0-0T01:00:00 +period = P0-0T04:00:00 +field_names = hice,cice,u,v +filename = polynya.diagnostic.nc [FluxConfiguredOcean] qio = 2.30 @@ -28,5 +29,5 @@ current_v = 0. Q_ia = 2 Q_ow = 200 dQia_dT = 0 -wind_u = 20 -wind_v = 0 +wind_u = 16 +wind_v = 12 From 805e1c0ca6c5f53653d043089aaeeafa9061c45a Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Einar=20=C3=96rn=20=C3=93lason?= Date: Fri, 14 Jun 2024 10:33:00 +0200 Subject: [PATCH 03/33] Create init file with CMake --- run/CMakeLists.txt | 1 + 1 file changed, 1 insertion(+) diff --git a/run/CMakeLists.txt b/run/CMakeLists.txt index b84cecd2b..ad60d7d49 100644 --- a/run/CMakeLists.txt +++ b/run/CMakeLists.txt @@ -4,3 +4,4 @@ execute_process(COMMAND "${Python_EXECUTABLE}" "${CMAKE_CURRENT_SOURCE_DIR}/make execute_process(COMMAND "${Python_EXECUTABLE}" "${CMAKE_CURRENT_SOURCE_DIR}/make_init_column.py" WORKING_DIRECTORY "${CMAKE_CURRENT_SOURCE_DIR}") execute_process(COMMAND "${Python_EXECUTABLE}" "${CMAKE_CURRENT_SOURCE_DIR}/make_init_para24x30.py" WORKING_DIRECTORY "${CMAKE_CURRENT_SOURCE_DIR}") execute_process(COMMAND "${Python_EXECUTABLE}" "${CMAKE_CURRENT_SOURCE_DIR}/make_init_benchmark.py" WORKING_DIRECTORY "${CMAKE_CURRENT_SOURCE_DIR}") +execute_process(COMMAND "${Python_EXECUTABLE}" "${CMAKE_CURRENT_SOURCE_DIR}/make_init_polynya.py" WORKING_DIRECTORY "${CMAKE_CURRENT_SOURCE_DIR}") From b99de51a24ca57c1d27f4328f34d6c0aeb7f1c7c Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Einar=20=C3=96rn=20=C3=93lason?= Date: Fri, 21 Feb 2025 15:08:59 +0100 Subject: [PATCH 04/33] Add open (Neumann zero) boundary conditions Adapts the dirichletFromEdge function to be neumannFromEdge and adds two neumannZero functions - one for DGVectors and one for CGVectors. Both are inspired by the existing dirichletZero function. The neumannZero function must be called after advection - which is a bit worse, but moving it into the DGTransport class is compicated. --- dynamics/src/CGDynamicsKernel.cpp | 53 ++++++++++++++++++- dynamics/src/ParametricMesh.cpp | 14 ++--- .../src/include/BrittleCGDynamicsKernel.hpp | 6 ++- dynamics/src/include/CGDynamicsKernel.hpp | 3 +- dynamics/src/include/DynamicsKernel.hpp | 44 +++++++++++++-- dynamics/src/include/ParametricMesh.hpp | 9 ++-- dynamics/test/ParametricMesh_test.cpp | 4 +- 7 files changed, 113 insertions(+), 20 deletions(-) diff --git a/dynamics/src/CGDynamicsKernel.cpp b/dynamics/src/CGDynamicsKernel.cpp index 2257ce037..d6427b698 100644 --- a/dynamics/src/CGDynamicsKernel.cpp +++ b/dynamics/src/CGDynamicsKernel.cpp @@ -1,7 +1,7 @@ /*! * @file CGDynamicsKernel.cpp * - * @date 14 Jan 2025 + * @date 21 Feb 2025 * @author Tim Spain */ @@ -418,10 +418,61 @@ void CGDynamicsKernel::dirichletZero(CGVector& v) const } } +template +void CGDynamicsKernel::neumannZero(CGVector& v) const +{ + // the four segments bottom, right, top, left, are each processed in parallel + for (size_t seg = 0; seg < 4; ++seg) { +#pragma omp parallel for + for (size_t i = 0; i < smesh->neumann[seg].size(); ++i) { + + const size_t eid = smesh->neumann[seg][i]; + const size_t ix = eid % smesh->nx; // compute coordinates of element + const size_t iy = eid / smesh->nx; + + switch (seg) { + case 0: // bottom <= top + for (size_t j = 0; j < CGdegree + 1; ++j) + v(iy * CGdegree * (CGdegree * smesh->nx + 1) + CGdegree * ix + j, 0) = v( + (iy + 1) * CGdegree * (CGdegree * smesh->nx + 1) + CGdegree * ix + j, 0); + break; + case 1: // right <= left + for (size_t j = 0; j < CGdegree + 1; ++j) + v(iy * CGdegree * (CGdegree * smesh->nx + 1) + CGdegree * ix + CGdegree + + (CGdegree * smesh->nx + 1) * j, + 0) + = v(iy * CGdegree * (CGdegree * smesh->nx + 1) + CGdegree * ix + + (CGdegree * smesh->nx + 1) * j, + 0); + break; + case 2: // top <= bottom + for (size_t j = 0; j < CGdegree + 1; ++j) + v((iy + 1) * CGdegree * (CGdegree * smesh->nx + 1) + CGdegree * ix + j, 0) + = v(iy * CGdegree * (CGdegree * smesh->nx + 1) + CGdegree * ix + j, 0); + break; + case 3: // left <= right + for (size_t j = 0; j < CGdegree + 1; ++j) + v(iy * CGdegree * (CGdegree * smesh->nx + 1) + CGdegree * ix + + (CGdegree * smesh->nx + 1) * j, + 0) + = v(iy * CGdegree * (CGdegree * smesh->nx + 1) + CGdegree * ix + CGdegree + + (CGdegree * smesh->nx + 1) * j, + 0); + break; + default: + std::cerr << "That should not have happened!" << std::endl; + abort(); + } + } + } +} + template void CGDynamicsKernel::applyBoundaries() { dirichletZero(u); dirichletZero(v); + neumannZero(u); + neumannZero(v); // TODO Periodic boundary conditions. } diff --git a/dynamics/src/ParametricMesh.cpp b/dynamics/src/ParametricMesh.cpp index 1bf15bd2d..3387be641 100644 --- a/dynamics/src/ParametricMesh.cpp +++ b/dynamics/src/ParametricMesh.cpp @@ -1,6 +1,6 @@ /*! * @file ParametricMesh.hpp - * @date 14 Jan 2025 + * @date 21 Feb 2025 * @author Thomas Richter */ @@ -8,6 +8,7 @@ #include "ParametricTools.hpp" +#include #include #include @@ -259,11 +260,11 @@ void ParametricMesh::dirichletFromMask() } /*! - * Add to the dirichlet arrays due to the domain edges according to an edge index. + * Add to the neumann arrays due to the domain edges according to an edge index. * - * @param edge index of the edge to add closed boundary conditions to. + * @param edge index of the edge to add open boundary conditions to. */ -void ParametricMesh::dirichletFromEdge(Edge edge) +void ParametricMesh::neumannFromEdge(Edge edge) { // BOTTOM, RIGHT, TOP, LEFT const std::array start = { 0, nx - 1, nelements - nx, 0 }; @@ -272,10 +273,11 @@ void ParametricMesh::dirichletFromEdge(Edge edge) for (size_t idx = start[edge]; idx < stop[edge]; idx += stride[edge]) { if (landmask[idx]) { - dirichlet[edge].push_back(idx); + neumann[edge].push_back(idx); } } - sortDirichlet(edge); + for (Edge edge : edges) + std::sort(neumann[edge].begin(), neumann[edge].end()); } /*! diff --git a/dynamics/src/include/BrittleCGDynamicsKernel.hpp b/dynamics/src/include/BrittleCGDynamicsKernel.hpp index c1ba25b71..fb20e0b94 100644 --- a/dynamics/src/include/BrittleCGDynamicsKernel.hpp +++ b/dynamics/src/include/BrittleCGDynamicsKernel.hpp @@ -1,7 +1,7 @@ /*! * @file BrittleCGDynamicsKernel.hpp * - * @date 06 Dec 2024 + * @date 21 Feb 2025 * @author Tim Spain * @author Einar Ólason */ @@ -12,7 +12,6 @@ #include "CGDynamicsKernel.hpp" #include "BBMParameters.hpp" -#include "ParametricMap.hpp" #include "StressUpdateStep.hpp" #include "include/constants.hpp" #include @@ -37,6 +36,7 @@ template class BrittleCGDynamicsKernel : public CGDynamicsKern using DynamicsKernel::applyBoundaries; using DynamicsKernel::advectionAndLimits; using DynamicsKernel::dgtransport; + using DynamicsKernel::neumannZero; using CGDynamicsKernel::u; using CGDynamicsKernel::v; @@ -109,6 +109,8 @@ template class BrittleCGDynamicsKernel : public CGDynamicsKern Nextsim::LimitMax(damage, 1.0); Nextsim::LimitMin(damage, 1e-12); + neumannZero(damage); + prepareIteration({ { hiceName, hice }, { ciceName, cice } }); // The timestep for the brittle solver is the solver subtimestep diff --git a/dynamics/src/include/CGDynamicsKernel.hpp b/dynamics/src/include/CGDynamicsKernel.hpp index 554f0c713..c5f902cec 100644 --- a/dynamics/src/include/CGDynamicsKernel.hpp +++ b/dynamics/src/include/CGDynamicsKernel.hpp @@ -1,7 +1,7 @@ /*! * @file CGDynamicsKernel.hpp * - * @date 06 Dec 2024 + * @date 21 Feb 2025 * @author Tim Spain */ @@ -69,6 +69,7 @@ class CGDynamicsKernel : public DynamicsKernel { protected: void addStressTensorCell(const size_t eid, const size_t cx, const size_t cy); void dirichletZero(CGVector&) const; + void neumannZero(CGVector&) const; // CG ice velocity CGVector u; CGVector v; diff --git a/dynamics/src/include/DynamicsKernel.hpp b/dynamics/src/include/DynamicsKernel.hpp index 08f9ecc19..446f7c80a 100644 --- a/dynamics/src/include/DynamicsKernel.hpp +++ b/dynamics/src/include/DynamicsKernel.hpp @@ -1,7 +1,7 @@ /*! * @file DynamicsKernel.hpp * - * @date Jan 5, 2024 + * @date 21 Feb 2025 * @author Tim Spain */ @@ -53,9 +53,9 @@ template class DynamicsKernel { smesh->RotatePoleToGreenland(); smesh->landmaskFromModelArray(mask); smesh->dirichletFromMask(); - // TODO: handle periodic and open edges - for (ParametricMesh::Edge edge : ParametricMesh::edges) { - smesh->dirichletFromEdge(edge); + // TODO: handle periodic edges + for (const ParametricMesh::Edge edge : ParametricMesh::edges) { + smesh->neumannFromEdge(edge); } //! Initialize transport @@ -180,6 +180,9 @@ template class DynamicsKernel { Nextsim::LimitMax(cice, 1.0); Nextsim::LimitMin(cice, 0.0); Nextsim::LimitMin(hice, 0.0); + + neumannZero(cice); + neumannZero(hice); } protected: @@ -232,6 +235,39 @@ template class DynamicsKernel { */ virtual void prepareAdvection() = 0; + template void neumannZero(DGVector& v) const + { + // the four segments bottom, right, top, left, are each processed in parallel + for (size_t seg = 0; seg < 4; ++seg) { +#pragma omp parallel for + for (size_t i = 0; i < smesh->neumann[seg].size(); ++i) { + + const size_t eid = smesh->neumann[seg][i]; + const size_t ix = eid % smesh->nx; // compute coordinates of element + const size_t iy = eid / smesh->nx; + + switch (seg) { + case 0: // bottom <= top + for (size_t j = 0; j < DG; ++j) + v(iy * smesh->nx + ix, j) = v((iy + 1) * smesh->nx + ix, j); + break; + case 1: // right <= left + for (size_t j = 0; j < DG; ++j) + v(iy * smesh->nx + ix, j) = v(iy * smesh->nx + ix - 1, j); + break; + case 2: // top <= bottom + for (size_t j = 0; j < DG; ++j) + v(iy * smesh->nx + ix, j) = v((iy - 1) * smesh->nx + ix, j); + break; + case 3: // left <= right + for (size_t j = 0; j < DG; ++j) + v(iy * smesh->nx + ix, j) = v(iy * smesh->nx + ix + 1, j); + break; + } + } + } + } + private: std::unordered_map> advectedFields; diff --git a/dynamics/src/include/ParametricMesh.hpp b/dynamics/src/include/ParametricMesh.hpp index 25e566cf1..c34c8f559 100644 --- a/dynamics/src/include/ParametricMesh.hpp +++ b/dynamics/src/include/ParametricMesh.hpp @@ -1,6 +1,6 @@ /*! * @file ParametricMesh.hpp - * @date 09 Jan 2025 + * @date 21 Feb 2025 * @author Thomas Richter */ @@ -56,9 +56,9 @@ class ParametricMesh { Eigen::Matrix vertices; // stores the /*! - * Stores the Dirichlet boundary information + * Stores the Dirichlet and Neumann boundary information * - * Dirichlet-Information is stored in 4 vectors + * Dirichlet- and Neumann-Information is stored in 4 vectors, e.g. * dirichlet[0]: List of elements that have a Dirichlet-Boundary on the bottom * dirichlet[1]: List of elements that have a Dirichlet-Boundary on the right * dirichlet[2]: List of elements that have a Dirichlet-Boundary on the top @@ -67,6 +67,7 @@ class ParametricMesh { * Each of the vectors is processed in parallel */ std::array, 4> dirichlet; + std::array, 4> neumann; /*! * Periodic: @@ -403,7 +404,7 @@ class ParametricMesh { * * @param edge index of the edge to add closed boundary conditions to. */ - void dirichletFromEdge(Edge edge); + void neumannFromEdge(Edge edge); /*! * Sort all the dirichlet arrays, so the element indices are ordered. diff --git a/dynamics/test/ParametricMesh_test.cpp b/dynamics/test/ParametricMesh_test.cpp index 28d1d243d..a16321636 100644 --- a/dynamics/test/ParametricMesh_test.cpp +++ b/dynamics/test/ParametricMesh_test.cpp @@ -3,7 +3,7 @@ * * @brief Test the ParametricMesh class, especially processing from ModelArray files. * - * @date Dec 15, 2023 + * @date 21 Feb 2025 * @author Tim Spain */ @@ -122,7 +122,7 @@ TEST_CASE("Compare readmesh and landmask reading") // Dirichlet conditions fromArrays.dirichletFromMask(); for (ParametricMesh::Edge edge : ParametricMesh::edges) { - fromArrays.dirichletFromEdge(edge); + fromArrays.neumannFromEdge(edge); } fromArrays.sortDirichlet(); for (ParametricMesh::Edge edge : ParametricMesh::edges) { From 60f876bd1d18c91cda483cb60dcf0392fef637bb Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Einar=20=C3=96rn=20=C3=93lason?= Date: Fri, 14 Mar 2025 09:46:15 +0100 Subject: [PATCH 05/33] Revert "Add open (Neumann zero) boundary conditions" This reverts commit b99de51a24ca57c1d27f4328f34d6c0aeb7f1c7c. --- dynamics/src/CGDynamicsKernel.cpp | 53 +------------------ dynamics/src/ParametricMesh.cpp | 14 +++-- .../src/include/BrittleCGDynamicsKernel.hpp | 6 +-- dynamics/src/include/CGDynamicsKernel.hpp | 3 +- dynamics/src/include/DynamicsKernel.hpp | 44 ++------------- dynamics/src/include/ParametricMesh.hpp | 9 ++-- dynamics/test/ParametricMesh_test.cpp | 4 +- 7 files changed, 20 insertions(+), 113 deletions(-) diff --git a/dynamics/src/CGDynamicsKernel.cpp b/dynamics/src/CGDynamicsKernel.cpp index d6427b698..2257ce037 100644 --- a/dynamics/src/CGDynamicsKernel.cpp +++ b/dynamics/src/CGDynamicsKernel.cpp @@ -1,7 +1,7 @@ /*! * @file CGDynamicsKernel.cpp * - * @date 21 Feb 2025 + * @date 14 Jan 2025 * @author Tim Spain */ @@ -418,61 +418,10 @@ void CGDynamicsKernel::dirichletZero(CGVector& v) const } } -template -void CGDynamicsKernel::neumannZero(CGVector& v) const -{ - // the four segments bottom, right, top, left, are each processed in parallel - for (size_t seg = 0; seg < 4; ++seg) { -#pragma omp parallel for - for (size_t i = 0; i < smesh->neumann[seg].size(); ++i) { - - const size_t eid = smesh->neumann[seg][i]; - const size_t ix = eid % smesh->nx; // compute coordinates of element - const size_t iy = eid / smesh->nx; - - switch (seg) { - case 0: // bottom <= top - for (size_t j = 0; j < CGdegree + 1; ++j) - v(iy * CGdegree * (CGdegree * smesh->nx + 1) + CGdegree * ix + j, 0) = v( - (iy + 1) * CGdegree * (CGdegree * smesh->nx + 1) + CGdegree * ix + j, 0); - break; - case 1: // right <= left - for (size_t j = 0; j < CGdegree + 1; ++j) - v(iy * CGdegree * (CGdegree * smesh->nx + 1) + CGdegree * ix + CGdegree - + (CGdegree * smesh->nx + 1) * j, - 0) - = v(iy * CGdegree * (CGdegree * smesh->nx + 1) + CGdegree * ix - + (CGdegree * smesh->nx + 1) * j, - 0); - break; - case 2: // top <= bottom - for (size_t j = 0; j < CGdegree + 1; ++j) - v((iy + 1) * CGdegree * (CGdegree * smesh->nx + 1) + CGdegree * ix + j, 0) - = v(iy * CGdegree * (CGdegree * smesh->nx + 1) + CGdegree * ix + j, 0); - break; - case 3: // left <= right - for (size_t j = 0; j < CGdegree + 1; ++j) - v(iy * CGdegree * (CGdegree * smesh->nx + 1) + CGdegree * ix - + (CGdegree * smesh->nx + 1) * j, - 0) - = v(iy * CGdegree * (CGdegree * smesh->nx + 1) + CGdegree * ix + CGdegree - + (CGdegree * smesh->nx + 1) * j, - 0); - break; - default: - std::cerr << "That should not have happened!" << std::endl; - abort(); - } - } - } -} - template void CGDynamicsKernel::applyBoundaries() { dirichletZero(u); dirichletZero(v); - neumannZero(u); - neumannZero(v); // TODO Periodic boundary conditions. } diff --git a/dynamics/src/ParametricMesh.cpp b/dynamics/src/ParametricMesh.cpp index 3387be641..1bf15bd2d 100644 --- a/dynamics/src/ParametricMesh.cpp +++ b/dynamics/src/ParametricMesh.cpp @@ -1,6 +1,6 @@ /*! * @file ParametricMesh.hpp - * @date 21 Feb 2025 + * @date 14 Jan 2025 * @author Thomas Richter */ @@ -8,7 +8,6 @@ #include "ParametricTools.hpp" -#include #include #include @@ -260,11 +259,11 @@ void ParametricMesh::dirichletFromMask() } /*! - * Add to the neumann arrays due to the domain edges according to an edge index. + * Add to the dirichlet arrays due to the domain edges according to an edge index. * - * @param edge index of the edge to add open boundary conditions to. + * @param edge index of the edge to add closed boundary conditions to. */ -void ParametricMesh::neumannFromEdge(Edge edge) +void ParametricMesh::dirichletFromEdge(Edge edge) { // BOTTOM, RIGHT, TOP, LEFT const std::array start = { 0, nx - 1, nelements - nx, 0 }; @@ -273,11 +272,10 @@ void ParametricMesh::neumannFromEdge(Edge edge) for (size_t idx = start[edge]; idx < stop[edge]; idx += stride[edge]) { if (landmask[idx]) { - neumann[edge].push_back(idx); + dirichlet[edge].push_back(idx); } } - for (Edge edge : edges) - std::sort(neumann[edge].begin(), neumann[edge].end()); + sortDirichlet(edge); } /*! diff --git a/dynamics/src/include/BrittleCGDynamicsKernel.hpp b/dynamics/src/include/BrittleCGDynamicsKernel.hpp index fb20e0b94..c1ba25b71 100644 --- a/dynamics/src/include/BrittleCGDynamicsKernel.hpp +++ b/dynamics/src/include/BrittleCGDynamicsKernel.hpp @@ -1,7 +1,7 @@ /*! * @file BrittleCGDynamicsKernel.hpp * - * @date 21 Feb 2025 + * @date 06 Dec 2024 * @author Tim Spain * @author Einar Ólason */ @@ -12,6 +12,7 @@ #include "CGDynamicsKernel.hpp" #include "BBMParameters.hpp" +#include "ParametricMap.hpp" #include "StressUpdateStep.hpp" #include "include/constants.hpp" #include @@ -36,7 +37,6 @@ template class BrittleCGDynamicsKernel : public CGDynamicsKern using DynamicsKernel::applyBoundaries; using DynamicsKernel::advectionAndLimits; using DynamicsKernel::dgtransport; - using DynamicsKernel::neumannZero; using CGDynamicsKernel::u; using CGDynamicsKernel::v; @@ -109,8 +109,6 @@ template class BrittleCGDynamicsKernel : public CGDynamicsKern Nextsim::LimitMax(damage, 1.0); Nextsim::LimitMin(damage, 1e-12); - neumannZero(damage); - prepareIteration({ { hiceName, hice }, { ciceName, cice } }); // The timestep for the brittle solver is the solver subtimestep diff --git a/dynamics/src/include/CGDynamicsKernel.hpp b/dynamics/src/include/CGDynamicsKernel.hpp index c5f902cec..554f0c713 100644 --- a/dynamics/src/include/CGDynamicsKernel.hpp +++ b/dynamics/src/include/CGDynamicsKernel.hpp @@ -1,7 +1,7 @@ /*! * @file CGDynamicsKernel.hpp * - * @date 21 Feb 2025 + * @date 06 Dec 2024 * @author Tim Spain */ @@ -69,7 +69,6 @@ class CGDynamicsKernel : public DynamicsKernel { protected: void addStressTensorCell(const size_t eid, const size_t cx, const size_t cy); void dirichletZero(CGVector&) const; - void neumannZero(CGVector&) const; // CG ice velocity CGVector u; CGVector v; diff --git a/dynamics/src/include/DynamicsKernel.hpp b/dynamics/src/include/DynamicsKernel.hpp index 446f7c80a..08f9ecc19 100644 --- a/dynamics/src/include/DynamicsKernel.hpp +++ b/dynamics/src/include/DynamicsKernel.hpp @@ -1,7 +1,7 @@ /*! * @file DynamicsKernel.hpp * - * @date 21 Feb 2025 + * @date Jan 5, 2024 * @author Tim Spain */ @@ -53,9 +53,9 @@ template class DynamicsKernel { smesh->RotatePoleToGreenland(); smesh->landmaskFromModelArray(mask); smesh->dirichletFromMask(); - // TODO: handle periodic edges - for (const ParametricMesh::Edge edge : ParametricMesh::edges) { - smesh->neumannFromEdge(edge); + // TODO: handle periodic and open edges + for (ParametricMesh::Edge edge : ParametricMesh::edges) { + smesh->dirichletFromEdge(edge); } //! Initialize transport @@ -180,9 +180,6 @@ template class DynamicsKernel { Nextsim::LimitMax(cice, 1.0); Nextsim::LimitMin(cice, 0.0); Nextsim::LimitMin(hice, 0.0); - - neumannZero(cice); - neumannZero(hice); } protected: @@ -235,39 +232,6 @@ template class DynamicsKernel { */ virtual void prepareAdvection() = 0; - template void neumannZero(DGVector& v) const - { - // the four segments bottom, right, top, left, are each processed in parallel - for (size_t seg = 0; seg < 4; ++seg) { -#pragma omp parallel for - for (size_t i = 0; i < smesh->neumann[seg].size(); ++i) { - - const size_t eid = smesh->neumann[seg][i]; - const size_t ix = eid % smesh->nx; // compute coordinates of element - const size_t iy = eid / smesh->nx; - - switch (seg) { - case 0: // bottom <= top - for (size_t j = 0; j < DG; ++j) - v(iy * smesh->nx + ix, j) = v((iy + 1) * smesh->nx + ix, j); - break; - case 1: // right <= left - for (size_t j = 0; j < DG; ++j) - v(iy * smesh->nx + ix, j) = v(iy * smesh->nx + ix - 1, j); - break; - case 2: // top <= bottom - for (size_t j = 0; j < DG; ++j) - v(iy * smesh->nx + ix, j) = v((iy - 1) * smesh->nx + ix, j); - break; - case 3: // left <= right - for (size_t j = 0; j < DG; ++j) - v(iy * smesh->nx + ix, j) = v(iy * smesh->nx + ix + 1, j); - break; - } - } - } - } - private: std::unordered_map> advectedFields; diff --git a/dynamics/src/include/ParametricMesh.hpp b/dynamics/src/include/ParametricMesh.hpp index c34c8f559..25e566cf1 100644 --- a/dynamics/src/include/ParametricMesh.hpp +++ b/dynamics/src/include/ParametricMesh.hpp @@ -1,6 +1,6 @@ /*! * @file ParametricMesh.hpp - * @date 21 Feb 2025 + * @date 09 Jan 2025 * @author Thomas Richter */ @@ -56,9 +56,9 @@ class ParametricMesh { Eigen::Matrix vertices; // stores the /*! - * Stores the Dirichlet and Neumann boundary information + * Stores the Dirichlet boundary information * - * Dirichlet- and Neumann-Information is stored in 4 vectors, e.g. + * Dirichlet-Information is stored in 4 vectors * dirichlet[0]: List of elements that have a Dirichlet-Boundary on the bottom * dirichlet[1]: List of elements that have a Dirichlet-Boundary on the right * dirichlet[2]: List of elements that have a Dirichlet-Boundary on the top @@ -67,7 +67,6 @@ class ParametricMesh { * Each of the vectors is processed in parallel */ std::array, 4> dirichlet; - std::array, 4> neumann; /*! * Periodic: @@ -404,7 +403,7 @@ class ParametricMesh { * * @param edge index of the edge to add closed boundary conditions to. */ - void neumannFromEdge(Edge edge); + void dirichletFromEdge(Edge edge); /*! * Sort all the dirichlet arrays, so the element indices are ordered. diff --git a/dynamics/test/ParametricMesh_test.cpp b/dynamics/test/ParametricMesh_test.cpp index a16321636..28d1d243d 100644 --- a/dynamics/test/ParametricMesh_test.cpp +++ b/dynamics/test/ParametricMesh_test.cpp @@ -3,7 +3,7 @@ * * @brief Test the ParametricMesh class, especially processing from ModelArray files. * - * @date 21 Feb 2025 + * @date Dec 15, 2023 * @author Tim Spain */ @@ -122,7 +122,7 @@ TEST_CASE("Compare readmesh and landmask reading") // Dirichlet conditions fromArrays.dirichletFromMask(); for (ParametricMesh::Edge edge : ParametricMesh::edges) { - fromArrays.neumannFromEdge(edge); + fromArrays.dirichletFromEdge(edge); } fromArrays.sortDirichlet(); for (ParametricMesh::Edge edge : ParametricMesh::edges) { From 0077c28823ec6ccc35f0ea053cde6753b0530375 Mon Sep 17 00:00:00 2001 From: Thomas Richter Date: Mon, 7 Apr 2025 08:40:25 +0200 Subject: [PATCH 06/33] New boundary handling. Dirichlet: all done using the cglandmask field Neumann: can only appear on the outer edges of the mesh, if the last elements are ice elements. --- dynamics/src/CGDynamicsKernel.cpp | 39 ++------- dynamics/src/DGTransport.cpp | 95 +++++++++++---------- dynamics/src/include/VPCGDynamicsKernel.hpp | 13 +++ run/config_polynya.cfg | 5 +- run/make_init_polynya.py | 9 +- 5 files changed, 78 insertions(+), 83 deletions(-) diff --git a/dynamics/src/CGDynamicsKernel.cpp b/dynamics/src/CGDynamicsKernel.cpp index 83726fc9d..d261a86db 100644 --- a/dynamics/src/CGDynamicsKernel.cpp +++ b/dynamics/src/CGDynamicsKernel.cpp @@ -383,40 +383,13 @@ template void CGDynamicsKernel::stressDivergence( template void CGDynamicsKernel::dirichletZero(CGVector& v) const { - // the four segments bottom, right, top, left, are each processed in parallel - for (size_t seg = 0; seg < 4; ++seg) { + // TR 07.04.2025: Dirichlet Zero (u=v=0) holds on land and on the boundary betreen + // land and ice. Hence on all elements with landmask = 0, or, on cg nodes with + // cglandmask = 0 #pragma omp parallel for - for (size_t i = 0; i < smesh->dirichlet[seg].size(); ++i) { - - const size_t eid = smesh->dirichlet[seg][i]; - const size_t ix = eid % smesh->nx; // compute coordinates of element - const size_t iy = eid / smesh->nx; - - if (seg == 0) // bottom - for (size_t j = 0; j < CGdegree + 1; ++j) - v(iy * CGdegree * (CGdegree * smesh->nx + 1) + CGdegree * ix + j, 0) = 0.0; - else if (seg == 1) // right - for (size_t j = 0; j < CGdegree + 1; ++j) - v(iy * CGdegree * (CGdegree * smesh->nx + 1) + CGdegree * ix + CGdegree - + (CGdegree * smesh->nx + 1) * j, - 0) - = 0.0; - else if (seg == 2) // top - for (size_t j = 0; j < CGdegree + 1; ++j) - v((iy + 1) * CGdegree * (CGdegree * smesh->nx + 1) + CGdegree * ix + j, 0) - = 0.0; - else if (seg == 3) // left - for (size_t j = 0; j < CGdegree + 1; ++j) - v(iy * CGdegree * (CGdegree * smesh->nx + 1) + CGdegree * ix - + (CGdegree * smesh->nx + 1) * j, - 0) - = 0.0; - else { - std::cerr << "That should not have happened!" << std::endl; - abort(); - } - } - } + for (size_t i = 0; i < pmap->cglandmask.rows(); ++i) + if (pmap->cglandmask(i) == 0) // land + v(i) = 0; } template void CGDynamicsKernel::applyBoundaries() diff --git a/dynamics/src/DGTransport.cpp b/dynamics/src/DGTransport.cpp index b42294c75..62cc8d13c 100644 --- a/dynamics/src/DGTransport.cpp +++ b/dynamics/src/DGTransport.cpp @@ -227,28 +227,23 @@ template void DGTransport::reinitnormalvelocity() } } - // Take care of the boundaries. Usually, the normal velocity is the average velocity - // from left and from the right. Hence, we get the factor 0.5 above. At boundaries, - // the normal is set only once, from the inside. These edges must be scaled with 2.0 - - for (size_t seg = 0; seg < 4; ++seg) // run over the 4 segments (bot, right, top, left) - { -#pragma omp parallel for - for (size_t i = 0; i < smesh.dirichlet[seg].size(); ++i) { - const size_t eid = smesh.dirichlet[seg][i]; //! The i of the boundary element - const size_t ix = eid % smesh.nx; //! x & y indices of the element - const size_t iy = eid / smesh.nx; - - if (seg == 0) // bottom - normalvel_X.row(smesh.nx * iy + ix) *= 2.0; - else if (seg == 1) // right - normalvel_Y.row((smesh.nx + 1) * iy + ix + 1) *= 2.0; - else if (seg == 2) // top - normalvel_X.row(smesh.nx * (iy + 1) + ix) *= 2.0; - else if (seg == 3) // left - normalvel_Y.row((smesh.nx + 1) * iy + ix) *= 2.0; - } - } + // TR 07.04.2025 + // correction of the normal velocity is needed on the mesh-boundary, where each edge + // only has one neighbor. Above, we sum both sides, weighted with 0.5. Therefore, + // we must scale the outer edges with 2. This is only used for inflow/outflow + // Neumann boundaries. +#pragma omp parallel for + for (size_t i=0;i::DGTransportOperator(const ParametricMesh& smesh, const dou } } - // Dirichlet - for (size_t seg = 0; seg < 4; ++seg) // run over the 4 segments (bot, right, top, left) - { + // TR 07.04.2025 + // + // Dirichlet: we do not need to do anything special in DGTransport. On Dirichlet + // boundaries, the normal velocity is zero. Hence, any flux is zero. + + // + // Neumann inflow / outflow + // This can only appear on outer mesh-boundaries, whenever the element at the edge + // is ice. #pragma omp parallel for - for (size_t i = 0; i < smesh.dirichlet[seg].size(); ++i) { - const size_t eid = smesh.dirichlet[seg][i]; - const size_t ix = eid % smesh.nx; // compute 'coordinate' of element - const size_t iy = eid / smesh.nx; - - if (seg == 0) // bottom - boundary_lower(smesh, dt, phiup, phi, normalvel_X, eid, smesh.nx * iy + ix); - else if (seg == 1) // right - boundary_right( - smesh, dt, phiup, phi, normalvel_Y, eid, (smesh.nx + 1) * iy + ix + 1); - else if (seg == 2) // top - boundary_upper(smesh, dt, phiup, phi, normalvel_X, eid, smesh.nx * (iy + 1) + ix); - else if (seg == 3) // left - boundary_left(smesh, dt, phiup, phi, normalvel_Y, eid, (smesh.nx + 1) * iy + ix); - else { - std::cerr << "Wrong Dirichlet boundary information in the mesh. Boundary side " - << smesh.dirichlet[seg][i] << " not valid" << std::endl; - abort(); - } - } - } + for (size_t i=0; i class VPCGDynamicsKernel : public CGDynamicsKernel::stepNumber%si==0) + { + int vtkn = DynamicsKernel::stepNumber/si; + VTK::write_dg("Polynya/hice", vtkn, hice, *smesh); + VTK::write_dg("Polynya/cice", vtkn, cice ,*smesh); + VTK::write_cg_velocity("Polynya/vel",vtkn, u,v,*smesh); + VTK::write_cg("Polynya/cgl",vtkn, pmap->cglandmask,*smesh); + } + + + // Let DynamicsKernel handle the advection step DynamicsKernel::advectionAndLimits(tst); diff --git a/run/config_polynya.cfg b/run/config_polynya.cfg index a65c97292..61fcdb6d6 100644 --- a/run/config_polynya.cfg +++ b/run/config_polynya.cfg @@ -2,7 +2,7 @@ init_file = init_polynya.nc start = 2023-01-01T00:00:00Z stop = 2023-01-03T00:00:00Z -time_step = P0-0T00:02:00 +time_step = P0-0T00:01:00 missing_value = 1e20 [Modules] @@ -10,7 +10,8 @@ DiagnosticOutputModule = Nextsim::ConfigOutput DynamicsModule = Nextsim::MEVPDynamics AtmosphereBoundaryModule = Nextsim::FluxConfiguredAtmosphere OceanBoundaryModule = Nextsim::FluxConfiguredOcean -IceThermodynamicsModule = Nextsim::ThermoWinton +#IceThermodynamicsModule = Nextsim::ThermoWinton +IceThermodynamicsModule = Nextsim::DummyIceThermodynamics [ConfigOutput] period = P0-0T04:00:00 diff --git a/run/make_init_polynya.py b/run/make_init_polynya.py index 8317ac353..852769243 100644 --- a/run/make_init_polynya.py +++ b/run/make_init_polynya.py @@ -3,10 +3,14 @@ # Creates initial conditions for the Bjornsson et al. (2001) polynya case # Domain size [km] -x = 100 -y = 50 +#x = 100 +#y = 50 +#res = 2 +x = 64 +y = 48 res = 2 + nfirst = int(y / res) nsecond = int(x / res) nLayers = 3 @@ -21,6 +25,7 @@ initializer.mask[0, :] = 0. initializer.mask[-1, :] = 0. initializer.mask[:, 0] = 0. +#initializer.mask[:, -1] = 0. ## right # Uniform concentration of 90% initializer.cice[:, :] = 0.9 From 02720e126d394a5a27376b0374f7ff72f4a9463a Mon Sep 17 00:00:00 2001 From: Thomas Richter Date: Mon, 7 Apr 2025 16:25:00 +0200 Subject: [PATCH 07/33] - Removed old Dirichlet-Info from the mesh. - fixed advection and mesh test cases --- dynamics/src/DGTransport.cpp | 6 + dynamics/src/ParametricMesh.cpp | 104 +--- dynamics/src/cgParametricMomentum.cpp | 557 ------------------ dynamics/src/include/DynamicsKernel.hpp | 5 - dynamics/src/include/ParametricMesh.hpp | 38 -- dynamics/src/include/cgParametricMomentum.hpp | 285 --------- dynamics/test/AdvectionPeriodicBC_test.cpp | 49 +- dynamics/test/Advection_test.cpp | 89 +-- dynamics/test/ParametricMesh_test.cpp | 28 - 9 files changed, 90 insertions(+), 1071 deletions(-) delete mode 100644 dynamics/src/cgParametricMomentum.cpp delete mode 100644 dynamics/src/include/cgParametricMomentum.hpp diff --git a/dynamics/src/DGTransport.cpp b/dynamics/src/DGTransport.cpp index 62cc8d13c..5e2b0ac64 100644 --- a/dynamics/src/DGTransport.cpp +++ b/dynamics/src/DGTransport.cpp @@ -484,6 +484,11 @@ void DGTransport::DGTransportOperator(const ParametricMesh& smesh, const dou // Neumann inflow / outflow // This can only appear on outer mesh-boundaries, whenever the element at the edge // is ice. + + // TO DISCUSS: there can't be in/outflow and periodic on the same edge. We either need a structure + // for in/outflow or we say: either there is an in/outflow boundary or periodic. But not both.. + if (smesh.periodic.size()==0) + { #pragma omp parallel for for (size_t i=0; i::DGTransportOperator(const ParametricMesh& smesh, const dou if (smesh.landmask[element_right] == 1) boundary_right(smesh, dt, phiup, phi, normalvel_Y, element_right, edge_right); } + } #pragma omp parallel for for (size_t eid = 0; eid < smesh.nelements; ++eid) diff --git a/dynamics/src/ParametricMesh.cpp b/dynamics/src/ParametricMesh.cpp index 1bf15bd2d..a38af8576 100644 --- a/dynamics/src/ParametricMesh.cpp +++ b/dynamics/src/ParametricMesh.cpp @@ -70,18 +70,6 @@ void ParametricMesh::readmesh(std::string fname) // Boundary if (version == "1.0") // all four boundaries are dirichlet { - for (size_t i = 0; i < nx; ++i) // lower - dirichlet[0].push_back(i); - - for (size_t i = 0; i < ny; ++i) // right - dirichlet[1].push_back(i * nx + nx - 1); - - for (size_t i = 0; i < nx; ++i) // upper - dirichlet[2].push_back(i + nx * (ny - 1)); - - for (size_t i = 0; i < ny; ++i) // left - dirichlet[3].push_back(i * nx); - landmask.resize(nx * ny, true); // set landmask } else if (version == "2.0") // landmask, dirichlet and periodic boundary as additional lists { @@ -128,25 +116,14 @@ void ParametricMesh::readmesh(std::string fname) << std::endl; abort(); } + std::cerr << "ParametricMesh V2.0 is not longer fully supported. Dirichlet now uses the landmask information and Dirichlet data within the mesh file is ignored" << std::endl; size_t nd; IN >> nd; - - if (statuslog > 0) - std::cout << "reading " << nd << " dirichlet segments" << std::endl; - - for (size_t i = 0; i < nd; ++i) { - size_t n0, n1; - IN >> n0 >> n1; // read the element and the side - assert(n0 < nx * ny); - assert(n0 >= 0); - - dirichlet[n1].push_back(n0); - - if (IN.eof()) { - std::cerr << "ParametricMesh :: Unexpected eof << " << fname << std::endl; - abort(); - } - } + double tmp; + for (int i=0;i> tmp >> tmp; + + IN >> status; if (status != "periodic") { @@ -229,75 +206,6 @@ void ParametricMesh::landmaskFromModelArray(const ModelArray& mask) } } -/*! - * Add to the dirichlet arrays according to the stored landmask. - */ -void ParametricMesh::dirichletFromMask() -{ - // Edges are accessed in the order: BOTTOM, RIGHT, TOP, LEFT. See also ParametricMesh::edges. - const std::array startX = { 0, 0, 0, 1 }; - const std::array stopX = { nx, nx - 1, nx, nx }; - const std::array startY = { 1, 0, 0, 0 }; - const std::array stopY = { ny, ny, ny - 1, ny }; - const std::array deltaIdx = { -static_cast(nx), 1, static_cast(nx), -1 }; - - // Loop over edges - for (Edge edge : edges) { - for (size_t j = startY[edge]; j < stopY[edge]; ++j) { - for (size_t i = startX[edge]; i < stopX[edge]; ++i) { - size_t idx = ModelArray::indexFromLocation(ModelArray::Type::H, { i, j }); - if (!landmask[idx]) - continue; - // mask(i, j) is ocean. Check the appropriate neighbour - if (!landmask[idx + deltaIdx[edge]]) { - dirichlet[edge].push_back(idx); - } - } - } - sortDirichlet(edge); - } -} - -/*! - * Add to the dirichlet arrays due to the domain edges according to an edge index. - * - * @param edge index of the edge to add closed boundary conditions to. - */ -void ParametricMesh::dirichletFromEdge(Edge edge) -{ - // BOTTOM, RIGHT, TOP, LEFT - const std::array start = { 0, nx - 1, nelements - nx, 0 }; - const std::array stop = { nx, nelements, nelements, nelements }; - const std::array stride = { 1, nx, 1, nx }; - - for (size_t idx = start[edge]; idx < stop[edge]; idx += stride[edge]) { - if (landmask[idx]) { - dirichlet[edge].push_back(idx); - } - } - sortDirichlet(edge); -} - -/*! - * Sort all the dirichlet arrays, so the element indices are ordered. - */ -void ParametricMesh::sortDirichlet() -{ - for (ParametricMesh::Edge edge : ParametricMesh::edges) { - sortDirichlet(edge); - } -} - -/*! - * Sort the dirichlet array of one particular edge. - * - * @param edge the edge to be sorted. - */ -void ParametricMesh::sortDirichlet(Edge edge) -{ - std::sort(dirichlet[edge].begin(), dirichlet[edge].end()); -} - /*! * returns minimum mesh size. * diff --git a/dynamics/src/cgParametricMomentum.cpp b/dynamics/src/cgParametricMomentum.cpp deleted file mode 100644 index 9bd465c88..000000000 --- a/dynamics/src/cgParametricMomentum.cpp +++ /dev/null @@ -1,557 +0,0 @@ -/*! - * @file ParametricMomentum.cpp - * @date 06 Dec 2024 - * @author Thomas Richter - */ - -#include "cgParametricMomentum.hpp" -#include "BBM.hpp" -#include "Interpolations.hpp" -#include "MEB.hpp" -#include "ParametricTools.hpp" -#include "VectorManipulations.hpp" -#include "dgVisu.hpp" -#include "mevp.hpp" - -namespace Nextsim { - -#define DGSTRESS(CG) ((CG == 1 ? 3 : (CG == 2 ? 8 : -1))) - -////////////////////////////////////////////////// Strain (ParametricMesh) - -template void CGParametricMomentum::ProjectCGVelocityToDGStrain() -{ - // !!! must still be converted to the spherical system!!! - - assert(static_cast((CG * smesh.nx + 1) * (CG * smesh.ny + 1)) == vx.rows()); - assert(static_cast((CG * smesh.nx + 1) * (CG * smesh.ny + 1)) == vy.rows()); - assert(static_cast(smesh.nx * smesh.ny) == E11.rows()); - assert(static_cast(smesh.nx * smesh.ny) == E12.rows()); - assert(static_cast(smesh.nx * smesh.ny) == E22.rows()); - - const int cgshift = CG * smesh.nx + 1; //!< Index shift for each row - - // parallelize over the rows -#pragma omp parallel for - for (size_t row = 0; row < smesh.ny; ++row) { - int dgi = smesh.nx * row; //!< Index of dg vector - int cgi = CG * cgshift * row; //!< Lower left index of cg vector - - for (size_t col = 0; col < smesh.nx; ++col, ++dgi, cgi += CG) { // loop over all elements - - if (smesh.landmask[dgi] == 0) // only on ice - continue; - - // get the 4 (cg1) 9 (cg2) local x/y - velocity coefficients on the element - Eigen::Matrix vx_local, vy_local; - if (CG == 1) { - vx_local << vx(cgi), vx(cgi + 1), vx(cgi + cgshift), vx(cgi + 1 + cgshift); - vy_local << vy(cgi), vy(cgi + 1), vy(cgi + cgshift), vy(cgi + 1 + cgshift); - } else if (CG == 2) { - vx_local << vx(cgi), vx(cgi + 1), vx(cgi + 2), vx(cgi + cgshift), - vx(cgi + 1 + cgshift), vx(cgi + 2 + cgshift), vx(cgi + 2 * cgshift), - vx(cgi + 1 + 2 * cgshift), vx(cgi + 2 + 2 * cgshift); - - vy_local << vy(cgi), vy(cgi + 1), vy(cgi + 2), vy(cgi + cgshift), - vy(cgi + 1 + cgshift), vy(cgi + 2 + cgshift), vy(cgi + 2 * cgshift), - vy(cgi + 1 + 2 * cgshift), vy(cgi + 2 + 2 * cgshift); - } else - abort(); - - // Solve (E, Psi) = (0.5(DV + DV^T), Psi) - // by integrating rhs and inverting with dG(stress) mass matrix - // - E11.row(dgi) = pmap.iMgradX[dgi] * vx_local; - E22.row(dgi) = pmap.iMgradY[dgi] * vy_local; - E12.row(dgi) = 0.5 * (pmap.iMgradX[dgi] * vy_local + pmap.iMgradY[dgi] * vx_local); - - if (smesh.CoordinateSystem == SPHERICAL) { - E11.row(dgi) -= pmap.iMM[dgi] * vy_local; - E12.row(dgi) += 0.5 * pmap.iMM[dgi] * vx_local; - } - } - } -} - -////////////////////////////////////////////////// STRESS Tensor -// Sasip-Mesh Interface -template -void CGParametricMomentum::DivergenceOfStress( - const double scale, CGVector& tx, CGVector& ty) const -{ -#pragma omp parallel for - for (size_t i = 0; i < tx.rows(); ++i) { - tx(i) = 0.0; - ty(i) = 0.0; - } - - // parallelization in stripes - for (size_t p = 0; p < 2; ++p) -#pragma omp parallel for schedule(static) - for (size_t cy = 0; cy < smesh.ny; ++cy) //!< loop over all cells of the mesh - { - if (cy % 2 == p) { - size_t c = smesh.nx * cy; - for (size_t cx = 0; cx < smesh.nx; ++cx, ++c) //!< loop over all cells of the mesh - if (smesh.landmask[c] == 1) // only on ice! - AddStressTensorCell(scale, c, cx, cy, tx, ty); - } - } - // set zero on the Dirichlet boundaries - DirichletZero(tx); - DirichletZero(ty); - // add the contributions on the periodic boundaries - VectorManipulations::CGAveragePeriodic(smesh, tx); - VectorManipulations::CGAveragePeriodic(smesh, ty); -} - -template void CGParametricMomentum::DirichletZero(CGVector& v) const -{ - // the four segments bottom, right, top, left, are each processed in parallel - for (size_t seg = 0; seg < 4; ++seg) { -#pragma omp parallel for - for (size_t i = 0; i < smesh.dirichlet[seg].size(); ++i) { - - const size_t eid = smesh.dirichlet[seg][i]; - const size_t ix = eid % smesh.nx; // compute 'coordinate' of element - const size_t iy = eid / smesh.nx; - - if (seg == 0) // bottom - for (size_t j = 0; j < CG + 1; ++j) - v(iy * CG * (CG * smesh.nx + 1) + CG * ix + j, 0) = 0.0; - else if (seg == 1) // right - for (size_t j = 0; j < CG + 1; ++j) - v(iy * CG * (CG * smesh.nx + 1) + CG * ix + CG + (CG * smesh.nx + 1) * j, 0) - = 0.0; - else if (seg == 2) // top - for (size_t j = 0; j < CG + 1; ++j) - v((iy + 1) * CG * (CG * smesh.nx + 1) + CG * ix + j, 0) = 0.0; - else if (seg == 3) // left - for (size_t j = 0; j < CG + 1; ++j) - v(iy * CG * (CG * smesh.nx + 1) + CG * ix + (CG * smesh.nx + 1) * j, 0) = 0.0; - else { - std::cerr << "That should not have happened!" << std::endl; - abort(); - } - } - } -} - -template void CGParametricMomentum::CheckPeriodicity(CGVector& v) -{ - // the two segments bottom, right, top, left, are each processed in parallel - for (size_t seg = 0; seg < smesh.periodic.size(); ++seg) { - // #pragma omp parallel for - for (size_t i = 0; i < smesh.periodic[seg].size(); ++i) { - - const size_t ptype = smesh.periodic[seg][i][0]; - const size_t eid_lb = smesh.periodic[seg][i][2]; - const size_t eid_rt = smesh.periodic[seg][i][1]; - - size_t ix_lb = eid_lb % smesh.nx; - size_t iy_lb = eid_lb / smesh.nx; - size_t i0_lb = (CG * smesh.nx + 1) * CG * iy_lb - + CG * ix_lb; // lower/left index in left/bottom element - size_t ix_rt = eid_rt % smesh.nx; - size_t iy_rt = eid_rt / smesh.nx; - size_t i0_rt = (CG * smesh.nx + 1) * CG * iy_rt - + CG * ix_rt; // lower/left index in right/top element - - std::cout << std::setprecision(16); - if (ptype == 0) // X-edge, bottom/top - { - for (size_t j = 0; j <= CG; ++j) { - double check = v(i0_lb + j) - v(i0_rt + CG * (CG * smesh.nx + 1) + j); - if (fabs(check) > 1.e-13) - std::cout << v(i0_lb + j) - v(i0_rt + CG * (CG * smesh.nx + 1) + j) - << std::endl; - } - } else if (ptype == 1) // Y-edge, left/right - { - for (size_t j = 0; j < CG; ++j) { - double check = v(i0_lb + j * (CG * smesh.nx + 1)) - - v(i0_rt + CG + j * (CG * smesh.nx + 1)); - if (fabs(check) > 1.e-13) - std::cout << v(i0_lb + j * (CG * smesh.nx + 1)) - - v(i0_rt + CG + j * (CG * smesh.nx + 1)) - << std::endl; - } - } else - abort(); - } - } -} - -// -------------------------------------------------- - -template -template -void CGParametricMomentum::prepareIteration(const DGVector& H, const DGVector& A) -{ - // copy old velocity - vx_mevp = vx; - vy_mevp = vy; - // interpolate ice height and concentration to local cg Variables - Interpolations::DG2CG(smesh, cg_A, A); - VectorManipulations::CGAveragePeriodic(smesh, cg_A); - Interpolations::DG2CG(smesh, cg_H, H); - VectorManipulations::CGAveragePeriodic(smesh, cg_H); - - // limit A to [0,1] and H to [0, ...) - cg_A = cg_A.cwiseMin(1.0); - cg_A = cg_A.cwiseMax(1.e-4); - cg_H = cg_H.cwiseMax(1.e-4); -} - -template -template -void CGParametricMomentum::mEVPStep(const VPParameters& params, const size_t NT_evp, - const double alpha, const double beta, double dt_adv, const DGVector& H, - const DGVector& A) -{ - - // Compute Strain Rate - - ProjectCGVelocityToDGStrain(); - - // Update the stresses according to the mEVP model - - std::cerr << "CGParametricMomentum::mEVPStep(): Fatal: StressUpdateHighOrder requires " - "ParametricMomentumMap with proper template parameter. See header file and use " - "kernel infrastructure." - << std::endl; - abort(); - /* NB! For this function to work, StressUpdateHighOrder must be called with the proper template - * parameters. This would require some work, but cgParametricMomentum is depricated and will be - * removed in future version. */ - // Nextsim::mEVP::StressUpdateHighOrder( - // params, pmap, smesh, S11, S12, S22, E11, E12, E22, H, A, alpha, beta); - - // Compute the divergence of the stress tensor - - double stressscale = 1.0; // 2nd-order Stress term has different scaling with the EarthRadius - if (smesh.CoordinateSystem == Nextsim::SPHERICAL) - stressscale = 1.0 / Nextsim::EarthRadius / Nextsim::EarthRadius; - - DivergenceOfStress(stressscale, tmpx, tmpy); // Compute divergence of stress tensor - - // Update the velocity - double SC = 1.0; ///(1.0-pow(1.0+1.0/beta,-1.0*NT_evp)); - - const double FOcean = params.COcean * params.rhoOcean; - const double FAtm = params.CAtm * params.rhoAtm; - - // update by a loop.. implicit parts and h-dependent -#pragma omp parallel for - for (int i = 0; i < vx.rows(); ++i) { - double absatm = sqrt(ax(i) * ax(i) + ay(i) * ay(i)); - double absocn = sqrt(SQR(vx(i) - ox(i)) + SQR(vy(i) - oy(i))); - - vx(i) = (1.0 - / (params.rhoIce * cg_H(i) / dt_adv * (1.0 + beta) // implicit parts - + cg_A(i) * FOcean * absocn) // implicit parts - * (params.rhoIce * cg_H(i) / dt_adv * (beta * vx(i) + vx_mevp(i)) // pseudo-timestepping - + cg_A(i) - * (FAtm * absatm * ax(i) + // atm forcing - FOcean * absocn * SC * ox(i)) // ocean forcing - + params.rhoIce * cg_H(i) * params.fc * (vy(i) - oy(i)) // cor + surface - + tmpx(i) / pmap.lumpedcgmass(i))); - vy(i) = (1.0 - / (params.rhoIce * cg_H(i) / dt_adv * (1.0 + beta) // implicit parts - + cg_A(i) * FOcean * absocn) // implicit parts - * (params.rhoIce * cg_H(i) / dt_adv * (beta * vy(i) + vy_mevp(i)) // pseudo-timestepping - + cg_A(i) - * (FAtm * absatm * ay(i) + // atm forcing - FOcean * absocn * SC * oy(i)) // ocean forcing - + params.rhoIce * cg_H(i) * params.fc * (ox(i) - vx(i)) - + tmpy(i) / pmap.lumpedcgmass(i))); - } - - DirichletZero(); - - // Landmask - const size_t inrow = CG * smesh.nx + 1; -#pragma omp parallel for - for (size_t eid = 0; eid < smesh.nelements; ++eid) - if (smesh.landmask[eid] == 0) { - const size_t ex = eid % smesh.nx; - const size_t ey = eid / smesh.nx; - for (int jy = 0; jy < CG + 1; ++jy) - for (int jx = 0; jx < CG + 1; ++jx) { - vx(inrow * (CG * ey + jy) + CG * ex + jx) = 0.0; - vy(inrow * (CG * ey + jy) + CG * ex + jx) = 0.0; - } - } -} -// -------------------------------------------------- -template -template -void CGParametricMomentum::prepareIteration( - const DGVector& H, const DGVector& A, const DGVector& D) -{ - - // set the average sub-iteration velocity to zero - avg_vx.setZero(); - avg_vy.setZero(); - - // interpolate ice height and concentration to local cg Variables - Interpolations::DG2CG(smesh, cg_A, A); - VectorManipulations::CGAveragePeriodic(smesh, cg_A); - Interpolations::DG2CG(smesh, cg_H, H); - VectorManipulations::CGAveragePeriodic(smesh, cg_H); - Interpolations::DG2CG(smesh, cg_D, D); - VectorManipulations::CGAveragePeriodic(smesh, cg_D); - - // limit A and D to [0,1] and H to [0, ...) - cg_A = cg_A.cwiseMin(1.0); - cg_A = cg_A.cwiseMax(1.e-4); - cg_H = cg_H.cwiseMax(1.e-4); - cg_D = cg_D.cwiseMin(1.0); - cg_D = cg_D.cwiseMax(0.0); -} - -/* This is Hunke and Dukowicz's solution to (22), multiplied - * with (\Delta t/m)^2 to ensure stability for c' = 0 - * - * This scheme includes Coriolis terms in an implicit way - */ -template -template -void CGParametricMomentum::BBMStep(const BBMParameters& params, const size_t NT_meb, - double dt_adv, const DGVector& H, const DGVector& A, DGVector& D) -{ - - double dt_mom = dt_adv / NT_meb; - - //! Compute Strain Rate - ProjectCGVelocityToDGStrain(); - - // TODO compute stress update with precomputed transformations - Nextsim::MEB::StressUpdateHighOrder( - params, smesh, S11, S12, S22, E11, E12, E22, H, A, D, dt_mom); - // Nextsim::MEB::StressUpdateHighOrder(params, ptrans, smesh, S11, S12, S22, E11, E12, E22, H, - // A, D, dt_mom); - - // Divergence of the stress tensor -#pragma omp parallel for - for (int i = 0; i < tmpx.rows(); ++i) - tmpx(i) = tmpy(i) = 0; - - // AddStressTensor(ptrans_stress, -1.0, tmpx, tmpy, S11, S12, S22); - // AddStressTensor(-1.0, tmpx, tmpy); - // Check first parameter, this scaling - DivergenceOfStress(1.0, tmpx, tmpy); // Compute divergence of stress tensor - - // FIXME: We're missing the gradient of the sea-surface slope (\nabla \eta) - - /* This is Hunke and Dukowicz's solution to (22), multiplied - * with (\Delta t/m)^2 to ensure stability for c' = 0 */ - double const cos_ocean_turning_angle = std::cos(params.oceanTurningAngle * M_PI / 180.); - double const sin_ocean_turning_angle = std::sin(params.oceanTurningAngle * M_PI / 180.); - - const double FOcean = params.COcean * params.rhoOcean; - const double FAtm = params.CAtm * params.rhoAtm; - -#pragma omp parallel for - for (int i = 0; i < vx.rows(); ++i) { - // FIXME: dte_over_mass should include snow (total mass) - double const dte_over_mass = dt_mom / (params.rhoIce * cg_H(i)); - double const uice = vx(i); - double const vice = vy(i); - - double const c_prime = cg_A(i) * FOcean * std::hypot(ox(i) - uice, oy(i) - vice); - - // FIXME: Need the grounding term: tau_b = C_bu[i]/(std::hypot(uice,vice)+u0); - double const tau_b = 0.; - double const alpha = 1. + dte_over_mass * (c_prime * cos_ocean_turning_angle + tau_b); - /* FIXME: We need latitude here. Then this becomes: - * double const beta = dt_mom*fc + - * dte_over_mass*c_prime*std::copysign(sin_ocean_turning_angle, lat[i]); */ - double const beta = dt_mom * params.fc + dte_over_mass * c_prime * sin_ocean_turning_angle; - double const rdenom = 1. / (alpha * alpha + beta * beta); - - double const drag_atm = cg_A(i) * FAtm * std::hypot(ax(i), ay(i)); - double const tau_x = drag_atm * ax(i) - + c_prime * (ox(i) * cos_ocean_turning_angle - oy(i) * sin_ocean_turning_angle); - /* FIXME: Need latitude here. Then This becomes: - * + c_prime*( ox(i)*cos_ocean_turning_angle - oy(i)*std::copysign(sin_ocean_turning_angle, - * lat[i]) ); */ - double const tau_y = drag_atm * ay(i) - + c_prime * (oy(i) * cos_ocean_turning_angle + ox(i) * sin_ocean_turning_angle); - /* FIXME: Need latitude here. Then This becomes: - * + c_prime*( oy(i)*cos_ocean_turning_angle + ox(i)*std::copysign(sin_ocean_turning_angle, - * lat[i]) ); */ - - // We need to divide the gradient terms with the lumped mass matrix term - double const grad_x = tmpx(i) / pmap.lumpedcgmass(i); - double const grad_y = tmpy(i) / pmap.lumpedcgmass(i); - - vx(i) = alpha * uice + beta * vice - + dte_over_mass * (alpha * (grad_x + tau_x) + beta * (grad_y + tau_y)); - vx(i) *= rdenom; - - vy(i) = alpha * vice - beta * uice - + dte_over_mass * (alpha * (grad_y + tau_y) + beta * (grad_x + tau_x)); - vy(i) *= rdenom; - } - - DirichletZero(); - - avg_vx += vx / NT_meb; - avg_vy += vy / NT_meb; -} - -template -template -void CGParametricMomentum::MEBStep(const BBMParameters& params, const size_t NT_meb, - double dt_adv, const DGVector& H, const DGVector& A, DGVector& D) -{ - - double dt_mom = dt_adv / NT_meb; - - //! Compute Strain Rate - ProjectCGVelocityToDGStrain(); - - // TODO compute stress update with precomputed transformations - Nextsim::MEB::StressUpdateHighOrder( - params, smesh, S11, S12, S22, E11, E12, E22, H, A, D, dt_mom); - // Nextsim::BBM::StressUpdateHighOrder(params, smesh, S11, S12, S22, E11, - // E12, E22, H, A, D, dt_mom); - - // Nextsim::MEB::StressUpdateHighOrder(params, ptrans, smesh, S11, S12, S22, E11, E12, E22, H, - // A, D, dt_mom); - -#pragma omp parallel for - for (int i = 0; i < tmpx.rows(); ++i) - tmpx(i) = tmpy(i) = 0; - - // AddStressTensor(-1.0, tmpx, tmpy); - DivergenceOfStress(1.0, tmpx, tmpy); // Compute divergence of stress tensor - - double const cos_ocean_turning_angle = std::cos(params.oceanTurningAngle * M_PI / 180.); - double const sin_ocean_turning_angle = std::sin(params.oceanTurningAngle * M_PI / 180.); - - const double FOcean = params.COcean * params.rhoOcean; - const double FAtm = params.CAtm * params.rhoAtm; - -#pragma omp parallel for - for (int i = 0; i < vx.rows(); ++i) { - - double absatm = sqrt(ax(i) * ax(i) + ay(i) * ay(i)); - - double corsurf_x = vx(i) - ox(i); - double corsurf_y = vy(i) - oy(i); - double absocn = sqrt(SQR(corsurf_x) + SQR(corsurf_y)); - - vx(i) = (1.0 - / (params.rhoIce * cg_H(i) / dt_mom // implicit parts - + cg_A(i) * FOcean * absocn) // implicit parts - * (params.rhoIce * cg_H(i) / dt_mom * vx(i) - + cg_A(i) - * (FAtm * absatm * ax(i) + // atm forcing - FOcean * absocn * ox(i)) // ocean forcing - + params.rhoIce * cg_H(i) * params.fc * corsurf_y)); // cor + surface - - vy(i) = (1.0 - / (params.rhoIce * cg_H(i) / dt_mom // implicit parts - + cg_A(i) * FOcean * absocn) // implicit parts - * (params.rhoIce * cg_H(i) / dt_mom * vy(i) - + cg_A(i) - * (FAtm * absatm * ay(i) + // atm forcing - FOcean * absocn * oy(i)) // ocean forcing - + params.rhoIce * cg_H(i) * params.fc * corsurf_x)); // cor + surface - - vx(i) += (1.0 - / (params.rhoIce * cg_H(i) / dt_mom // implicit parts - + cg_A(i) * FOcean * absocn) // implicit parts - * tmpx(i)) - / pmap.lumpedcgmass(i); - ; - - vy(i) += (1.0 - / (params.rhoIce * cg_H(i) / dt_mom // implicit parts - + cg_A(i) * FOcean * absocn) // implicit parts - * tmpy(i)) - / pmap.lumpedcgmass(i); - ; - } - - DirichletZero(); - - avg_vx += vx / NT_meb; - avg_vy += vy / NT_meb; -} -// -------------------------------------------------- - -template class CGParametricMomentum<1>; -template class CGParametricMomentum<2>; - -template void CGParametricMomentum<1>::prepareIteration(const DGVector<1>& H, const DGVector<1>& A); -template void CGParametricMomentum<1>::prepareIteration(const DGVector<3>& H, const DGVector<3>& A); -template void CGParametricMomentum<1>::prepareIteration(const DGVector<6>& H, const DGVector<6>& A); -template void CGParametricMomentum<2>::prepareIteration(const DGVector<1>& H, const DGVector<1>& A); -template void CGParametricMomentum<2>::prepareIteration(const DGVector<3>& H, const DGVector<3>& A); -template void CGParametricMomentum<2>::prepareIteration(const DGVector<6>& H, const DGVector<6>& A); - -template void CGParametricMomentum<1>::prepareIteration( - const DGVector<1>& H, const DGVector<1>& A, const DGVector<1>& D); -template void CGParametricMomentum<1>::prepareIteration( - const DGVector<3>& H, const DGVector<3>& A, const DGVector<3>& D); -template void CGParametricMomentum<1>::prepareIteration( - const DGVector<6>& H, const DGVector<6>& A, const DGVector<6>& D); -template void CGParametricMomentum<2>::prepareIteration( - const DGVector<1>& H, const DGVector<1>& A, const DGVector<1>& D); -template void CGParametricMomentum<2>::prepareIteration( - const DGVector<3>& H, const DGVector<3>& A, const DGVector<3>& D); -template void CGParametricMomentum<2>::prepareIteration( - const DGVector<6>& H, const DGVector<6>& A, const DGVector<6>& D); - -// -------------------------------------------------- - -template void CGParametricMomentum<1>::mEVPStep(const VPParameters& params, size_t NT_evp, - double alpha, double beta, double dt_adv, const DGVector<1>& H, const DGVector<1>& A); -template void CGParametricMomentum<1>::mEVPStep(const VPParameters& params, size_t NT_evp, - double alpha, double beta, double dt_adv, const DGVector<3>& H, const DGVector<3>& A); -template void CGParametricMomentum<1>::mEVPStep(const VPParameters& params, size_t NT_evp, - double alpha, double beta, double dt_adv, const DGVector<6>& H, const DGVector<6>& A); - -template void CGParametricMomentum<2>::mEVPStep(const VPParameters& params, size_t NT_evp, - double alpha, double beta, double dt_adv, const DGVector<1>& H, const DGVector<1>& A); -template void CGParametricMomentum<2>::mEVPStep(const VPParameters& params, size_t NT_evp, - double alpha, double beta, double dt_adv, const DGVector<3>& H, const DGVector<3>& A); -template void CGParametricMomentum<2>::mEVPStep(const VPParameters& params, size_t NT_evp, - double alpha, double beta, double dt_adv, const DGVector<6>& H, const DGVector<6>& A); - -// -------------------------------------------------- - -template void CGParametricMomentum<1>::MEBStep(const BBMParameters& params, size_t NT_evp, - double dt_adv, const DGVector<1>& H, const DGVector<1>& A, DGVector<1>& D); -template void CGParametricMomentum<1>::MEBStep(const BBMParameters& params, size_t NT_evp, - double dt_adv, const DGVector<3>& H, const DGVector<3>& A, DGVector<3>& D); -template void CGParametricMomentum<1>::MEBStep(const BBMParameters& params, size_t NT_evp, - double dt_adv, const DGVector<6>& H, const DGVector<6>& A, DGVector<6>& D); - -template void CGParametricMomentum<2>::MEBStep(const BBMParameters& params, size_t NT_evp, - double dt_adv, const DGVector<1>& H, const DGVector<1>& A, DGVector<1>& D); -template void CGParametricMomentum<2>::MEBStep(const BBMParameters& params, size_t NT_evp, - double dt_adv, const DGVector<3>& H, const DGVector<3>& A, DGVector<3>& D); -template void CGParametricMomentum<2>::MEBStep(const BBMParameters& params, size_t NT_evp, - double dt_adv, const DGVector<6>& H, const DGVector<6>& A, DGVector<6>& D); - -// -------------------------------------------------- - -template void CGParametricMomentum<1>::BBMStep(const BBMParameters& params, size_t NT_evp, - double dt_adv, const DGVector<1>& H, const DGVector<1>& A, DGVector<1>& D); -template void CGParametricMomentum<1>::BBMStep(const BBMParameters& params, size_t NT_evp, - double dt_adv, const DGVector<3>& H, const DGVector<3>& A, DGVector<3>& D); -template void CGParametricMomentum<1>::BBMStep(const BBMParameters& params, size_t NT_evp, - double dt_adv, const DGVector<6>& H, const DGVector<6>& A, DGVector<6>& D); - -template void CGParametricMomentum<2>::BBMStep(const BBMParameters& params, size_t NT_evp, - double dt_adv, const DGVector<1>& H, const DGVector<1>& A, DGVector<1>& D); -template void CGParametricMomentum<2>::BBMStep(const BBMParameters& params, size_t NT_evp, - double dt_adv, const DGVector<3>& H, const DGVector<3>& A, DGVector<3>& D); -template void CGParametricMomentum<2>::BBMStep(const BBMParameters& params, size_t NT_evp, - double dt_adv, const DGVector<6>& H, const DGVector<6>& A, DGVector<6>& D); - -} /* namespace Nextsim */ diff --git a/dynamics/src/include/DynamicsKernel.hpp b/dynamics/src/include/DynamicsKernel.hpp index 08f9ecc19..a48979c7c 100644 --- a/dynamics/src/include/DynamicsKernel.hpp +++ b/dynamics/src/include/DynamicsKernel.hpp @@ -52,11 +52,6 @@ template class DynamicsKernel { if (isSpherical) smesh->RotatePoleToGreenland(); smesh->landmaskFromModelArray(mask); - smesh->dirichletFromMask(); - // TODO: handle periodic and open edges - for (ParametricMesh::Edge edge : ParametricMesh::edges) { - smesh->dirichletFromEdge(edge); - } //! Initialize transport dgtransport = std::make_unique>(*smesh); diff --git a/dynamics/src/include/ParametricMesh.hpp b/dynamics/src/include/ParametricMesh.hpp index 25e566cf1..b100d3628 100644 --- a/dynamics/src/include/ParametricMesh.hpp +++ b/dynamics/src/include/ParametricMesh.hpp @@ -55,19 +55,6 @@ class ParametricMesh { Eigen::Matrix vertices; // stores the - /*! - * Stores the Dirichlet boundary information - * - * Dirichlet-Information is stored in 4 vectors - * dirichlet[0]: List of elements that have a Dirichlet-Boundary on the bottom - * dirichlet[1]: List of elements that have a Dirichlet-Boundary on the right - * dirichlet[2]: List of elements that have a Dirichlet-Boundary on the top - * dirichlet[3]: List of elements that have a Dirichlet-Boundary on the left - * - * Each of the vectors is processed in parallel - */ - std::array, 4> dirichlet; - /*! * Periodic: * [0]=0 for X-edge (bottom/top) and [0]=1 for Y-edge (left/right) @@ -102,8 +89,6 @@ class ParametricMesh { ny = -1; nnodes = -1; nelements = -1; - for (size_t i = 0; i < 4; ++i) - dirichlet[i].clear(); periodic.clear(); landmask.clear(); } @@ -393,29 +378,6 @@ class ParametricMesh { */ void landmaskFromModelArray(const ModelArray& mask); - /*! - * Add to the dirichlet arrays according to the stored landmask. - */ - void dirichletFromMask(); - - /*! - * Add to the dirichlet arrays due to the domain edges according to an edge index. - * - * @param edge index of the edge to add closed boundary conditions to. - */ - void dirichletFromEdge(Edge edge); - - /*! - * Sort all the dirichlet arrays, so the element indices are ordered. - */ - void sortDirichlet(); - - /*! - * Sort the dirichlet array of one particular edge. - * - * @param edge the edge to be sorted. - */ - void sortDirichlet(Edge edge); // Global access functions diff --git a/dynamics/src/include/cgParametricMomentum.hpp b/dynamics/src/include/cgParametricMomentum.hpp deleted file mode 100644 index 6284f9a18..000000000 --- a/dynamics/src/include/cgParametricMomentum.hpp +++ /dev/null @@ -1,285 +0,0 @@ -/*! - * @file cgParametricMomentum.hpp - * @date 06 Dec 2024 - * @author Thomas Richter - */ - -#ifndef __CGPARAMETRICMOMENTUM_HPP -#define __CGPARAMETRICMOMENTUM_HPP - -#include "BBMParameters.hpp" -#include "ParametricMap.hpp" -#include "ParametricTools.hpp" -#include "VPParameters.hpp" -#include "cgVector.hpp" -#include "codeGenerationCGToDG.hpp" -#include "codeGenerationCGinGauss.hpp" -#include "codeGenerationDGinGauss.hpp" -#include "dgVector.hpp" -#include "dgVisu.hpp" - -namespace Nextsim { - -template class CGParametricMomentum { -private: - const ParametricMesh& smesh; //!< const-reference to the mesh - - /*! - * Stores precomputed values for efficient numerics on transformed mesh - * accelerates numerics but substantial memory effort! - */ - static constexpr int precompute_matrices = 1; - -public: - //! CGParametricMomentum should not be used but removed. ParametricMomentumMap needs DGadvection - //! but this is not provided here. - ParametricMomentumMap pmap; - - //! vectors storing the velocity (node-wise) - CGVector vx, vy; - - //! vectors storing the average velocity in sub-iterations (node-wise) - CGVector avg_vx, avg_vy; - - //! vectors storing ocean and atm velocity (node-wise) - CGVector ox, oy, ax, ay; - - ////// Internal variables - - //! temporary vectors. Maybe we can remove them? - CGVector tmpx, tmpy; - - //! old velocities. They are required during temporary vectors. Maybe we can remove them? - CGVector vx_mevp, vy_mevp; - - //! Vector to store the CG-Version of ice concentration and ice height - CGVector cg_A, cg_H, cg_D; - - //! Vectors storing strain and strss - DGVector E11, E12, E22; - DGVector S11, S12, S22; - -public: - CGParametricMomentum(const ParametricMesh& sm) - : smesh(sm) - , pmap(sm) - { - std::cerr << "CGParametricMomentum can't be used any more. It is replaced by the various " - "dynamics kernels. ParametricMomentumMap requires the DGadvection degree as " - "template but this is not available in CGParametricMomentum. Please use the " - "Kernel-infrastructure instead." - << std::endl; - - if (!(smesh.nelements > 0)) { - std::cerr << "CGParametricMomentum: The mesh has to be initialized first!" << std::endl; - abort(); - } - - // just simple sanity checks - assert(smesh.nelements == smesh.nx * smesh.ny); - assert(smesh.nnodes == (smesh.nx + 1) * (smesh.ny + 1)); - - // resize the vectors - vx.resize_by_mesh(smesh); - vy.resize_by_mesh(smesh); - vx.setZero(); - vy.setZero(); - - avg_vx.resize_by_mesh(smesh); - avg_vy.resize_by_mesh(smesh); - avg_vx.setZero(); - avg_vy.setZero(); - - cg_A.resize_by_mesh(smesh); - cg_H.resize_by_mesh(smesh); - cg_D.resize_by_mesh(smesh); - - ax.resize_by_mesh(smesh); - ay.resize_by_mesh(smesh); - ox.resize_by_mesh(smesh); - oy.resize_by_mesh(smesh); - - tmpx.resize_by_mesh(smesh); - tmpy.resize_by_mesh(smesh); - - E11.resize_by_mesh(smesh); - E12.resize_by_mesh(smesh); - E22.resize_by_mesh(smesh); - S11.resize_by_mesh(smesh); - S12.resize_by_mesh(smesh); - S22.resize_by_mesh(smesh); - - /*! - * initialize the lumped mass - * At periodic boundaries, the values must be added from both sides - */ - pmap.InitializeLumpedCGMassMatrix(); - - /*! - * Compute matrices for performing mEVP / BBM / MEB etc. stress updates - */ - pmap.InitializeDivSMatrices(); - } - - // Access to members - const CGVector& GetVx() const { return vx; } - const CGVector& GetVy() const { return vy; } - CGVector& GetVx() { return vx; } - CGVector& GetVy() { return vy; } - - const CGVector& GetAvgSubiterVx() const { return avg_vx; } - const CGVector& GetAvgSubiterVy() const { return avg_vy; } - CGVector& GetAvgSubiterVx() { return avg_vx; } - CGVector& GetAvgSubiterVy() { return avg_vy; } - - const CGVector& GetOceanx() const { return ox; } - const CGVector& GetOceany() const { return oy; } - CGVector& GetOceanx() { return ox; } - CGVector& GetOceany() { return oy; } - const CGVector& GetAtmx() const { return ax; } - const CGVector& GetAtmy() const { return ay; } - CGVector& GetAtmx() { return ax; } - CGVector& GetAtmy() { return ay; } - - const DGVector& GetE11() const { return E11; } - const DGVector& GetE12() const { return E12; } - const DGVector& GetE22() const { return E22; } - - const DGVector& GetS11() const { return S11; } - const DGVector& GetS12() const { return S12; } - const DGVector& GetS22() const { return S22; } - DGVector& GetS11() { return S11; } - DGVector& GetS12() { return S12; } - DGVector& GetS22() { return S22; } - - const CGVector& GetcgH() const { return cg_H; } - const CGVector& GetcgA() const { return cg_A; } - const CGVector& GetcgD() const { return cg_D; } - CGVector& GetcgH() { return cg_H; } - CGVector& GetcgA() { return cg_A; } - CGVector& GetcgD() { return cg_D; } - - // High level Functions - - /*! - * prepare the subcucling iteration: - * - store old velocity - * - interpoalte ice height & concentration ( & damage) to cg - */ - template void prepareIteration(const DGVector& H, const DGVector& A); - template - void prepareIteration(const DGVector& H, const DGVector& A, const DGVector& D); - - //! performs one complete mEVP cycle with NT_evp subiterations - template - void mEVPStep(const VPParameters& vpparameters, size_t NT_evp, double alpha, double beta, - double dt_adv, const DGVector& H, const DGVector& A); - - //! performs one complete MEB timestep with NT_meb subiterations - template - void MEBStep(const BBMParameters& vpparameters, size_t NT_meb, double dt_adv, - const DGVector& H, const DGVector& A, DGVector& D); - - //! performs one complete BBMStep timestep with NT_meb subiterations - template - void BBMStep(const BBMParameters& vpparameters, size_t NT_meb, double dt_adv, - const DGVector& H, const DGVector& A, DGVector& D); - - /*! - * The following functions take care of the interpolation and projection - * between CG and DG functions - */ - //! Projects the symmetric gradient of the CG velocity into the DG space - void ProjectCGVelocityToDGStrain(); - - /*! - * Evaluates (S, nabla phi) and writes it in the tx/ty - Vector - */ - void DivergenceOfStress(const double scale, CGVector& tx, CGVector& ty) const; - - void AddStressTensorCell(const double scale, const size_t c, const size_t cx, const size_t cy, - CGVector& tx, CGVector& ty) const; - - //! Sets the velocity vector to zero along the boundary - void DirichletZero() - { - DirichletZero(vx); - DirichletZero(vy); - } - void DirichletZero(CGVector& v) const; - - /*! - * AddPeriodic is to be called, after (sigma, Nabla Phi) is computed - * On periodic boundaries, the contributions from both sides must be added - */ - void AddPeriodic(CGVector& v); - /*! - * AveragePeriodic replaces the values on both sides by - * the average of them - */ - void AveragePeriodic(CGVector& v); - void CheckPeriodicity(CGVector& v); -}; - -template -void CGParametricMomentum::AddStressTensorCell(const double scale, const size_t eid, - const size_t cx, const size_t cy, CGVector& tmpx, CGVector& tmpy) const -{ - // pick the number of Gauss points according to the degree - -#define NGP (CG == 1 ? 2 : 3) - - Eigen::Vector tx = scale - * (pmap.divS1[eid] * S11.row(eid).transpose() + pmap.divS2[eid] * S12.row(eid).transpose()); - Eigen::Vector ty = scale - * (pmap.divS1[eid] * S12.row(eid).transpose() + pmap.divS2[eid] * S22.row(eid).transpose()); - - if (smesh.CoordinateSystem - == SPHERICAL) // In spherical coordinates there is the additional 'derivative term' arising - // from the derivative of the units - { - tx += scale * pmap.divM[eid] * S12.row(eid).transpose(); - ty -= scale * pmap.divM[eid] * S11.row(eid).transpose(); - } - - const size_t CGROW = CG * smesh.nx + 1; - const size_t cg_i = CG * CGROW * cy + CG * cx; //!< lower left CG-index in element (cx,cy) - - if (CG == 1) { - tmpx(cg_i + 0) += -tx(0); - tmpx(cg_i + 1) += -tx(1); - tmpx(cg_i + 0 + CGROW) += -tx(2); - tmpx(cg_i + 1 + CGROW) += -tx(3); - - tmpy(cg_i + 0) += -ty(0); - tmpy(cg_i + 1) += -ty(1); - tmpy(cg_i + 0 + CGROW) += -ty(2); - tmpy(cg_i + 1 + CGROW) += -ty(3); - } else if (CG == 2) { - tmpx(cg_i + 0) += -tx(0); - tmpx(cg_i + 1) += -tx(1); - tmpx(cg_i + 2) += -tx(2); - tmpx(cg_i + 0 + CGROW) += -tx(3); - tmpx(cg_i + 1 + CGROW) += -tx(4); - tmpx(cg_i + 2 + CGROW) += -tx(5); - tmpx(cg_i + 0 + CGROW * 2) += -tx(6); - tmpx(cg_i + 1 + CGROW * 2) += -tx(7); - tmpx(cg_i + 2 + CGROW * 2) += -tx(8); - - tmpy(cg_i + 0) += -ty(0); - tmpy(cg_i + 1) += -ty(1); - tmpy(cg_i + 2) += -ty(2); - tmpy(cg_i + 0 + CGROW) += -ty(3); - tmpy(cg_i + 1 + CGROW) += -ty(4); - tmpy(cg_i + 2 + CGROW) += -ty(5); - tmpy(cg_i + 0 + CGROW * 2) += -ty(6); - tmpy(cg_i + 1 + CGROW * 2) += -ty(7); - tmpy(cg_i + 2 + CGROW * 2) += -ty(8); - } else - abort(); -#undef NGP -} - -} /* namespace Nextsim */ - -#endif /* __CGMOMENTUM_HPP */ diff --git a/dynamics/test/AdvectionPeriodicBC_test.cpp b/dynamics/test/AdvectionPeriodicBC_test.cpp index dd56df878..45bdc28e8 100644 --- a/dynamics/test/AdvectionPeriodicBC_test.cpp +++ b/dynamics/test/AdvectionPeriodicBC_test.cpp @@ -45,7 +45,7 @@ const size_t NT0 = 50; size_t NT = NT0; } -bool WRITE_VTK = true; //!< set to true for vtk output +bool WRITE_VTK = false; //!< set to true for vtk output double TOL = 1.e-7; //!< tolerance for checking test results @@ -232,25 +232,22 @@ void create_rectanglemesh(Nextsim::ParametricMesh& smesh) // landmask: set all to ice smesh.landmask.resize(smesh.nelements, 1); - //// Boundary. Dirichlet on bottom / top, periodic left/right - smesh.dirichlet[0].resize(smesh.nx); // bottom boundary + //// Boundary. Dirichlet on bottom (by landmask) / top, periodic left/right, for (size_t i = 0; i < smesh.nx; ++i) - smesh.dirichlet[0][i] = i; - - smesh.dirichlet[2].resize(smesh.nx); // top boundary + smesh.landmask[i] = 0; // bot for (size_t i = 0; i < smesh.nx; ++i) - smesh.dirichlet[2][i] = smesh.nx * (smesh.ny - 1) + i; + smesh.landmask[i + smesh.nx * (smesh.ny-1) ] = 0; smesh.periodic.resize(1); - smesh.periodic[0].resize(smesh.ny); - for (size_t i = 0; i < smesh.ny; ++i) { + smesh.periodic[0].resize(smesh.ny-2); // leave out top and bottom element + for (size_t i = 0; i < smesh.ny-2; ++i) { smesh.periodic[0][i][0] = 1; // periodic on left/right boundary - smesh.periodic[0][i][1] = (i + 1) * smesh.nx + smesh.periodic[0][i][1] = (i+1 + 1) * smesh.nx - 1; // element on the left of the edge (so the most right one in the row) smesh.periodic[0][i][2] - = i * smesh.nx; // element right of the edge (so the very left one) as... + = (i+1) * smesh.nx; // element right of the edge (so the very left one) as... smesh.periodic[0][i][3] - = i * (smesh.nx + 1) + i; // the edge is the most right edge in the row + = (i+2) * (smesh.nx+1) - 1; // the edge is the most right edge in the row } } @@ -278,17 +275,25 @@ template void run(const std::array, 3>& exact) TEST_SUITE_BEGIN("Advection Periodic Boundary Conditions"); TEST_CASE("Advection Periodic Boundary Conditions") { + // std::array, 3> exact = // Exact values taken 26/06/2024 + // { std::array( + // { 1.0338503986019776e+00, 1.1451366598186583e+00, 1.0681593193338459e+00, + // 9.5252231195653603e-01, 8.1458581892610615e-01, 6.8950068528265651e-01 }), + // std::array( + // { 1.0949680727313791e+00, 8.3851748134278226e-01, 5.8501741891364523e-01, + // 4.0006092999256893e-01, 2.8493698219799068e-01, 2.0801424342840474e-01 }), + // std::array( + // { 7.1704413076190610e-01, 4.6135258391583606e-01, 3.2770311091970727e-01, + // 2.3424633076579465e-01, 1.7038068017215552e-01, 1.2486932724366147e-01 }) }; + + // New values 07/04/25 due to new way of handling boundaries. We now include a land-layer + // at the top and bottom. The slightly changed domain leads to changed values std::array, 3> exact = // Exact values taken 26/06/2024 - { std::array( - { 1.0338503986019776e+00, 1.1451366598186583e+00, 1.0681593193338459e+00, - 9.5252231195653603e-01, 8.1458581892610615e-01, 6.8950068528265651e-01 }), - std::array( - { 1.0949680727313791e+00, 8.3851748134278226e-01, 5.8501741891364523e-01, - 4.0006092999256893e-01, 2.8493698219799068e-01, 2.0801424342840474e-01 }), - std::array( - { 7.1704413076190610e-01, 4.6135258391583606e-01, 3.2770311091970727e-01, - 2.3424633076579465e-01, 1.7038068017215552e-01, 1.2486932724366147e-01 }) }; - + { std::array({1.0338503986019776e+00,1.1451366598186583e+00,1.0681593193338459e+00,9.5252231195653614e-01,8.1458581892610615e-01,6.8950068528265651e-01}), + std::array({8.7678853597311990e-01,8.2389349440313742e-01,5.8562737114004948e-01,4.0006113872450388e-01,2.8493698219799091e-01,2.0801424342840488e-01}), + std::array({6.0231338733029061e-01,4.6967937778236674e-01,3.2769774426726950e-01,2.3427321065537998e-01,1.7038111756667601e-01,1.2487338791917807e-01}) + }; + std::cout << "DG\tNT\tNX\tmass\t\terror\t\texact" << std::endl; std::cout << std::setprecision(4) << std::scientific; run<1>(exact); diff --git a/dynamics/test/Advection_test.cpp b/dynamics/test/Advection_test.cpp index 4eebab45f..bf096fed9 100644 --- a/dynamics/test/Advection_test.cpp +++ b/dynamics/test/Advection_test.cpp @@ -218,29 +218,19 @@ void create_rectanglemesh(Nextsim::ParametricMesh& smesh, const double Lx, const } //// Boundary - // landmask: set all to ice + // landmask: set all to ice but outer layer along the domain smesh.landmask.resize(smesh.nelements, 1); - - smesh.dirichlet[0].resize(smesh.nx); // bottom boundary for (size_t i = 0; i < smesh.nx; ++i) - smesh.dirichlet[0][i] = i; - - smesh.dirichlet[1].resize(smesh.ny); // right boundary - for (size_t i = 0; i < smesh.ny; ++i) - smesh.dirichlet[1][i] = i * Nx + Nx - 1; - - smesh.dirichlet[2].resize(smesh.nx); // top boundary + smesh.landmask[i] = 0; for (size_t i = 0; i < smesh.nx; ++i) - smesh.dirichlet[2][i] = Nx * (Ny - 1) + i; - - smesh.dirichlet[3].resize(smesh.ny); // left boundary + smesh.landmask[smesh.nx * (smesh.ny-1) + i] = 0; for (size_t i = 0; i < smesh.ny; ++i) - smesh.dirichlet[3][i] = i * Nx; - - // no other boundaries. + smesh.landmask[i * smesh.nx] = 0; + for (size_t i = 0; i < smesh.ny; ++i) + smesh.landmask[(i+1) * smesh.nx-1 ] = 0; } -template void run(double distort, const std::array, 3>& exact) +template void run(double distort, const std::array, 3>& exact) { Nextsim::ParametricMesh smesh(Nextsim::CARTESIAN); // 0 means no output @@ -271,17 +261,29 @@ TEST_CASE("Advection") // Exact values taken 25/06/2024 after some plausibility check of the results // and by checking the theoretical order of convergence to be expected - std::array, 3> exact - = { std::array( - { 5.1256149074257538e-02, 4.8288256703303903e-02, 4.3105635248809886e-02, - 3.5793986049422598e-02, 2.6937487824223016e-02, 1.8456775604583933e-02 }), - std::array( - { 2.9450967798560313e-02, 1.3427281824939470e-02, 5.8574889800512000e-03, - 2.2756813704797301e-03, 7.3564079504566610e-04, 1.9177169540867740e-04 }), - std::array( - { 9.9340386651904228e-03, 4.0274889287136816e-03, 1.2400186092664943e-03, - 2.8460130790583933e-04, 4.5459555769938981e-05, 4.6851425365137733e-06 }) }; - + // std::array, 3> exact + // = { std::array( + // { 5.1256149074257538e-02, 4.8288256703303903e-02, 4.3105635248809886e-02, + // 3.5793986049422598e-02, 2.6937487824223016e-02, 1.8456775604583933e-02 }), + // std::array( + // { 2.9450967798560313e-02, 1.3427281824939470e-02, 5.8574889800512000e-03, + // 2.2756813704797301e-03, 7.3564079504566610e-04, 1.9177169540867740e-04 }), + // std::array( + // { 9.9340386651904228e-03, 4.0274889287136816e-03, 1.2400186092664943e-03, + // 2.8460130790583933e-04, 4.5459555769938981e-05, 4.6851425365137733e-06 }) }; + + + // New values 07/04/2025 after reimplementation of boundary data. Dirichlet data is now + // realized using the landmaks info. Therefore, the physical domain shrinks by one layer all around + // and this slightly changes the exact values. + std::array, 3> exact + = { std::array( + { 5.0568601321622546e-02,4.8785010758681462e-02,4.3807523095968269e-02,3.6079575213776381e-02 }), + std::array( + { 2.8463796850208754e-02,1.4000140482658283e-02,5.8605850177331619e-03,2.2757018442945147e-03 }), + std::array( + { 1.0146940898974235e-02,4.0853715304369851e-03,1.2400975498406005e-03,2.8460132054714237e-04}) }; + std::cout << std::endl << "DG\tNT\tNX\tmass loss\terror\t\texact" << std::endl; run<1>(0.0, exact); run<3>(0.0, exact); @@ -289,16 +291,27 @@ TEST_CASE("Advection") } TEST_CASE("Distorted Mesh") { - std::array, 3> exact - = { std::array( - { 5.0958748236875594e-02, 4.8461087243594887e-02, 4.3918572621094686e-02, - 3.7038043078562323e-02, 2.8334464207979679e-02, 1.9666196904181983e-02 }), - std::array( - { 3.2033423984965226e-02, 1.5639052766007147e-02, 6.7926997203524983e-03, - 2.7369916393700151e-03, 9.2109964949992139e-04, 2.5128064874203957e-04 }), - std::array( - { 1.1830071142946669e-02, 4.9207680503709564e-03, 1.5831512504162868e-03, - 3.8869595113351363e-04, 6.7162895595772516e-05, 7.5853786764529606e-06 }) }; + // std::array, 3> exact + // = { std::array( + // { 5.0958748236875594e-02, 4.8461087243594887e-02, 4.3918572621094686e-02, + // 3.7038043078562323e-02, 2.8334464207979679e-02, 1.9666196904181983e-02 }), + // std::array( + // { 3.2033423984965226e-02, 1.5639052766007147e-02, 6.7926997203524983e-03, + // 2.7369916393700151e-03, 9.2109964949992139e-04, 2.5128064874203957e-04 }), + // std::array( + // { 1.1830071142946669e-02, 4.9207680503709564e-03, 1.5831512504162868e-03, + // 3.8869595113351363e-04, 6.7162895595772516e-05, 7.5853786764529606e-06 }) }; + + // New values 07/04/2025 after reimplementation of boundary data. Dirichlet data is now + // realized using the landmaks info. Therefore, the physical domain shrinks by one layer all around + // and this slightly changes the exact values. + std::array, 3> exact + = { std::array( + { 4.8195942137802428e-02,4.9345757529052611e-02,4.4903895856787640e-02,3.7467928845481245e-02 }), + std::array( + { 3.2422647655365219e-02,1.7062569520319999e-02,6.8213111669173628e-03,2.7372092446537699e-03 }), + std::array( + { 1.2603978369132178e-02,5.1434197544894689e-03,1.5845944976351275e-03,3.8869603127842235e-04}) }; std::cout << std::endl << "Distorted mesh" << std::endl; std::cout << "DG\tNT\tNX\tmass loss\terror\t\texact" << std::endl; diff --git a/dynamics/test/ParametricMesh_test.cpp b/dynamics/test/ParametricMesh_test.cpp index 28d1d243d..03c277bae 100644 --- a/dynamics/test/ParametricMesh_test.cpp +++ b/dynamics/test/ParametricMesh_test.cpp @@ -65,19 +65,6 @@ TEST_CASE("Test readmesh") REQUIRE(smesh.landmask[nx - 1]); REQUIRE(!smesh.landmask[(ny - 1) * nx]); - // Dirichlet conditions - // Element 5 comes first, and is closed on edges 0 & 3 - REQUIRE(smesh.dirichlet[0][0] == 5); - REQUIRE(smesh.dirichlet[3][0] == 5); - // Element 7 is the first with a closed edge 1 - REQUIRE(smesh.dirichlet[1][0] == 7); - // Element 48 is the first with a closed edge 2 - REQUIRE(smesh.dirichlet[2][0] == 48); - // sizes of the Dirichlet arrays - REQUIRE(smesh.dirichlet[0].size() == 369); - REQUIRE(smesh.dirichlet[1].size() == 384); - REQUIRE(smesh.dirichlet[0].size() == smesh.dirichlet[2].size()); - REQUIRE(smesh.dirichlet[1].size() == smesh.dirichlet[3].size()); // No periodic boundary conditions REQUIRE(smesh.periodic.size() == 0); @@ -119,21 +106,6 @@ TEST_CASE("Compare readmesh and landmask reading") REQUIRE(fromArrays.landmask[idx] == fromFile.landmask[idx]); } - // Dirichlet conditions - fromArrays.dirichletFromMask(); - for (ParametricMesh::Edge edge : ParametricMesh::edges) { - fromArrays.dirichletFromEdge(edge); - } - fromArrays.sortDirichlet(); - for (ParametricMesh::Edge edge : ParametricMesh::edges) { - REQUIRE(fromArrays.dirichlet[edge].size() == fromFile.dirichlet[edge].size()); - } - - for (ParametricMesh::Edge edge : ParametricMesh::edges) { - for (size_t idx = 0; idx < fromArrays.dirichlet[edge].size(); ++idx) { - REQUIRE(fromArrays.dirichlet[edge][idx] == fromFile.dirichlet[edge][idx]); - } - } // No periodic boundary conditions REQUIRE(fromArrays.periodic.size() == 0); } From 1a3c692315fcf660d4fa3acfd54bcb4e26c25d29 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Einar=20=C3=96rn=20=C3=93lason?= Date: Wed, 30 Apr 2025 09:01:21 +0200 Subject: [PATCH 08/33] Remove VTK output from VPCGDynamicsKernel It was there for debugging purposes only. --- dynamics/src/include/VPCGDynamicsKernel.hpp | 15 +-------------- 1 file changed, 1 insertion(+), 14 deletions(-) diff --git a/dynamics/src/include/VPCGDynamicsKernel.hpp b/dynamics/src/include/VPCGDynamicsKernel.hpp index 4aa5fbe4a..5b62d61e0 100644 --- a/dynamics/src/include/VPCGDynamicsKernel.hpp +++ b/dynamics/src/include/VPCGDynamicsKernel.hpp @@ -1,7 +1,7 @@ /*! * @file VPCGDynamicsKernel.hpp * - * @date 27 Mar 2025 + * @date 30 Apr 2025 * @author Tim Spain * @author Robert Jendersie */ @@ -65,19 +65,6 @@ template class VPCGDynamicsKernel : public CGDynamicsKernel::stepNumber%si==0) - { - int vtkn = DynamicsKernel::stepNumber/si; - VTK::write_dg("Polynya/hice", vtkn, hice, *smesh); - VTK::write_dg("Polynya/cice", vtkn, cice ,*smesh); - VTK::write_cg_velocity("Polynya/vel",vtkn, u,v,*smesh); - VTK::write_cg("Polynya/cgl",vtkn, pmap->cglandmask,*smesh); - } - - - // Let DynamicsKernel handle the advection step DynamicsKernel::advectionAndLimits(tst); From be756c3d2884add73d237c6987e013b6d22be6b3 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Einar=20=C3=96rn=20=C3=93lason?= Date: Wed, 30 Apr 2025 09:12:46 +0200 Subject: [PATCH 09/33] Add u and v wind direction to ConfiguredAtmosphere This makes comparison with the published works easier, as they use something similar to ConfiguredAtmosphere, not FluxConfiguredAtmosphere. --- .../ConfiguredAtmosphere.cpp | 27 +++++++++++++------ .../include/ConfiguredAtmosphere.hpp | 8 +++--- 2 files changed, 24 insertions(+), 11 deletions(-) diff --git a/physics/src/modules/AtmosphereBoundaryModule/ConfiguredAtmosphere.cpp b/physics/src/modules/AtmosphereBoundaryModule/ConfiguredAtmosphere.cpp index d124ce397..5ddb2f0af 100644 --- a/physics/src/modules/AtmosphereBoundaryModule/ConfiguredAtmosphere.cpp +++ b/physics/src/modules/AtmosphereBoundaryModule/ConfiguredAtmosphere.cpp @@ -1,7 +1,7 @@ /*! * @file ConfiguredAtmosphere.cpp * - * @date 20 Nov 2024 + * @date 30 Apr 2025 * @author Tim Spain */ @@ -10,6 +10,8 @@ #include "include/Finalizer.hpp" #include "include/NextsimModule.hpp" +#include + namespace Nextsim { double ConfiguredAtmosphere::tair0 = 0; // Freezing @@ -19,7 +21,8 @@ double ConfiguredAtmosphere::sw0 = 0; // Night double ConfiguredAtmosphere::lw0 = 315.637; // Stefan-Boltzmann at 0˚C double ConfiguredAtmosphere::snowfall0 = 0; // No snow double ConfiguredAtmosphere::rain0 = 0; // No rain -double ConfiguredAtmosphere::windspeed0 = 0; // Still +double ConfiguredAtmosphere::uWind0 = 0; // Still +double ConfiguredAtmosphere::vWind0 = 0; // Still static const std::string pfx = "ConfiguredAtmosphere"; static const std::string tKey = pfx + ".t_air"; @@ -29,7 +32,8 @@ static const std::string swKey = pfx + ".sw_in"; static const std::string lwKey = pfx + ".lw_in"; static const std::string snowKey = pfx + ".snow"; static const std::string rainKey = pfx + ".rainfall"; -static const std::string windKey = pfx + ".wind_speed"; +static const std::string uWindKey = pfx + ".wind_u"; +static const std::string vWindKey = pfx + ".wind_v"; const static std::map keyMap = { { ConfiguredAtmosphere::TAIR_KEY, tKey }, @@ -39,7 +43,8 @@ const static std::map keyMap = { { ConfiguredAtmosphere::LW_KEY, lwKey }, { ConfiguredAtmosphere::SNOW_KEY, snowKey }, { ConfiguredAtmosphere::RAIN_KEY, rainKey }, - { ConfiguredAtmosphere::WIND_KEY, windKey }, + { ConfiguredAtmosphere::UWIND_KEY, uWindKey }, + { ConfiguredAtmosphere::VWIND_KEY, vWindKey }, }; ConfiguredAtmosphere::ConfiguredAtmosphere() @@ -70,8 +75,10 @@ ConfigurationHelp::HelpMap& ConfiguredAtmosphere::getHelpRecursive(HelpMap& map, "", "Snowfall mass flux (kg s⁻¹ m⁻²)." }, { rainKey, ConfigType::NUMERIC, { "0", "2.998e8" }, ConfigurationHelp::toString(rain0), "", "Rainfall mass flux (kg s⁻¹ m⁻²)." }, - { windKey, ConfigType::NUMERIC, { "0", "2.998e8" }, ConfigurationHelp::toString(windspeed0), - "", "Windspeed (m s⁻¹)." }, + { uWindKey, ConfigType::NUMERIC, { "0", "2.998e8" }, ConfigurationHelp::toString(uWind0), + "m/s", "Component of wind in the x (eastward) direction (m s⁻¹)." }, + { vWindKey, ConfigType::NUMERIC, { "0", "2.998e8" }, ConfigurationHelp::toString(vWind0), + "m/s", "Component of wind in the y (northward) direction (m s⁻¹)." }, }; Module::getHelpRecursive(map, getAll); @@ -87,7 +94,8 @@ void ConfiguredAtmosphere::configure() lw0 = Configured::getConfiguration(keyMap.at(LW_KEY), lw0); snowfall0 = Configured::getConfiguration(keyMap.at(SNOW_KEY), snowfall0); rain0 = Configured::getConfiguration(keyMap.at(RAIN_KEY), rain0); - windspeed0 = Configured::getConfiguration(keyMap.at(WIND_KEY), windspeed0); + uWind0 = Configured::getConfiguration(keyMap.at(UWIND_KEY), uWind0); + vWind0 = Configured::getConfiguration(keyMap.at(VWIND_KEY), vWind0); Finalizer::registerUnique(Module::finalize); fluxImpl = &Module::getImplementation(); @@ -112,7 +120,10 @@ void ConfiguredAtmosphere::setData(const ModelState::DataMap& dm) lw_in = lw0; snow = snowfall0; rain = rain0; - wind = windspeed0; + wind = std::hypot(uWind0, vWind0); + + uwind = uWind0; + vwind = vWind0; fluxImpl->setData(dm); } diff --git a/physics/src/modules/AtmosphereBoundaryModule/include/ConfiguredAtmosphere.hpp b/physics/src/modules/AtmosphereBoundaryModule/include/ConfiguredAtmosphere.hpp index e8b93939d..10cd144f5 100644 --- a/physics/src/modules/AtmosphereBoundaryModule/include/ConfiguredAtmosphere.hpp +++ b/physics/src/modules/AtmosphereBoundaryModule/include/ConfiguredAtmosphere.hpp @@ -1,7 +1,7 @@ /*! * @file ConfiguredAtmosphere.hpp * - * @date Aug 31, 2022 + * @date 30 Apr 2025 * @author Tim Spain */ @@ -28,7 +28,8 @@ class ConfiguredAtmosphere : public IAtmosphereBoundary, public Configured Date: Wed, 30 Apr 2025 09:14:02 +0200 Subject: [PATCH 10/33] Initialise damage to one And don't warn about zeros any more. The default rheology is still mEVP (as this is directly comparable to the published work). --- run/make_init_polynya.py | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/run/make_init_polynya.py b/run/make_init_polynya.py index 852769243..7704a01df 100644 --- a/run/make_init_polynya.py +++ b/run/make_init_polynya.py @@ -18,7 +18,7 @@ fname = f"init_polynya.nc" # The model expects everything in metres -initializer = initMaker(fname, nfirst, nsecond, nLayers, res*1e3) +initializer = initMaker(fname, nfirst, nsecond, nLayers, res*1e3, checkZeros=False) # Ice everywhere and all boundaries closed, except the x = 100 km end initializer.mask[:, :] = 1. @@ -33,6 +33,9 @@ # Uniform thickness of 20 cm initializer.hice[:, :] = 0.2 +# Undamaged ice +initializer.damage[:, :] = 1. + # Ice and ocean temperature and salinity at the freezing point ice_salinity = 5 # should match Ice::s in constants.hpp mu: float = -0.055 # should match Water::mu in constants.hpp From 21cf14e15c68de4ddbf3bf048349f1280336d6fc Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Einar=20=C3=96rn=20=C3=93lason?= Date: Wed, 30 Apr 2025 11:35:52 +0200 Subject: [PATCH 11/33] Fix open-water flux in FluxConfiguredAtmosphere The open water heat flux is reset by the thermodynamics, that the (slab) ocean doesn't super cool. This commit sets it back to the initial value on every FluxConfiguredAtmosphere::update call. --- .../FluxConfiguredAtmosphere.cpp | 9 ++++++++- .../include/FluxConfiguredAtmosphere.hpp | 6 +++--- 2 files changed, 11 insertions(+), 4 deletions(-) diff --git a/physics/src/modules/AtmosphereBoundaryModule/FluxConfiguredAtmosphere.cpp b/physics/src/modules/AtmosphereBoundaryModule/FluxConfiguredAtmosphere.cpp index 80b6eacb2..51b3f9fff 100644 --- a/physics/src/modules/AtmosphereBoundaryModule/FluxConfiguredAtmosphere.cpp +++ b/physics/src/modules/AtmosphereBoundaryModule/FluxConfiguredAtmosphere.cpp @@ -1,7 +1,7 @@ /*! * @file FluxConfiguredAtmosphere.cpp * - * @date 20 Nov 2024 + * @date 30 Apr 2025 * @author Tim Spain */ @@ -93,4 +93,11 @@ void FluxConfiguredAtmosphere::setData(const ModelState::DataMap& dm) vwind = v0; } +void FluxConfiguredAtmosphere::update(const TimestepTime& tst) +{ + /* The open water heat flux is reset by the thermodynamics, so that the (slab) ocean doesn't + * super cool. Therefore, we have to reset it here on every update */ + qow = qow0; +} + } /* namespace Nextsim */ diff --git a/physics/src/modules/AtmosphereBoundaryModule/include/FluxConfiguredAtmosphere.hpp b/physics/src/modules/AtmosphereBoundaryModule/include/FluxConfiguredAtmosphere.hpp index 084e76ebc..7edc7280f 100644 --- a/physics/src/modules/AtmosphereBoundaryModule/include/FluxConfiguredAtmosphere.hpp +++ b/physics/src/modules/AtmosphereBoundaryModule/include/FluxConfiguredAtmosphere.hpp @@ -1,7 +1,7 @@ /*! * @file FluxConfiguredAtmosphere.hpp * - * @date Sep 29, 2022 + * @date 30 Apr 2025 * @author Tim Spain */ @@ -41,8 +41,8 @@ class FluxConfiguredAtmosphere : public IAtmosphereBoundary, void configure() override; protected: - //! Performs the implementation specific updates. Does nothing. - void update(const TimestepTime&) override { } + //! Performs the implementation specific updates. + void update(const TimestepTime&) override; private: static double qia0; From 4c6a41aeaed89126e65974491a82468143a1627f Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Einar=20=C3=96rn=20=C3=93lason?= Date: Wed, 30 Apr 2025 11:37:38 +0200 Subject: [PATCH 12/33] Update config_polynya.cfg Run for the full five days (as previous publications did), use Winton thermodynamics (again) and tweak the ice-atmosphere and open-water heat fluxes. --- run/config_polynya.cfg | 9 ++++----- 1 file changed, 4 insertions(+), 5 deletions(-) diff --git a/run/config_polynya.cfg b/run/config_polynya.cfg index 61fcdb6d6..b05b4e1c3 100644 --- a/run/config_polynya.cfg +++ b/run/config_polynya.cfg @@ -1,7 +1,7 @@ [model] init_file = init_polynya.nc start = 2023-01-01T00:00:00Z -stop = 2023-01-03T00:00:00Z +stop = 2023-01-05T00:00:00Z time_step = P0-0T00:01:00 missing_value = 1e20 @@ -10,8 +10,7 @@ DiagnosticOutputModule = Nextsim::ConfigOutput DynamicsModule = Nextsim::MEVPDynamics AtmosphereBoundaryModule = Nextsim::FluxConfiguredAtmosphere OceanBoundaryModule = Nextsim::FluxConfiguredOcean -#IceThermodynamicsModule = Nextsim::ThermoWinton -IceThermodynamicsModule = Nextsim::DummyIceThermodynamics +IceThermodynamicsModule = Nextsim::ThermoWinton [ConfigOutput] period = P0-0T04:00:00 @@ -27,8 +26,8 @@ current_u = 0. current_v = 0. [FluxConfiguredAtmosphere] -Q_ia = 2 -Q_ow = 200 +Q_ia = 30 +Q_ow = 300 dQia_dT = 0 wind_u = 16 wind_v = 12 From 768c4d4cd4a15d85f2940e68530993a02f030c82 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Einar=20=C3=96rn=20=C3=93lason?= Date: Wed, 30 Apr 2025 11:39:21 +0200 Subject: [PATCH 13/33] Minor tweaks and fixes Added a lot of const-s, changed one equation to be more readable (for me at least), modified the limits to include fixing negative concentration and thickness values, and removed redundant qualifiers. --- physics/src/IceGrowth.cpp | 18 +++++++-------- .../FluxConfiguredAtmosphere.cpp | 18 +++++++-------- .../FiniteElementFluxes.cpp | 22 +++++++++---------- 3 files changed, 29 insertions(+), 29 deletions(-) diff --git a/physics/src/IceGrowth.cpp b/physics/src/IceGrowth.cpp index d1417fda8..263c56d46 100644 --- a/physics/src/IceGrowth.cpp +++ b/physics/src/IceGrowth.cpp @@ -1,7 +1,7 @@ /*! * @file IceGrowth.cpp * - * @date 20 Nov 2024 + * @date 30 Apr 2025 * @author Tim Spain * @author Einar Ólason */ @@ -203,22 +203,22 @@ void IceGrowth::newIceFormation(size_t i, const TimestepTime& tst) { // Flux cooling the ocean from open water // TODO Add assimilation fluxes here - double coolingFlux = qow[i]; + const double coolingFlux = qow[i]; // Temperature change of the mixed layer during this timestep - double deltaTml = -coolingFlux / mixedLayerBulkHeatCapacity[i] * tst.step; + const double deltaTml = -coolingFlux / mixedLayerBulkHeatCapacity[i] * tst.step; // Initial temperature - double t0 = sst[i]; + const double t0 = sst[i]; // Freezing point temperature - double tf0 = tf[i]; + const double tf0 = tf[i]; // Final temperature - double t1 = t0 + deltaTml; + const double t1 = t0 + deltaTml; // deal with cooling below the freezing point if (t1 < tf0) { // Heat lost cooling the mixed layer to freezing point - double sensibleFlux = (tf0 - t0) / deltaTml * coolingFlux; + const double sensibleFlux = -(tf0 - t0) * mixedLayerBulkHeatCapacity[i] / tst.step; // Any heat beyond that is latent heat forming new ice - double latentFlux = coolingFlux - sensibleFlux; + const double latentFlux = coolingFlux - sensibleFlux; qow[i] = sensibleFlux; newice[i] = latentFlux * tst.step * (1 - cice[i]) / (Ice::Lf * Ice::rho); @@ -260,7 +260,7 @@ void IceGrowth::lateralIceSpread(size_t i, const TimestepTime& tstep) void IceGrowth::applyLimits(size_t i, const TimestepTime& tstep) { - if ((0. < cice[i] && cice[i] < IceMinima::c()) || (0. < hice[i] && hice[i] < IceMinima::h())) { + if (cice[i] < IceMinima::c() || hice[i] < IceMinima::h()) { qow[i] += cice[i] * Water::Lf * (hice[i] * Ice::rho + hsnow[i] * Ice::rhoSnow) / tstep.step; hice[i] = 0; cice[i] = 0; diff --git a/physics/src/modules/AtmosphereBoundaryModule/FluxConfiguredAtmosphere.cpp b/physics/src/modules/AtmosphereBoundaryModule/FluxConfiguredAtmosphere.cpp index 51b3f9fff..2365fdad4 100644 --- a/physics/src/modules/AtmosphereBoundaryModule/FluxConfiguredAtmosphere.cpp +++ b/physics/src/modules/AtmosphereBoundaryModule/FluxConfiguredAtmosphere.cpp @@ -68,15 +68,15 @@ ConfigurationHelp::HelpMap& FluxConfiguredAtmosphere::getHelpRecursive(HelpMap& } void FluxConfiguredAtmosphere::configure() { - qia0 = Configured::getConfiguration(keyMap.at(QIA_KEY), qia0); - dqia_dt0 = Configured::getConfiguration(keyMap.at(DQIA_DT_KEY), dqia_dt0); - qow0 = Configured::getConfiguration(keyMap.at(QOW_KEY), qow0); - subl0 = Configured::getConfiguration(keyMap.at(SUBL_KEY), subl0); - snowfall0 = Configured::getConfiguration(keyMap.at(SNOW_KEY), snowfall0); - rain0 = Configured::getConfiguration(keyMap.at(RAIN_KEY), rain0); - evap0 = Configured::getConfiguration(keyMap.at(EVAP_KEY), evap0); - u0 = Configured::getConfiguration(keyMap.at(WINDU_KEY), u0); - v0 = Configured::getConfiguration(keyMap.at(WINDV_KEY), v0); + qia0 = getConfiguration(keyMap.at(QIA_KEY), qia0); + dqia_dt0 = getConfiguration(keyMap.at(DQIA_DT_KEY), dqia_dt0); + qow0 = getConfiguration(keyMap.at(QOW_KEY), qow0); + subl0 = getConfiguration(keyMap.at(SUBL_KEY), subl0); + snowfall0 = getConfiguration(keyMap.at(SNOW_KEY), snowfall0); + rain0 = getConfiguration(keyMap.at(RAIN_KEY), rain0); + evap0 = getConfiguration(keyMap.at(EVAP_KEY), evap0); + u0 = getConfiguration(keyMap.at(WINDU_KEY), u0); + v0 = getConfiguration(keyMap.at(WINDV_KEY), v0); } void FluxConfiguredAtmosphere::setData(const ModelState::DataMap& dm) diff --git a/physics/src/modules/FluxCalculationModule/FiniteElementFluxes.cpp b/physics/src/modules/FluxCalculationModule/FiniteElementFluxes.cpp index 358d8d2e9..d8dfbf7ef 100644 --- a/physics/src/modules/FluxCalculationModule/FiniteElementFluxes.cpp +++ b/physics/src/modules/FluxCalculationModule/FiniteElementFluxes.cpp @@ -1,7 +1,7 @@ /*! * @file FiniteElementFluxes.cpp * - * @date 11 Feb 2025 + * @date 30 Apr 2025 * @author Tim Spain */ @@ -140,6 +140,8 @@ inline double FiniteElementFluxes::dragOcean_m(double windSpeed) void FiniteElementFluxes::calculateIce(size_t i, const TimestepTime& tst) { + // An alias for surface temperature to make things more readable + const double& Tsurf = tice.zIndexAndLayer(i, 0); // Mass flux ice subl[i] = dragIce_t * rho_air[i] * windSpeed[i] * (sh_ice[i] - sh_air[i]); @@ -147,28 +149,26 @@ void FiniteElementFluxes::calculateIce(size_t i, const TimestepTime& tst) // Heat flux ice-atmosphere // Latent heat from sublimation - Q_lh_ia[i] = subl[i] * latentHeatIce(tice.zIndexAndLayer(i, 0)); - double dmdot_dT = dragIce_t * rho_air[i] * windSpeed[i] * dshice_dT[i]; - double dQlh_dT = latentHeatIce(tice.zIndexAndLayer(i, 0)) * dmdot_dT; + Q_lh_ia[i] = subl[i] * latentHeatIce(Tsurf); + const double dmdot_dT = dragIce_t * rho_air[i] * windSpeed[i] * dshice_dT[i]; + const double dQlh_dT = latentHeatIce(Tsurf) * dmdot_dT; // Sensible heat flux - Q_sh_ia[i] = dragIce_t * rho_air[i] * cp_air[i] * windSpeed[i] - * (tice.zIndexAndLayer(i, 0) - t_air[i]); - double dQsh_dT = dragIce_t * rho_air[i] * cp_air[i] * windSpeed[i]; + Q_sh_ia[i] = dragIce_t * rho_air[i] * cp_air[i] * windSpeed[i] * (Tsurf - t_air[i]); + const double dQsh_dT = dragIce_t * rho_air[i] * cp_air[i] * windSpeed[i]; // Shortwave flux double albedoValue, i0; std::tie(albedoValue, i0) - = iIceAlbedoImpl->surfaceShortWaveBalance(tice.zIndexAndLayer(i, 0), h_snow_true[i], m_I0); + = iIceAlbedoImpl->surfaceShortWaveBalance(Tsurf, h_snow_true[i], m_I0); Q_sw_ia[i] = -sw_in[i] * (1. - albedoValue) * (1. - i0); const double extinction = 0.; // TODO: Replace with de Beer's law or a module penSW[i] = sw_in[i] * (1. - albedoValue) * i0 * (1. - extinction); Q_sw_base[i] = sw_in[i] * (1. - albedoValue) * i0 * extinction; // Longwave flux - Q_lw_ia[i] = stefanBoltzmannLaw(tice.zIndexAndLayer(i, 0)) - lw_in[i]; - double dQlw_dT - = 4 / kelvin(tice.zIndexAndLayer(i, 0)) * stefanBoltzmannLaw(tice.zIndexAndLayer(i, 0)); + Q_lw_ia[i] = stefanBoltzmannLaw(Tsurf) - lw_in[i]; + const double dQlw_dT = 4 / kelvin(Tsurf) * stefanBoltzmannLaw(Tsurf); // Total flux qia[i] = Q_lh_ia[i] + Q_sh_ia[i] + Q_sw_ia[i] + Q_lw_ia[i]; From e10ccbfc3564c832c3f431d8d9e67c66516e6e35 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Einar=20=C3=96rn=20=C3=93lason?= Date: Wed, 30 Apr 2025 11:48:56 +0200 Subject: [PATCH 14/33] Clang-format fixes The test failed, so I ran clang-format on the files. --- dynamics/src/CGDynamicsKernel.cpp | 14 ++--- dynamics/src/DGTransport.cpp | 78 ++++++++++++------------- dynamics/src/ParametricMesh.cpp | 14 ++--- dynamics/src/include/ParametricMesh.hpp | 3 +- 4 files changed, 52 insertions(+), 57 deletions(-) diff --git a/dynamics/src/CGDynamicsKernel.cpp b/dynamics/src/CGDynamicsKernel.cpp index 2ea5e9088..012890fe0 100644 --- a/dynamics/src/CGDynamicsKernel.cpp +++ b/dynamics/src/CGDynamicsKernel.cpp @@ -1,7 +1,7 @@ /*! * @file CGDynamicsKernel.cpp * - * @date 27 Mar 2025 + * @date 30 Apr 2025 * @author Tim Spain * @author Robert Jendersie */ @@ -397,13 +397,13 @@ template void CGDynamicsKernel::stressDivergence( template void CGDynamicsKernel::dirichletZero(CGVector& v) const { - // TR 07.04.2025: Dirichlet Zero (u=v=0) holds on land and on the boundary betreen - // land and ice. Hence on all elements with landmask = 0, or, on cg nodes with - // cglandmask = 0 + // TR 07.04.2025: Dirichlet Zero (u=v=0) holds on land and on the boundary betreen + // land and ice. Hence on all elements with landmask = 0, or, on cg nodes with + // cglandmask = 0 #pragma omp parallel for - for (size_t i = 0; i < pmap->cglandmask.rows(); ++i) - if (pmap->cglandmask(i) == 0) // land - v(i) = 0; + for (size_t i = 0; i < pmap->cglandmask.rows(); ++i) + if (pmap->cglandmask(i) == 0) // land + v(i) = 0; } template void CGDynamicsKernel::applyBoundaries() diff --git a/dynamics/src/DGTransport.cpp b/dynamics/src/DGTransport.cpp index 5e2b0ac64..12b1311a2 100644 --- a/dynamics/src/DGTransport.cpp +++ b/dynamics/src/DGTransport.cpp @@ -1,6 +1,6 @@ /*! * @file DGTransport.cpp - * @date July 10, 2022 + * @date 30 Apr 2025 * @author Thomas Richter */ @@ -232,18 +232,16 @@ template void DGTransport::reinitnormalvelocity() // only has one neighbor. Above, we sum both sides, weighted with 0.5. Therefore, // we must scale the outer edges with 2. This is only used for inflow/outflow // Neumann boundaries. -#pragma omp parallel for - for (size_t i=0;i::DGTransportOperator(const ParametricMesh& smesh, const dou } // TR 07.04.2025 - // + // // Dirichlet: we do not need to do anything special in DGTransport. On Dirichlet // boundaries, the normal velocity is zero. Hence, any flux is zero. @@ -485,35 +483,33 @@ void DGTransport::DGTransportOperator(const ParametricMesh& smesh, const dou // This can only appear on outer mesh-boundaries, whenever the element at the edge // is ice. - // TO DISCUSS: there can't be in/outflow and periodic on the same edge. We either need a structure - // for in/outflow or we say: either there is an in/outflow boundary or periodic. But not both.. - if (smesh.periodic.size()==0) - { + // TO DISCUSS: there can't be in/outflow and periodic on the same edge. We either need a + // structure for in/outflow or we say: either there is an in/outflow boundary or periodic. But + // not both.. + if (smesh.periodic.size() == 0) { #pragma omp parallel for - for (size_t i=0; i */ @@ -116,14 +116,14 @@ void ParametricMesh::readmesh(std::string fname) << std::endl; abort(); } - std::cerr << "ParametricMesh V2.0 is not longer fully supported. Dirichlet now uses the landmask information and Dirichlet data within the mesh file is ignored" << std::endl; + std::cerr << "ParametricMesh V2.0 is not longer fully supported. Dirichlet now uses the " + "landmask information and Dirichlet data within the mesh file is ignored" + << std::endl; size_t nd; IN >> nd; - double tmp; - for (int i=0;i> tmp >> tmp; - - + double tmp; + for (int i = 0; i < nd; ++i) + IN >> tmp >> tmp; IN >> status; if (status != "periodic") { diff --git a/dynamics/src/include/ParametricMesh.hpp b/dynamics/src/include/ParametricMesh.hpp index b100d3628..42a451908 100644 --- a/dynamics/src/include/ParametricMesh.hpp +++ b/dynamics/src/include/ParametricMesh.hpp @@ -1,6 +1,6 @@ /*! * @file ParametricMesh.hpp - * @date 09 Jan 2025 + * @date 30 Apr 2025 * @author Thomas Richter */ @@ -378,7 +378,6 @@ class ParametricMesh { */ void landmaskFromModelArray(const ModelArray& mask); - // Global access functions double hmin() const; //! returns the minimum mesh size From 2e06954b201c8a6c9c4542801843aaa309fcb149 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Einar=20=C3=96rn=20=C3=93lason?= Date: Wed, 30 Apr 2025 11:52:47 +0200 Subject: [PATCH 15/33] Increase the domain size This is the same now as in Bjornsson et al. --- run/make_init_polynya.py | 7 ++----- 1 file changed, 2 insertions(+), 5 deletions(-) diff --git a/run/make_init_polynya.py b/run/make_init_polynya.py index 7704a01df..40079566c 100644 --- a/run/make_init_polynya.py +++ b/run/make_init_polynya.py @@ -3,11 +3,8 @@ # Creates initial conditions for the Bjornsson et al. (2001) polynya case # Domain size [km] -#x = 100 -#y = 50 -#res = 2 -x = 64 -y = 48 +x = 100 +y = 50 res = 2 From 6bcb88d0fec2e5bc0ce972a724fd6f5f553eae64 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Einar=20=C3=96rn=20=C3=93lason?= Date: Wed, 30 Apr 2025 11:57:33 +0200 Subject: [PATCH 16/33] Update PrognosticData_test.cpp ConfiguredAtmosphere doesn't accept wind_speed as an option anymore, only wind_u and wind_v, so I set them to 3 and 4 to get a wind speed of 5. Also removed an unused header file. --- core/test/PrognosticData_test.cpp | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/core/test/PrognosticData_test.cpp b/core/test/PrognosticData_test.cpp index 4b5ce225d..d1869a9f6 100644 --- a/core/test/PrognosticData_test.cpp +++ b/core/test/PrognosticData_test.cpp @@ -1,7 +1,7 @@ /*! * @file PrognosticData_test.cpp * - * @date 7 Sep 2023 + * @date 30 Apr 2025 * @author Tim Spain */ @@ -17,7 +17,6 @@ #include "include/constants.hpp" #include -#include extern template class Module::Module; @@ -41,7 +40,8 @@ TEST_CASE("PrognosticData call order test") config << "lw_in = 330" << std::endl; config << "snow = 0" << std::endl; config << "rainfall = 0" << std::endl; - config << "wind_speed = 5" << std::endl; + config << "wind_u = 3" << std::endl; + config << "wind_v = 4" << std::endl; std::unique_ptr pcstream(new std::stringstream(config.str())); Configurator::addStream(std::move(pcstream)); From 437d6d1ed7294b36eb598f09051b86d3f29a2439 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Einar=20=C3=96rn=20=C3=93lason?= Date: Wed, 30 Apr 2025 12:01:14 +0200 Subject: [PATCH 17/33] Remove commented-out code It's not needed any more! Also some clang-format changes. --- dynamics/test/AdvectionPeriodicBC_test.cpp | 42 ++++++-------- dynamics/test/Advection_test.cpp | 65 +++++++--------------- 2 files changed, 38 insertions(+), 69 deletions(-) diff --git a/dynamics/test/AdvectionPeriodicBC_test.cpp b/dynamics/test/AdvectionPeriodicBC_test.cpp index b3e5b6e5c..a53d47cbf 100644 --- a/dynamics/test/AdvectionPeriodicBC_test.cpp +++ b/dynamics/test/AdvectionPeriodicBC_test.cpp @@ -1,6 +1,6 @@ /*! * @file AdvectionPeriodicBC_test.cpp - * @date 19 August 2024 + * @date 30 Apr 2025 * @author Thomas Richter */ @@ -237,22 +237,22 @@ void create_rectanglemesh(Nextsim::ParametricMesh& smesh) // landmask: set all to ice smesh.landmask.resize(smesh.nelements, 1); - //// Boundary. Dirichlet on bottom (by landmask) / top, periodic left/right, + //// Boundary. Dirichlet on bottom (by landmask) / top, periodic left/right, for (size_t i = 0; i < smesh.nx; ++i) smesh.landmask[i] = 0; // bot for (size_t i = 0; i < smesh.nx; ++i) - smesh.landmask[i + smesh.nx * (smesh.ny-1) ] = 0; + smesh.landmask[i + smesh.nx * (smesh.ny - 1)] = 0; smesh.periodic.resize(1); - smesh.periodic[0].resize(smesh.ny-2); // leave out top and bottom element - for (size_t i = 0; i < smesh.ny-2; ++i) { + smesh.periodic[0].resize(smesh.ny - 2); // leave out top and bottom element + for (size_t i = 0; i < smesh.ny - 2; ++i) { smesh.periodic[0][i][0] = 1; // periodic on left/right boundary - smesh.periodic[0][i][1] = (i+1 + 1) * smesh.nx + smesh.periodic[0][i][1] = (i + 1 + 1) * smesh.nx - 1; // element on the left of the edge (so the most right one in the row) smesh.periodic[0][i][2] - = (i+1) * smesh.nx; // element right of the edge (so the very left one) as... + = (i + 1) * smesh.nx; // element right of the edge (so the very left one) as... smesh.periodic[0][i][3] - = (i+2) * (smesh.nx+1) - 1; // the edge is the most right edge in the row + = (i + 2) * (smesh.nx + 1) - 1; // the edge is the most right edge in the row } } @@ -280,25 +280,19 @@ template void run(const std::array, 3>& exact) TEST_SUITE_BEGIN("Advection Periodic Boundary Conditions"); TEST_CASE("Advection Periodic Boundary Conditions") { - // std::array, 3> exact = // Exact values taken 26/06/2024 - // { std::array( - // { 1.0338503986019776e+00, 1.1451366598186583e+00, 1.0681593193338459e+00, - // 9.5252231195653603e-01, 8.1458581892610615e-01, 6.8950068528265651e-01 }), - // std::array( - // { 1.0949680727313791e+00, 8.3851748134278226e-01, 5.8501741891364523e-01, - // 4.0006092999256893e-01, 2.8493698219799068e-01, 2.0801424342840474e-01 }), - // std::array( - // { 7.1704413076190610e-01, 4.6135258391583606e-01, 3.2770311091970727e-01, - // 2.3424633076579465e-01, 1.7038068017215552e-01, 1.2486932724366147e-01 }) }; - // New values 07/04/25 due to new way of handling boundaries. We now include a land-layer // at the top and bottom. The slightly changed domain leads to changed values std::array, 3> exact = // Exact values taken 26/06/2024 - { std::array({1.0338503986019776e+00,1.1451366598186583e+00,1.0681593193338459e+00,9.5252231195653614e-01,8.1458581892610615e-01,6.8950068528265651e-01}), - std::array({8.7678853597311990e-01,8.2389349440313742e-01,5.8562737114004948e-01,4.0006113872450388e-01,2.8493698219799091e-01,2.0801424342840488e-01}), - std::array({6.0231338733029061e-01,4.6967937778236674e-01,3.2769774426726950e-01,2.3427321065537998e-01,1.7038111756667601e-01,1.2487338791917807e-01}) - }; - + { std::array( + { 1.0338503986019776e+00, 1.1451366598186583e+00, 1.0681593193338459e+00, + 9.5252231195653614e-01, 8.1458581892610615e-01, 6.8950068528265651e-01 }), + std::array( + { 8.7678853597311990e-01, 8.2389349440313742e-01, 5.8562737114004948e-01, + 4.0006113872450388e-01, 2.8493698219799091e-01, 2.0801424342840488e-01 }), + std::array( + { 6.0231338733029061e-01, 4.6967937778236674e-01, 3.2769774426726950e-01, + 2.3427321065537998e-01, 1.7038111756667601e-01, 1.2487338791917807e-01 }) }; + std::cout << "DG\tNT\tNX\tmass\t\terror\t\texact" << std::endl; std::cout << std::setprecision(4) << std::scientific; run<1>(exact); diff --git a/dynamics/test/Advection_test.cpp b/dynamics/test/Advection_test.cpp index 6f95edc0c..feadb5c7e 100644 --- a/dynamics/test/Advection_test.cpp +++ b/dynamics/test/Advection_test.cpp @@ -1,6 +1,6 @@ /*! * @file Advection_test.cpp - * @date 19 August 2024 + * @date 30 Apr 2025 * @author Thomas Richter */ @@ -228,11 +228,11 @@ void create_rectanglemesh(Nextsim::ParametricMesh& smesh, const double Lx, const for (size_t i = 0; i < smesh.nx; ++i) smesh.landmask[i] = 0; for (size_t i = 0; i < smesh.nx; ++i) - smesh.landmask[smesh.nx * (smesh.ny-1) + i] = 0; + smesh.landmask[smesh.nx * (smesh.ny - 1) + i] = 0; for (size_t i = 0; i < smesh.ny; ++i) smesh.landmask[i * smesh.nx] = 0; for (size_t i = 0; i < smesh.ny; ++i) - smesh.landmask[(i+1) * smesh.nx-1 ] = 0; + smesh.landmask[(i + 1) * smesh.nx - 1] = 0; } template void run(double distort, const std::array, 3>& exact) @@ -264,31 +264,17 @@ TEST_CASE("Advection") { std::cout << std::setprecision(4) << std::scientific; - // Exact values taken 25/06/2024 after some plausibility check of the results - // and by checking the theoretical order of convergence to be expected - // std::array, 3> exact - // = { std::array( - // { 5.1256149074257538e-02, 4.8288256703303903e-02, 4.3105635248809886e-02, - // 3.5793986049422598e-02, 2.6937487824223016e-02, 1.8456775604583933e-02 }), - // std::array( - // { 2.9450967798560313e-02, 1.3427281824939470e-02, 5.8574889800512000e-03, - // 2.2756813704797301e-03, 7.3564079504566610e-04, 1.9177169540867740e-04 }), - // std::array( - // { 9.9340386651904228e-03, 4.0274889287136816e-03, 1.2400186092664943e-03, - // 2.8460130790583933e-04, 4.5459555769938981e-05, 4.6851425365137733e-06 }) }; - - // New values 07/04/2025 after reimplementation of boundary data. Dirichlet data is now - // realized using the landmaks info. Therefore, the physical domain shrinks by one layer all around - // and this slightly changes the exact values. + // realized using the landmaks info. Therefore, the physical domain shrinks by one layer all + // around and this slightly changes the exact values. std::array, 3> exact - = { std::array( - { 5.0568601321622546e-02,4.8785010758681462e-02,4.3807523095968269e-02,3.6079575213776381e-02 }), - std::array( - { 2.8463796850208754e-02,1.4000140482658283e-02,5.8605850177331619e-03,2.2757018442945147e-03 }), - std::array( - { 1.0146940898974235e-02,4.0853715304369851e-03,1.2400975498406005e-03,2.8460132054714237e-04}) }; - + = { std::array({ 5.0568601321622546e-02, 4.8785010758681462e-02, + 4.3807523095968269e-02, 3.6079575213776381e-02 }), + std::array({ 2.8463796850208754e-02, 1.4000140482658283e-02, + 5.8605850177331619e-03, 2.2757018442945147e-03 }), + std::array({ 1.0146940898974235e-02, 4.0853715304369851e-03, + 1.2400975498406005e-03, 2.8460132054714237e-04 }) }; + std::cout << std::endl << "DG\tNT\tNX\tmass loss\terror\t\texact" << std::endl; run<1>(0.0, exact); run<3>(0.0, exact); @@ -296,27 +282,16 @@ TEST_CASE("Advection") } TEST_CASE("Distorted Mesh") { - // std::array, 3> exact - // = { std::array( - // { 5.0958748236875594e-02, 4.8461087243594887e-02, 4.3918572621094686e-02, - // 3.7038043078562323e-02, 2.8334464207979679e-02, 1.9666196904181983e-02 }), - // std::array( - // { 3.2033423984965226e-02, 1.5639052766007147e-02, 6.7926997203524983e-03, - // 2.7369916393700151e-03, 9.2109964949992139e-04, 2.5128064874203957e-04 }), - // std::array( - // { 1.1830071142946669e-02, 4.9207680503709564e-03, 1.5831512504162868e-03, - // 3.8869595113351363e-04, 6.7162895595772516e-05, 7.5853786764529606e-06 }) }; - // New values 07/04/2025 after reimplementation of boundary data. Dirichlet data is now - // realized using the landmaks info. Therefore, the physical domain shrinks by one layer all around - // and this slightly changes the exact values. + // realized using the landmaks info. Therefore, the physical domain shrinks by one layer all + // around and this slightly changes the exact values. std::array, 3> exact - = { std::array( - { 4.8195942137802428e-02,4.9345757529052611e-02,4.4903895856787640e-02,3.7467928845481245e-02 }), - std::array( - { 3.2422647655365219e-02,1.7062569520319999e-02,6.8213111669173628e-03,2.7372092446537699e-03 }), - std::array( - { 1.2603978369132178e-02,5.1434197544894689e-03,1.5845944976351275e-03,3.8869603127842235e-04}) }; + = { std::array({ 4.8195942137802428e-02, 4.9345757529052611e-02, + 4.4903895856787640e-02, 3.7467928845481245e-02 }), + std::array({ 3.2422647655365219e-02, 1.7062569520319999e-02, + 6.8213111669173628e-03, 2.7372092446537699e-03 }), + std::array({ 1.2603978369132178e-02, 5.1434197544894689e-03, + 1.5845944976351275e-03, 3.8869603127842235e-04 }) }; std::cout << std::endl << "Distorted mesh" << std::endl; std::cout << "DG\tNT\tNX\tmass loss\terror\t\texact" << std::endl; From 5781dbfbb98d10afe3561d225c7c896eab474d12 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Einar=20=C3=96rn=20=C3=93lason?= Date: Tue, 20 May 2025 16:02:03 +0200 Subject: [PATCH 18/33] Minor fixes to get the model to compile and polynya simulation to start But it still crashes with unstabilities on the boundary after a few time steps. --- .../ConfiguredAtmosphere.cpp | 5 ++- .../FluxConfiguredAtmosphere.cpp | 38 +++++++++---------- run/make_init_polynya.py | 4 +- 3 files changed, 25 insertions(+), 22 deletions(-) diff --git a/physics/src/modules/AtmosphereBoundaryModule/ConfiguredAtmosphere.cpp b/physics/src/modules/AtmosphereBoundaryModule/ConfiguredAtmosphere.cpp index 9b49ede5b..2a3fed067 100644 --- a/physics/src/modules/AtmosphereBoundaryModule/ConfiguredAtmosphere.cpp +++ b/physics/src/modules/AtmosphereBoundaryModule/ConfiguredAtmosphere.cpp @@ -1,7 +1,7 @@ /*! * @file ConfiguredAtmosphere.cpp * - * @date 30 Apr 2025 + * @date 20 May 2025 * @author Tim Spain */ @@ -112,7 +112,8 @@ ConfigMap ConfiguredAtmosphere::getConfiguration() const { keyMap.at(LW_KEY), lw0 }, { keyMap.at(SNOW_KEY), snowfall0 }, { keyMap.at(RAIN_KEY), rain0 }, - { keyMap.at(WIND_KEY), windspeed0 }, + { keyMap.at(UWIND_KEY), uWind0 }, + { keyMap.at(VWIND_KEY), vWind0 }, }; } diff --git a/physics/src/modules/AtmosphereBoundaryModule/FluxConfiguredAtmosphere.cpp b/physics/src/modules/AtmosphereBoundaryModule/FluxConfiguredAtmosphere.cpp index 2b8db8ff0..1562b1309 100644 --- a/physics/src/modules/AtmosphereBoundaryModule/FluxConfiguredAtmosphere.cpp +++ b/physics/src/modules/AtmosphereBoundaryModule/FluxConfiguredAtmosphere.cpp @@ -1,7 +1,7 @@ /*! * @file FluxConfiguredAtmosphere.cpp * - * @date 30 Apr 2025 + * @date 20 May 2025 * @author Tim Spain */ @@ -68,29 +68,29 @@ ConfigurationHelp::HelpMap& FluxConfiguredAtmosphere::getHelpRecursive(HelpMap& } void FluxConfiguredAtmosphere::configure() { - qia0 = getConfiguration(keyMap.at(QIA_KEY), qia0); - dqia_dt0 = getConfiguration(keyMap.at(DQIA_DT_KEY), dqia_dt0); - qow0 = getConfiguration(keyMap.at(QOW_KEY), qow0); - subl0 = getConfiguration(keyMap.at(SUBL_KEY), subl0); - snowfall0 = getConfiguration(keyMap.at(SNOW_KEY), snowfall0); - rain0 = getConfiguration(keyMap.at(RAIN_KEY), rain0); - evap0 = getConfiguration(keyMap.at(EVAP_KEY), evap0); - u0 = getConfiguration(keyMap.at(WINDU_KEY), u0); - v0 = getConfiguration(keyMap.at(WINDV_KEY), v0); + qia0 = Configured::getConfiguration(keyMap.at(QIA_KEY), qia0); + dqia_dt0 = Configured::getConfiguration(keyMap.at(DQIA_DT_KEY), dqia_dt0); + qow0 = Configured::getConfiguration(keyMap.at(QOW_KEY), qow0); + subl0 = Configured::getConfiguration(keyMap.at(SUBL_KEY), subl0); + snowfall0 = Configured::getConfiguration(keyMap.at(SNOW_KEY), snowfall0); + rain0 = Configured::getConfiguration(keyMap.at(RAIN_KEY), rain0); + evap0 = Configured::getConfiguration(keyMap.at(EVAP_KEY), evap0); + u0 = Configured::getConfiguration(keyMap.at(WINDU_KEY), u0); + v0 = Configured::getConfiguration(keyMap.at(WINDV_KEY), v0); } ConfigMap FluxConfiguredAtmosphere::getConfiguration() const { return { - { keyMap.at(QIA_KEY), qia0 }, - { keyMap.at(DQIA_DT_KEY), dqia_dt0 }, - { keyMap.at(QOW_KEY), qow0 }, - { keyMap.at(SUBL_KEY), subl0 }, - { keyMap.at(SNOW_KEY), snowfall0 }, - { keyMap.at(RAIN_KEY), rain0 }, - { keyMap.at(EVAP_KEY), evap0 }, - { keyMap.at(WINDU_KEY), u0 }, - { keyMap.at(WINDV_KEY), v0 }, + { keyMap.at(QIA_KEY), qia0 }, + { keyMap.at(DQIA_DT_KEY), dqia_dt0 }, + { keyMap.at(QOW_KEY), qow0 }, + { keyMap.at(SUBL_KEY), subl0 }, + { keyMap.at(SNOW_KEY), snowfall0 }, + { keyMap.at(RAIN_KEY), rain0 }, + { keyMap.at(EVAP_KEY), evap0 }, + { keyMap.at(WINDU_KEY), u0 }, + { keyMap.at(WINDV_KEY), v0 }, }; } diff --git a/run/make_init_polynya.py b/run/make_init_polynya.py index 40079566c..6ce4c2848 100644 --- a/run/make_init_polynya.py +++ b/run/make_init_polynya.py @@ -41,7 +41,9 @@ initializer.sss[:, :] = ocean_salinity initializer.sst[:, :] = ocean_temperature -initializer.tice[:, :, :] = ice_salinity * mu +initializer.tsurf[:, :] = ice_salinity * mu +initializer.tbott = initializer.tsurf +initializer.tintr = initializer.tsurf # All other variables are zero or not needed From 0ad9d3237a2bcbd342b3472bf7f23633020faefa Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Einar=20=C3=96rn=20=C3=93lason?= Date: Fri, 30 May 2025 08:02:01 +0200 Subject: [PATCH 19/33] Revert "Minor tweaks and fixes" This reverts commit 768c4d4c --- physics/src/IceGrowth.cpp | 18 +++++++++--------- 1 file changed, 9 insertions(+), 9 deletions(-) diff --git a/physics/src/IceGrowth.cpp b/physics/src/IceGrowth.cpp index 14408e25d..8f0d2feac 100644 --- a/physics/src/IceGrowth.cpp +++ b/physics/src/IceGrowth.cpp @@ -1,7 +1,7 @@ /*! * @file IceGrowth.cpp * - * @date 30 Apr 2025 + * @date 20 Nov 2024 * @author Tim Spain * @author Einar Ólason */ @@ -207,22 +207,22 @@ void IceGrowth::newIceFormation(size_t i, const TimestepTime& tst) { // Flux cooling the ocean from open water // TODO Add assimilation fluxes here - const double coolingFlux = qow[i]; + double coolingFlux = qow[i]; // Temperature change of the mixed layer during this timestep - const double deltaTml = -coolingFlux / mixedLayerBulkHeatCapacity[i] * tst.step; + double deltaTml = -coolingFlux / mixedLayerBulkHeatCapacity[i] * tst.step; // Initial temperature - const double t0 = sst[i]; + double t0 = sst[i]; // Freezing point temperature - const double tf0 = tf[i]; + double tf0 = tf[i]; // Final temperature - const double t1 = t0 + deltaTml; + double t1 = t0 + deltaTml; // deal with cooling below the freezing point if (t1 < tf0) { // Heat lost cooling the mixed layer to freezing point - const double sensibleFlux = -(tf0 - t0) * mixedLayerBulkHeatCapacity[i] / tst.step; + double sensibleFlux = (tf0 - t0) / deltaTml * coolingFlux; // Any heat beyond that is latent heat forming new ice - const double latentFlux = coolingFlux - sensibleFlux; + double latentFlux = coolingFlux - sensibleFlux; qow[i] = sensibleFlux; newice[i] = latentFlux * tst.step * (1 - cice[i]) / (Ice::Lf * Ice::rho); @@ -264,7 +264,7 @@ void IceGrowth::lateralIceSpread(size_t i, const TimestepTime& tstep) void IceGrowth::applyLimits(size_t i, const TimestepTime& tstep) { - if (cice[i] < IceMinima::c() || hice[i] < IceMinima::h()) { + if ((0. < cice[i] && cice[i] < IceMinima::c()) || (0. < hice[i] && hice[i] < IceMinima::h())) { qow[i] += cice[i] * Water::Lf * (hice[i] * Ice::rho + hsnow[i] * Ice::rhoSnow) / tstep.step; hice[i] = 0; cice[i] = 0; From 95c2a86b65fefc2f053a0a8b392c9ab4ecfdf29a Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Einar=20=C3=96rn=20=C3=93lason?= Date: Fri, 6 Jun 2025 09:07:55 +0200 Subject: [PATCH 20/33] Change back dates on several files With the latest merges, the date is the only difference between this branch and develop ... and that doesn't make sense. --- dynamics/src/CGDynamicsKernel.cpp | 2 +- dynamics/src/DGTransport.cpp | 2 +- dynamics/src/ParametricMesh.cpp | 2 +- dynamics/src/include/ParametricMesh.hpp | 2 +- 4 files changed, 4 insertions(+), 4 deletions(-) diff --git a/dynamics/src/CGDynamicsKernel.cpp b/dynamics/src/CGDynamicsKernel.cpp index 012890fe0..6938a6ee8 100644 --- a/dynamics/src/CGDynamicsKernel.cpp +++ b/dynamics/src/CGDynamicsKernel.cpp @@ -1,7 +1,7 @@ /*! * @file CGDynamicsKernel.cpp * - * @date 30 Apr 2025 + * @date 27 Mar 2025 * @author Tim Spain * @author Robert Jendersie */ diff --git a/dynamics/src/DGTransport.cpp b/dynamics/src/DGTransport.cpp index 12b1311a2..10475f300 100644 --- a/dynamics/src/DGTransport.cpp +++ b/dynamics/src/DGTransport.cpp @@ -1,6 +1,6 @@ /*! * @file DGTransport.cpp - * @date 30 Apr 2025 + * @date July 10, 2025 * @author Thomas Richter */ diff --git a/dynamics/src/ParametricMesh.cpp b/dynamics/src/ParametricMesh.cpp index 5637062c3..1e234be44 100644 --- a/dynamics/src/ParametricMesh.cpp +++ b/dynamics/src/ParametricMesh.cpp @@ -1,6 +1,6 @@ /*! * @file ParametricMesh.hpp - * @date 30 Apr 2025 + * @date 14 Jan 2025 * @author Thomas Richter */ diff --git a/dynamics/src/include/ParametricMesh.hpp b/dynamics/src/include/ParametricMesh.hpp index 42a451908..d78d34b81 100644 --- a/dynamics/src/include/ParametricMesh.hpp +++ b/dynamics/src/include/ParametricMesh.hpp @@ -1,6 +1,6 @@ /*! * @file ParametricMesh.hpp - * @date 30 Apr 2025 + * @date 09 Jan 2025 * @author Thomas Richter */ From fad93a1e1fd647c092bd3fcd24b16c2780b6ed34 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Einar=20=C3=96rn=20=C3=93lason?= Date: Fri, 6 Jun 2025 09:09:05 +0200 Subject: [PATCH 21/33] Fix the date fixing Got the year wrong on one of the files --- dynamics/src/DGTransport.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/dynamics/src/DGTransport.cpp b/dynamics/src/DGTransport.cpp index 10475f300..2c2661b4f 100644 --- a/dynamics/src/DGTransport.cpp +++ b/dynamics/src/DGTransport.cpp @@ -1,6 +1,6 @@ /*! * @file DGTransport.cpp - * @date July 10, 2025 + * @date July 10, 2022 * @author Thomas Richter */ From a783559ffa671407d55efd12bd5d751145a8f5dc Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Einar=20=C3=96rn=20=C3=93lason?= Date: Sun, 29 Jun 2025 00:10:36 +0200 Subject: [PATCH 22/33] Fix polynya initialisation The call to initMaker in make_init_polynya.py was still using nLayers, and not isWinton. --- run/make_init_polynya.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/run/make_init_polynya.py b/run/make_init_polynya.py index 6ce4c2848..2079501db 100644 --- a/run/make_init_polynya.py +++ b/run/make_init_polynya.py @@ -15,7 +15,7 @@ fname = f"init_polynya.nc" # The model expects everything in metres -initializer = initMaker(fname, nfirst, nsecond, nLayers, res*1e3, checkZeros=False) +initializer = initMaker(fname, nfirst, nsecond, res*1e3, isWinton=True, checkZeros=False) # Ice everywhere and all boundaries closed, except the x = 100 km end initializer.mask[:, :] = 1. From 7d4da1fbd4a475e0ef29841d2e1eeae6ecb59ec2 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Einar=20=C3=96rn=20=C3=93lason?= Date: Sun, 29 Jun 2025 07:08:21 +0200 Subject: [PATCH 23/33] Make an integration test out of the polynya case I reduce the resolution, increase the time step, and run it for five days. This results in a steady state. I then simply compare the max, min, and mean of hice, cice, u, and v with "known good values". --- .github/workflows/test_suite.yml | 2 + test/PolynyaIntegration_test.py | 196 +++++++++++++++++++++++++++++++ 2 files changed, 198 insertions(+) create mode 100644 test/PolynyaIntegration_test.py diff --git a/.github/workflows/test_suite.yml b/.github/workflows/test_suite.yml index 138821507..909714bae 100644 --- a/.github/workflows/test_suite.yml +++ b/.github/workflows/test_suite.yml @@ -70,6 +70,7 @@ jobs: mv ../test . cd test python3 ThermoIntegration_test.py + python3 PolynyaIntegration_test.py ./run_test_jan2010_integration_test.sh python3 test-ubuntu-mpi-noxios: @@ -178,4 +179,5 @@ jobs: mv ../test . cd test python ThermoIntegration_test.py + python3 PolynyaIntegration_test.py ./run_test_jan2010_integration_test.sh python diff --git a/test/PolynyaIntegration_test.py b/test/PolynyaIntegration_test.py new file mode 100644 index 000000000..783df6b68 --- /dev/null +++ b/test/PolynyaIntegration_test.py @@ -0,0 +1,196 @@ +import os +import subprocess +import sys +import unittest + +import netCDF4 +import numpy as np + + +class Polynya(unittest.TestCase): + # A few useful global variables for the class + executable = "../nextsim" + + init_file = "init_polynya.nc" + config_file = "PolynyaIntegration.cfg" + diagnostics_file = "polynya.diagnostic.nc" + + # Load the data once and re-use for all tests + hice = np.array([]) + cice = np.array([]) + uice = np.array([]) + vice = np.array([]) + + @classmethod + def setUpClass(cls): + """ + A set-up class which, + - Creates the initialisation file, using make_init_column.py + - Runs the model + - Loads the neccesary variables from the output file + """ + + # Make the init column + cls.__make_init_column() + + # Create the config file + cls.__make_cfg_file() + + # Run the model + subprocess.run(cls.executable + " --config-file " + cls.config_file, shell=True, check=True) + + # Load the basic variables + root = netCDF4.Dataset(cls.diagnostics_file, "r", format="NETCDF4") + cls.cice = np.squeeze(np.array(root.groups["data"].variables["cice"][:].data)) + cls.hice = np.squeeze(np.array(root.groups["data"].variables["hice"][:].data)) + cls.uice = np.array(root.groups["data"].variables["u"][:].data) + cls.vice = np.array(root.groups["data"].variables["v"][:].data) + + @classmethod + def __make_cfg_file(cls): + cfg = open(cls.config_file, "w") + cfg.write(""" +[model] +init_file = init_polynya.nc +start = 2023-01-01T00:00:00Z +stop = 2023-01-05T00:00:00Z +time_step = P0-0T00:30:00 +missing_value = 1e20 + +[Modules] +DiagnosticOutputModule = Nextsim::ConfigOutput +DynamicsModule = Nextsim::BBMDynamics +AtmosphereBoundaryModule = Nextsim::FluxConfiguredAtmosphere +OceanBoundaryModule = Nextsim::FluxConfiguredOcean +IceThermodynamicsModule = Nextsim::ThermoWinton + +[ConfigOutput] +period = P0-0T00:01:00 +field_names = hice,cice,u,v +filename = polynya.diagnostic.nc + +[FluxConfiguredOcean] +qio = 2.30 +sss = 28 +sst = -1.54 +mld = 10. +current_u = 0. +current_v = 0. + +[FluxConfiguredAtmosphere] +Q_ia = 30 +Q_ow = 300 +dQia_dT = 0 +wind_u = 16 +wind_v = 12 + """) + cfg.close() + + @classmethod + def __make_init_column(cls): + sys.path.append('../run') + from make_init_base import initMaker + + # Creates initial conditions for the Bjornsson et al. (2001) polynya case + + # Domain size [km] + x = 100 + y = 50 + res = 4 + + + nfirst = int(y / res) + nsecond = int(x / res) + + fname = f"init_polynya.nc" + + # The model expects everything in metres + initializer = initMaker(fname, nfirst, nsecond, res*1e3, isWinton=True, checkZeros=False) + + # Ice everywhere and all boundaries closed, except the x = 100 km end + initializer.mask[:, :] = 1. + initializer.mask[0, :] = 0. + initializer.mask[-1, :] = 0. + initializer.mask[:, 0] = 0. + #initializer.mask[:, -1] = 0. ## right + + # Uniform concentration of 90% + initializer.cice[:, :] = 0.9 + + # Uniform thickness of 20 cm + initializer.hice[:, :] = 0.2 + + # Undamaged ice + initializer.damage[:, :] = 1. + + # Ice and ocean temperature and salinity at the freezing point + ice_salinity = 5 # should match Ice::s in constants.hpp + mu: float = -0.055 # should match Water::mu in constants.hpp + ocean_temperature = -1.54 + ocean_salinity = ocean_temperature / mu + + initializer.sss[:, :] = ocean_salinity + initializer.sst[:, :] = ocean_temperature + initializer.tsurf[:, :] = ice_salinity * mu + initializer.tbott = initializer.tsurf + initializer.tintr = initializer.tsurf + + # All other variables are zero or not needed + + # The file is written when initializer goes out of scope + + @classmethod + def tearDownClass(cls): + """ + A tear-down class that deletes the netCDF output and temporary files + """ + + if os.path.isfile(cls.diagnostics_file): + os.remove(cls.diagnostics_file) + + if os.path.isfile(cls.init_file): + os.remove(cls.init_file) + + if os.path.isfile(cls.config_file): + os.remove(cls.config_file) + + def test_iceThickness(self): + """ + Test the ice thickness against standard max, min, and mean values + """ + + mean = 0.0301 + max = 2.3082 + min = -1.7963 + self.assertAlmostEqual(max, self.hice.max(), 4, "Max ice thickness not ~= " + str(max) + " m (full DG field)") + self.assertAlmostEqual(min, self.hice.min(), 4, "Min ice thickness not ~= " + str(min) + " m (full DG field)") + self.assertAlmostEqual(mean, self.hice.mean(), 4, "Mean ice thickness not ~= " + str(mean) + " m (full DG field)") + + def test_concentration(self): + """ + Test the ice concentration against standard max, min, and mean values + """ + + mean = 0.0931 + max = 2.2597 + min = -1.6293 + self.assertAlmostEqual(max, self.cice.max(), 4, "Max conentration not ~= " + str(max) + " (full DG field)") + self.assertAlmostEqual(min, self.cice.min(), 4, "Min concentration not ~= " + str(min) + " (full DG field)") + self.assertAlmostEqual(mean, self.cice.mean(), 4, "Mean concentration not ~= " + str(mean) + " (full DG field)") + + def test_veloities(self): + """ + Test the ice velocities against standard max, min, and mean values + """ + + mean = [0.2237, 0.0548] + max = [0.3910, 0.1305] + min = [-0.0700, -0.1559] + for (var, i, name) in zip([self.uice, self.vice], [0,1], ["uice", "vice"]): + self.assertAlmostEqual(max[i], var.max(), 4, "Max "+name+" velocity not ~= " + str(max[i]) + " m/s (full DG field)") + self.assertAlmostEqual(min[i], var.min(), 4, "Min "+name+" velocity not ~= " + str(min[i]) + " m/s (full DG field)") + self.assertAlmostEqual(mean[i], var.mean(), 4, "Mean "+name+" velocity not ~= " + str(mean[i]) + " m/s (full DG field)") + + +if __name__ == '__main__': + unittest.main() From dd519cddafe5b7528f938f86b50213a0d1f10354 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Einar=20=C3=96rn=20=C3=93lason?= Date: Sun, 29 Jun 2025 07:17:53 +0200 Subject: [PATCH 24/33] Additional python path for the CI The setup is different over there. --- test/PolynyaIntegration_test.py | 1 + 1 file changed, 1 insertion(+) diff --git a/test/PolynyaIntegration_test.py b/test/PolynyaIntegration_test.py index 783df6b68..b09607691 100644 --- a/test/PolynyaIntegration_test.py +++ b/test/PolynyaIntegration_test.py @@ -89,6 +89,7 @@ def __make_cfg_file(cls): @classmethod def __make_init_column(cls): sys.path.append('../run') + sys.path.append('../../run') from make_init_base import initMaker # Creates initial conditions for the Bjornsson et al. (2001) polynya case From 54919ec607030ad6f2a001f91ad76d15522f934d Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Einar=20=C3=96rn=20=C3=93lason?= Date: Sun, 29 Jun 2025 07:25:16 +0200 Subject: [PATCH 25/33] Fix wrong python command on macOS Copy/paste error and python vs python3 --- .github/workflows/test_suite.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.github/workflows/test_suite.yml b/.github/workflows/test_suite.yml index 909714bae..02fefb64c 100644 --- a/.github/workflows/test_suite.yml +++ b/.github/workflows/test_suite.yml @@ -179,5 +179,5 @@ jobs: mv ../test . cd test python ThermoIntegration_test.py - python3 PolynyaIntegration_test.py + python PolynyaIntegration_test.py ./run_test_jan2010_integration_test.sh python From 58d8277b1dae0b201219911ef13a06d317a9e98b Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Einar=20=C3=96rn=20=C3=93lason?= Date: Sun, 29 Jun 2025 07:36:52 +0200 Subject: [PATCH 26/33] Reduce accuracy of polynya integration test This is much more sensitive than I expected. So, I'm only using the DG0 component and not checking the velocities. Running it on macOS vs Linux gives fairly different maximum values for velocities and higher order DG fields. --- test/PolynyaIntegration_test.py | 40 ++++++++++++--------------------- 1 file changed, 14 insertions(+), 26 deletions(-) diff --git a/test/PolynyaIntegration_test.py b/test/PolynyaIntegration_test.py index b09607691..3bdcb8e60 100644 --- a/test/PolynyaIntegration_test.py +++ b/test/PolynyaIntegration_test.py @@ -160,38 +160,26 @@ def test_iceThickness(self): Test the ice thickness against standard max, min, and mean values """ - mean = 0.0301 - max = 2.3082 - min = -1.7963 - self.assertAlmostEqual(max, self.hice.max(), 4, "Max ice thickness not ~= " + str(max) + " m (full DG field)") - self.assertAlmostEqual(min, self.hice.min(), 4, "Min ice thickness not ~= " + str(min) + " m (full DG field)") - self.assertAlmostEqual(mean, self.hice.mean(), 4, "Mean ice thickness not ~= " + str(mean) + " m (full DG field)") + mean = 0.1317 + max = 1.1250 + min = 0.0000 + hice = self.hice[:,:,:,0] + self.assertAlmostEqual(max, hice.max(), 4, "Max ice thickness not ~= " + str(max) + " m") + self.assertAlmostEqual(min, hice.min(), 4, "Min ice thickness not ~= " + str(min) + " m") + self.assertAlmostEqual(mean, hice.mean(), 4, "Mean ice thickness not ~= " + str(mean) + " m") def test_concentration(self): """ Test the ice concentration against standard max, min, and mean values """ - mean = 0.0931 - max = 2.2597 - min = -1.6293 - self.assertAlmostEqual(max, self.cice.max(), 4, "Max conentration not ~= " + str(max) + " (full DG field)") - self.assertAlmostEqual(min, self.cice.min(), 4, "Min concentration not ~= " + str(min) + " (full DG field)") - self.assertAlmostEqual(mean, self.cice.mean(), 4, "Mean concentration not ~= " + str(mean) + " (full DG field)") - - def test_veloities(self): - """ - Test the ice velocities against standard max, min, and mean values - """ - - mean = [0.2237, 0.0548] - max = [0.3910, 0.1305] - min = [-0.0700, -0.1559] - for (var, i, name) in zip([self.uice, self.vice], [0,1], ["uice", "vice"]): - self.assertAlmostEqual(max[i], var.max(), 4, "Max "+name+" velocity not ~= " + str(max[i]) + " m/s (full DG field)") - self.assertAlmostEqual(min[i], var.min(), 4, "Min "+name+" velocity not ~= " + str(min[i]) + " m/s (full DG field)") - self.assertAlmostEqual(mean[i], var.mean(), 4, "Mean "+name+" velocity not ~= " + str(mean[i]) + " m/s (full DG field)") - + mean = 0.4402 + max = 1.0000 + min = 0.0000 + cice = self.cice[:,:,:,0] + self.assertAlmostEqual(max, cice.max(), 4, "Max conentration not ~= " + str(max)) + self.assertAlmostEqual(min, cice.min(), 4, "Min concentration not ~= " + str(min)) + self.assertAlmostEqual(mean, cice.mean(), 4, "Mean concentration not ~= " + str(mean)) if __name__ == '__main__': unittest.main() From ee55629c35327783af3e3652806e2b76dd338b7c Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Einar=20=C3=96rn=20=C3=93lason?= Date: Sun, 29 Jun 2025 07:47:07 +0200 Subject: [PATCH 27/33] Reduce the accuracy even further But it's still useful --- test/PolynyaIntegration_test.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/test/PolynyaIntegration_test.py b/test/PolynyaIntegration_test.py index 3bdcb8e60..f084d472a 100644 --- a/test/PolynyaIntegration_test.py +++ b/test/PolynyaIntegration_test.py @@ -158,19 +158,19 @@ def tearDownClass(cls): def test_iceThickness(self): """ Test the ice thickness against standard max, min, and mean values + The maximum varies by (at least) 3 cm between platforms, making a test of it useless """ mean = 0.1317 - max = 1.1250 min = 0.0000 hice = self.hice[:,:,:,0] - self.assertAlmostEqual(max, hice.max(), 4, "Max ice thickness not ~= " + str(max) + " m") self.assertAlmostEqual(min, hice.min(), 4, "Min ice thickness not ~= " + str(min) + " m") self.assertAlmostEqual(mean, hice.mean(), 4, "Mean ice thickness not ~= " + str(mean) + " m") def test_concentration(self): """ Test the ice concentration against standard max, min, and mean values + The min and max are easy, but the mean is sensitive """ mean = 0.4402 @@ -179,7 +179,7 @@ def test_concentration(self): cice = self.cice[:,:,:,0] self.assertAlmostEqual(max, cice.max(), 4, "Max conentration not ~= " + str(max)) self.assertAlmostEqual(min, cice.min(), 4, "Min concentration not ~= " + str(min)) - self.assertAlmostEqual(mean, cice.mean(), 4, "Mean concentration not ~= " + str(mean)) + self.assertAlmostEqual(mean, cice.mean(), 3, "Mean concentration not ~= " + str(mean)) if __name__ == '__main__': unittest.main() From d4be55e2dd160b7b777afac00135e64305c1dcb4 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Einar=20=C3=96rn=20=C3=93lason?= Date: Mon, 21 Jul 2025 09:03:13 +0200 Subject: [PATCH 28/33] Remove redundant f-string from make_init_polynya.py As suggested by Joe on github. Co-authored-by: Joe Wallwork <22053413+jwallwork23@users.noreply.github.com> --- run/make_init_polynya.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/run/make_init_polynya.py b/run/make_init_polynya.py index 2079501db..911e4b776 100644 --- a/run/make_init_polynya.py +++ b/run/make_init_polynya.py @@ -12,7 +12,7 @@ nsecond = int(x / res) nLayers = 3 -fname = f"init_polynya.nc" +fname = "init_polynya.nc" # The model expects everything in metres initializer = initMaker(fname, nfirst, nsecond, res*1e3, isWinton=True, checkZeros=False) From 17f55a72968045a5bbb42c68657c593a3aa26bc9 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Einar=20=C3=96rn=20=C3=93lason?= Date: Mon, 21 Jul 2025 09:16:46 +0200 Subject: [PATCH 29/33] Use integer division, instead of type conversion As suggested by Joe on github --- run/make_init_polynya.py | 5 ++--- test/PolynyaIntegration_test.py | 5 ++--- 2 files changed, 4 insertions(+), 6 deletions(-) diff --git a/run/make_init_polynya.py b/run/make_init_polynya.py index 911e4b776..ae91d3628 100644 --- a/run/make_init_polynya.py +++ b/run/make_init_polynya.py @@ -7,9 +7,8 @@ y = 50 res = 2 - -nfirst = int(y / res) -nsecond = int(x / res) +nfirst = y // res +nsecond = x // res nLayers = 3 fname = "init_polynya.nc" diff --git a/test/PolynyaIntegration_test.py b/test/PolynyaIntegration_test.py index f084d472a..889994c2d 100644 --- a/test/PolynyaIntegration_test.py +++ b/test/PolynyaIntegration_test.py @@ -99,9 +99,8 @@ def __make_init_column(cls): y = 50 res = 4 - - nfirst = int(y / res) - nsecond = int(x / res) + nfirst = y // res + nsecond = x // res fname = f"init_polynya.nc" From 00270c89985d13b35356ed0f30bfaec36355bf9a Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Einar=20=C3=96rn=20=C3=93lason?= Date: Mon, 21 Jul 2025 09:17:34 +0200 Subject: [PATCH 30/33] Use the class variable for init_file in PolynyaIntegration_test.py That's cleaner than using a fixed string. --- test/PolynyaIntegration_test.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/PolynyaIntegration_test.py b/test/PolynyaIntegration_test.py index 889994c2d..0f7e2beb1 100644 --- a/test/PolynyaIntegration_test.py +++ b/test/PolynyaIntegration_test.py @@ -102,7 +102,7 @@ def __make_init_column(cls): nfirst = y // res nsecond = x // res - fname = f"init_polynya.nc" + fname = cls.init_file # The model expects everything in metres initializer = initMaker(fname, nfirst, nsecond, res*1e3, isWinton=True, checkZeros=False) From 11ffd45d66db0e716f408179d13fd94c49cc9814 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Einar=20=C3=96rn=20=C3=93lason?= Date: Mon, 21 Jul 2025 09:23:20 +0200 Subject: [PATCH 31/33] Document the PolynyaIntegration_test.py with a docstring As suggested by Joe on github --- test/PolynyaIntegration_test.py | 7 +++++++ 1 file changed, 7 insertions(+) diff --git a/test/PolynyaIntegration_test.py b/test/PolynyaIntegration_test.py index 0f7e2beb1..d3a60b508 100644 --- a/test/PolynyaIntegration_test.py +++ b/test/PolynyaIntegration_test.py @@ -8,6 +8,13 @@ class Polynya(unittest.TestCase): + """ + A test class based on the Bjornsson et al. (2001) polynya case. We run the model for 5 days and check the output + against known values. The main purpose of this test is to check that the model is running and that the output is + reasonable. + See Bjornsson et al. (2001) doi:10.3402/tellusa.v53i2.12184 for details. + """ + # A few useful global variables for the class executable = "../nextsim" From 8c5962bc7c80e234e4df3e2679f26320291057cf Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Einar=20=C3=96rn=20=C3=93lason?= Date: Mon, 28 Jul 2025 15:07:03 +0200 Subject: [PATCH 32/33] Update PolynyaIntegration_test.py "known good values" This test seems almost too sensitive to small changes in the model. Changing the advection scheme caused minor changes in the results; they aren't physically important, but the mean thickness and concentration changed a bit. I don't think this invalidates the test. The mean concentration also changes by 0.001 between platforms. --- test/PolynyaIntegration_test.py | 9 +++++++-- 1 file changed, 7 insertions(+), 2 deletions(-) diff --git a/test/PolynyaIntegration_test.py b/test/PolynyaIntegration_test.py index d3a60b508..d3f697e53 100644 --- a/test/PolynyaIntegration_test.py +++ b/test/PolynyaIntegration_test.py @@ -12,6 +12,11 @@ class Polynya(unittest.TestCase): A test class based on the Bjornsson et al. (2001) polynya case. We run the model for 5 days and check the output against known values. The main purpose of this test is to check that the model is running and that the output is reasonable. + + NB! This test is very sensitive to small perturbations in the model. If the test fails, run the full-resolution + version, cited below, and compare the results. You may need to adjust the "known good values" of the test + accordingly. + See Bjornsson et al. (2001) doi:10.3402/tellusa.v53i2.12184 for details. """ @@ -167,7 +172,7 @@ def test_iceThickness(self): The maximum varies by (at least) 3 cm between platforms, making a test of it useless """ - mean = 0.1317 + mean = 0.1323 min = 0.0000 hice = self.hice[:,:,:,0] self.assertAlmostEqual(min, hice.min(), 4, "Min ice thickness not ~= " + str(min) + " m") @@ -185,7 +190,7 @@ def test_concentration(self): cice = self.cice[:,:,:,0] self.assertAlmostEqual(max, cice.max(), 4, "Max conentration not ~= " + str(max)) self.assertAlmostEqual(min, cice.min(), 4, "Min concentration not ~= " + str(min)) - self.assertAlmostEqual(mean, cice.mean(), 3, "Mean concentration not ~= " + str(mean)) + self.assertAlmostEqual(mean, cice.mean(), 2, "Mean concentration not ~= " + str(mean)) if __name__ == '__main__': unittest.main() From f76e2cfe860d0b87e4ccd3fefdb0fd2c1ae7855f Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Einar=20=C3=96rn=20=C3=93lason?= Date: Tue, 19 Aug 2025 08:04:31 +0200 Subject: [PATCH 33/33] Update the polynya integration test script Several changes in develop require coding changes in the script. In particular, there are no netCDF groups in the outputs any more, and temperature is not explicitly initialised when using initMaker. --- test/PolynyaIntegration_test.py | 13 +++++-------- 1 file changed, 5 insertions(+), 8 deletions(-) diff --git a/test/PolynyaIntegration_test.py b/test/PolynyaIntegration_test.py index d3f697e53..e5b1b1db4 100644 --- a/test/PolynyaIntegration_test.py +++ b/test/PolynyaIntegration_test.py @@ -53,10 +53,10 @@ def setUpClass(cls): # Load the basic variables root = netCDF4.Dataset(cls.diagnostics_file, "r", format="NETCDF4") - cls.cice = np.squeeze(np.array(root.groups["data"].variables["cice"][:].data)) - cls.hice = np.squeeze(np.array(root.groups["data"].variables["hice"][:].data)) - cls.uice = np.array(root.groups["data"].variables["u"][:].data) - cls.vice = np.array(root.groups["data"].variables["v"][:].data) + cls.cice = np.squeeze(np.array(root.variables["cice"][:].data)) + cls.hice = np.squeeze(np.array(root.variables["hice"][:].data)) + cls.uice = np.array(root.variables["u"][:].data) + cls.vice = np.array(root.variables["v"][:].data) @classmethod def __make_cfg_file(cls): @@ -117,7 +117,7 @@ def __make_init_column(cls): fname = cls.init_file # The model expects everything in metres - initializer = initMaker(fname, nfirst, nsecond, res*1e3, isWinton=True, checkZeros=False) + initializer = initMaker(fname, nfirst, nsecond, res*1e3, checkZeros=False) # Ice everywhere and all boundaries closed, except the x = 100 km end initializer.mask[:, :] = 1. @@ -143,9 +143,6 @@ def __make_init_column(cls): initializer.sss[:, :] = ocean_salinity initializer.sst[:, :] = ocean_temperature - initializer.tsurf[:, :] = ice_salinity * mu - initializer.tbott = initializer.tsurf - initializer.tintr = initializer.tsurf # All other variables are zero or not needed