Merge pull request #61 from nextsimdg/issue60_physics_module

Issue60 Physics Module

Author: einola

core/src/include/ElementData.hpp

@@ -7,8 +7,11 @@
 #ifndef SRC_INCLUDE_ELEMENTDATA_HPP
 #define SRC_INCLUDE_ELEMENTDATA_HPP
 
+#include "Configured.hpp"
 #include "ExternalData.hpp"
+#include "ModuleLoader.hpp"
 #include "PrognosticData.hpp"
+#include "include/IPhysics1d.hpp"
 #include "include/PhysicsData.hpp"
 
 namespace Nextsim {
@@ -22,31 +25,44 @@ class UnusedData : public BaseElementData {
  * @brief The class which holds all the data for a single element of the model.
  *
  * @details Inherits from PrognosticData, PhysicsData and ExternalData. The
- * physics implementation is provided as a template argument.
+ * physics implementation is provided as a per-instance module.
  */
-template <class Phys>
 class ElementData : public PrognosticData,
                     public PhysicsData,
                     public ExternalData,
-                    public Phys,
                     public UnusedData,
-                    public Configured<ElementData<Phys>> {
+                    public Configured<ElementData> {
 public:
-    ElementData() = default;
+    ElementData()
+    {
+        m_physicsImplData = std::move(ModuleLoader::getLoader().getInstance<IPhysics1d>());
+    }
     ~ElementData() = default;
 
     using PrognosticData::operator=;
     using PhysicsData::operator=;
     using ExternalData::operator=;
-    using Phys::operator=;
 
     //! Configures the PrognosticData and physics implementation aspects of the
     //!  object.
     void configure() override
     {
         PrognosticData::configure();
-        Phys::configure();
+        Nextsim::tryConfigure(&ModuleLoader::getLoader().getImplementation<IPhysics1d>());
+    }
+
+    void updateDerivedData(const PrognosticData& prog, const ExternalData& exter, PhysicsData& phys)
+    {
+        m_physicsImplData->updateDerivedData(prog, exter, phys);
     }
+
+    void calculate(const PrognosticData& prog, const ExternalData& exter, PhysicsData& phys)
+    {
+        m_physicsImplData->calculate(prog, exter, phys);
+    }
+
+private:
+    std::unique_ptr<IPhysics1d> m_physicsImplData;
 };
 
 } /* namespace Nextsim */

core/test/CMakeLists.txt

@@ -88,4 +88,23 @@ add_executable(testScopedTimer
 target_link_libraries(testScopedTimer PRIVATE Catch2::Catch2)
 target_include_directories(testScopedTimer PRIVATE "${SRC_DIR}")
 
+set(PhysicsDir "${PROJECT_SOURCE_DIR}/physics/src")
+set(PhysicsModulesDir "${PhysicsDir}/modules")
+add_executable(testElementData
+    "ElementData_test.cpp"
+    "${SRC_DIR}/ElementData.cpp"
+    "${SRC_DIR}/ModuleLoader.cpp"
+    "${SRC_DIR}/Configurator.cpp"
+    "${SRC_DIR}/ConfiguredModule.cpp"
+    "${SRC_DIR}/PrognosticData.cpp"
+    "${PhysicsModulesDir}/NextsimPhysics.cpp"
+    "${PhysicsModulesDir}/SMUIceAlbedo.cpp"
+    "${PhysicsModulesDir}/CCSMIceAlbedo.cpp"
+    "${PhysicsModulesDir}/SMU2IceAlbedo.cpp"
+    "${PhysicsModulesDir}/BasicIceOceanHeatFlux.cpp"
+    "${PhysicsModulesDir}/HiblerConcentration.cpp"
+    "${PhysicsModulesDir}/ThermoIce0.cpp"
+    )
 
+    target_include_directories(testElementData PRIVATE "${ModuleLoaderIppTargetDirectory}" "${SRC_DIR}" "${PhysicsDir}" "${PhysicsModulesDir}")
+    target_link_libraries(testElementData PRIVATE "${Boost_LIBRARIES}" Catch2::Catch2)

core/test/ElementData_test.cpp

@@ -4,8 +4,77 @@
  * @author Tim Spain <[email protected]>
  */
 
-#include <ElementData.hpp>
+#define CATCH_CONFIG_MAIN
+#include <catch2/catch.hpp>
+#include <sstream>
+
+#include "include/ConfiguredModule.hpp"
+#include "include/ElementData.hpp"
+#include "include/IIceAlbedo.hpp"
+#include "include/ModuleLoader.hpp"
+#include "include/NextsimPhysics.hpp"
 
 namespace Nextsim {
 
+TEST_CASE("Physics test using NextsimPhysics", "[ElementData]")
+{
+    // Largely copied from the "Melting conditions" test of NextsimPhysics
+    Configurator::clear();
+    std::stringstream config;
+    config << "[Modules]" << std::endl;
+    config << "Nextsim::IFreezingPoint = Nextsim::UnescoFreezing" << std::endl;
+    config << "Nextsim::IIceAlbedo = Nextsim::CCSMIceAlbedo" << std::endl;
+    config << "Nextsim::IPhysics1d = Nextsim::NextsimPhysics" << std::endl;
+    config << std::endl;
+    config << "[CCSMIceAlbedo]" << std::endl;
+    config << "iceAlbedo = 0.63" << std::endl;
+    config << "snowAlbedo = 0.88" << std::endl;
+
+    std::unique_ptr<std::istream> pcstream(new std::stringstream(config.str()));
+    Configurator::addStream(std::move(pcstream));
+
+    double tair = 3; //˚C
+    double tdew = 2; //˚C
+    double pair = 100000; // Pa, slightly low pressure
+    double sst = -1; //˚C
+    double sss = 32; // PSU
+    std::array<double, N_ICE_TEMPERATURES> tice = { -1., -1., -1. }; //˚C
+    double hice = 0.1; // m
+    double cice = 0.5;
+    double hsnow = 0.01; // m
+    double dml = 10.; // m
+
+    ModuleLoader::getLoader().setAllDefaults();
+    ConfiguredModule::parseConfigurator();
+    tryConfigure(ModuleLoader::getLoader().getImplementation<IIceAlbedo>());
+
+    ElementData data;
+    data.configure(); // Configure with the UNESCO freezing point
+
+    data = PrognosticData::generate(hice, cice, sst, sss, hsnow, tice);
+    data.setTimestep(600.); // s. Very long TS to get below freezing
+
+    data.airTemperature() = tair;
+    data.dewPoint2m() = tdew;
+    data.airPressure() = pair;
+    data.mixedLayerDepth() = dml;
+    data.incomingLongwave() = 330;
+    data.incomingShortwave() = 50;
+    data.snowfall() = 0;
+
+    data.windSpeed() = 5;
+
+    data.updateDerivedData(data, data, data);
+    REQUIRE(1.26488 == Approx(data.airDensity()).epsilon(1e-4));
+    REQUIRE(0.00441973 == Approx(data.specificHumidityAir()).epsilon(1e-4));
+    REQUIRE(0.0035214 == Approx(data.specificHumidityIce()).epsilon(1e-4));
+    REQUIRE(1012.86 == Approx(data.heatCapacityWetAir()).epsilon(1e-4));
+
+    data.calculate(data, data, data);
+
+    REQUIRE(0.12846 == Approx(data.updatedIceTrueThickness()).epsilon(1e-4));
+    REQUIRE(0.01957732 == Approx(data.updatedSnowTrueThickness()).epsilon(1e-4));
+    REQUIRE(0.368269 == Approx(data.updatedIceConcentration()).epsilon(1e-4));
+    REQUIRE(0.0 == Approx(data.updatedIceSurfaceTemperature()).epsilon(1e-4));
+}
 } /* namespace Nextsim */

physics/src/Physics1dBase.cpp

@@ -8,19 +8,9 @@
 
 namespace Nextsim {
 
-template <class Phys> void Physics1dBase::physics1d(ElementData<Phys>& data)
+void Physics1dBase::physics1d(ElementData& data)
 {
-    Phys::updateDerivedData(data);
-    Phys::massFluxOpenWater(data);
-    Phys::momentumFluxOpenWater(data);
-    Phys::heatFluxOpenWater(data);
-
-    Phys::massFluxIceAtmosphere(data);
-    // Ice momentum fluxes are handled by the dynamics
-    Phys::heatFluxIceAtmosphere(data);
-
-    Phys::massFluxIceOcean(data);
-    // Ice momentum fluxes are handled by the dynamics
-    Phys::heatFluxIceOcean(data);
+    data.updateDerivedData(data, data, data);
+    data.calculate(data, data, data);
 }
 } /* namespace Nextsim */

physics/src/include/Physics1dBase.hpp

@@ -24,7 +24,7 @@ class Physics1dBase {
     virtual ~Physics1dBase() = default;
 
     //! Performs the one dimensional physics calculations.
-    template <class Phys> void physics1d(ElementData<Phys>&);
+    void physics1d(ElementData&);
 };
 
 } /* namespace Nextsim */

physics/src/modules/NextsimPhysics.cpp

@@ -172,12 +172,12 @@ void NextsimPhysics::heatFluxIceAtmosphere(
     m_Qlhi = m_subl * latentHeatIce(prog.iceTemperatures()[0]);
     double dmdot_dT = dragIce_t * phys.airDensity() * phys.windSpeed()
         * specHumIce.dq_dT(prog.iceTemperatures()[0], exter.airPressure());
-    m_dQ_dT += latentHeatIce(prog.iceTemperatures()[0]) * dmdot_dT;
+    double dQlh_dT = latentHeatIce(prog.iceTemperatures()[0]) * dmdot_dT;
 
     // Sensible heat flux
     m_Qshi = dragIce_t * phys.airDensity() * phys.heatCapacityWetAir() * phys.windSpeed()
         * (prog.iceTemperatures()[0] - exter.airTemperature());
-    m_dQ_dT += dragIce_t * phys.airDensity() * phys.heatCapacityWetAir() * phys.windSpeed();
+    double dQsh_dT = dragIce_t * phys.airDensity() * phys.heatCapacityWetAir() * phys.windSpeed();
 
     // Shortwave flux
     double albedoValue = iIceAlbedoImpl->albedo(prog.iceTemperatures()[0],
@@ -186,11 +186,13 @@ void NextsimPhysics::heatFluxIceAtmosphere(
 
     // Longwave flux
     m_Qlwi = stefanBoltzmannLaw(prog.iceTemperatures()[0]) - exter.incomingLongwave();
-    m_dQ_dT
-        += 4 / kelvin(prog.iceTemperatures()[0]) * stefanBoltzmannLaw(prog.iceTemperatures()[0]);
+    double dQlw_dT
+        = 4 / kelvin(prog.iceTemperatures()[0]) * stefanBoltzmannLaw(prog.iceTemperatures()[0]);
 
     // Total flux
     m_Qia = m_Qlhi + m_Qshi + m_Qlwi + m_Qswi;
+    // Overall temperature dependence of flux
+    m_dQ_dT = dQlh_dT + dQsh_dT + dQlw_dT;
 }
 
 void NextsimPhysics::massFluxIceOcean(

physics/src/modules/include/NextsimPhysics.hpp

@@ -25,7 +25,7 @@ class IThermodynamics;
 class IIceOceanHeatFlux;
 class IConcentrationModel;
 
-template <class Phys> class ElementData;
+class ElementData;
 
 class NextsimPhysics;
 

physics/src/modules/modules.json

@@ -31,5 +31,11 @@
         "implementations": [
             "Nextsim::ThermoIce0"
         ]
+    },
+    {
+        "name": "Nextsim::IPhysics1d",
+        "implementations": [
+            "Nextsim::NextsimPhysics"
+        ]
     }
 ]