Implemented continental geotherm from Chapman (1986)

include/world_builder/features/continental_plate_models/temperature/chapman.h

@@ -0,0 +1,106 @@
+/*
+  Copyright (C) 2018-2024 by the authors of the World Builder code.
+
+  This file is part of the World Builder.
+
+   This program is free software: you can redistribute it and/or modify
+   it under the terms of the GNU Lesser General Public License as published
+   by the Free Software Foundation, either version 2 of the License, or
+   (at your option) any later version.
+
+   This program is distributed in the hope that it will be useful,
+   but WITHOUT ANY WARRANTY; without even the implied warranty of
+   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+   GNU Lesser General Public License for more details.
+
+   You should have received a copy of the GNU Lesser General Public License
+   along with this program.  If not, see <https://www.gnu.org/licenses/>.
+*/
+
+/*
+  Typical continental geotherm,
+  from D.S. Chapman (1986), "Thermal gradients in the continental crust"
+  Geological Society of London Special Publications, v.24, p.63-70.
+*/
+
+#ifndef WORLD_BUILDER_FEATURES_CONTINENTAL_PLATE_MODELS_TEMPERATURE_CHAPMAN_H
+#define WORLD_BUILDER_FEATURES_CONTINENTAL_PLATE_MODELS_TEMPERATURE_CHAPMAN_H
+
+
+#include "world_builder/features/continental_plate_models/temperature/interface.h"
+#include "world_builder/features/feature_utilities.h"
+#include "world_builder/objects/surface.h"
+
+
+namespace WorldBuilder
+{
+
+  namespace Features
+  {
+    using namespace FeatureUtilities;
+    namespace ContinentalPlateModels
+    {
+      namespace Temperature
+      {
+        /**
+         * This class represents a continental plate and can implement submodules
+         * for temperature and composition. These submodules determine what
+         * the returned temperature or composition of the temperature and composition
+         * functions of this class will be.
+         */
+        class Chapman final: public Interface
+        {
+          public:
+            /**
+             * constructor
+             */
+            Chapman(WorldBuilder::World *world);
+
+            /**
+             * Destructor
+             */
+            ~Chapman() override final;
+
+            /**
+             * declare and read in the world builder file into the parameters class
+             */
+            static
+            void declare_entries(Parameters &prm, const std::string &parent_name = "");
+
+            /**
+             * declare and read in the world builder file into the parameters class
+             */
+            void parse_entries(Parameters &prm, const std::vector<Point<2>> &coordinates) override final;
+
+
+            /**
+             * Returns a temperature based on the given position, depth in the model,
+             * gravity and current temperature.
+             */
+            double get_temperature(const Point<3> &position,
+                                   const Objects::NaturalCoordinate &position_in_natural_coordinates,
+                                   const double depth,
+                                   const double gravity,
+                                   double temperature,
+                                   const double feature_min_depth,
+                                   const double feature_max_depth) const override final;
+
+
+          private:
+            double top_temperature;
+            double top_heat_flux;
+            double thermal_conductivity;
+            double heat_production_per_unit_volume;
+            double min_depth;
+            Objects::Surface min_depth_surface;
+            double max_depth;
+            Objects::Surface max_depth_surface;
+            Operations operation;
+
+        };
+      } // namespace Temperature
+    } // namespace ContinentalPlateModels
+  } // namespace Features
+} // namespace WorldBuilder
+
+#endif

source/world_builder/features/continental_plate_models/temperature/chapman.cc

@@ -0,0 +1,159 @@
+/*
+  Copyright (C) 2018-2024 by the authors of the World Builder code.
+
+  This file is part of the World Builder.
+
+   This program is free software: you can redistribute it and/or modify
+   it under the terms of the GNU Lesser General Public License as published
+   by the Free Software Foundation, either version 2 of the License, or
+   (at your option) any later version.
+
+   This program is distributed in the hope that it will be useful,
+   but WITHOUT ANY WARRANTY; without even the implied warranty of
+   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+   GNU Lesser General Public License for more details.
+
+   You should have received a copy of the GNU Lesser General Public License
+   along with this program.  If not, see <https://www.gnu.org/licenses/>.
+*/
+
+#include "world_builder/features/continental_plate_models/temperature/chapman.h"
+
+
+#include "world_builder/nan.h"
+#include "world_builder/types/array.h"
+#include "world_builder/types/double.h"
+#include "world_builder/types/object.h"
+#include "world_builder/types/one_of.h"
+#include "world_builder/types/value_at_points.h"
+#include "world_builder/world.h"
+
+
+namespace WorldBuilder
+{
+  using namespace Utilities;
+
+  namespace Features
+  {
+    using namespace FeatureUtilities;
+    namespace ContinentalPlateModels
+    {
+      namespace Temperature
+      {
+        Chapman::Chapman(WorldBuilder::World *world_)
+          :
+          top_temperature(NaN::DSNAN),
+          top_heat_flux(NaN::DSNAN),
+          thermal_conductivity(NaN::DSNAN),
+          heat_production_per_unit_volume(NaN::DSNAN),
+          min_depth(NaN::DSNAN),
+          max_depth(NaN::DSNAN),
+          operation(Operations::REPLACE)
+        {
+          this->world = world_;
+          this->name = "chapman";
+        }
+
+        Chapman::~Chapman()
+          = default;
+
+        void
+        Chapman::declare_entries(Parameters &prm, const std::string & /*unused*/)
+        {
+          // Declare entries of this plugin
+          prm.declare_entry("top temperature", Types::Double(293.15),
+                            "The temperature at the top surface in K of this feature."
+                            "If the value is below zero, then an adiabatic temperature is used.");
+
+          prm.declare_entry("top heat flux", Types::Double(0.055)),
+                            "The heat flux at the top surface in W m^(-2) of this feature."
+                            "The default value is 0.055.");
+
+          prm.declare_entry("thermal conductivity", Types::Double(2.5)),
+          "The thermal conductivity in W m^(-1) K^(-1) of this feature."
+          "The default value is 2.5.");
+
+          prm.declare_entry("heat generation per unit volume", Types::Double(1.e-6)),
+          "The heat generation per unit volume in W m^(-3) of this feature."
+          "The default value is 1e-6.");
+
+          prm.declare_entry("min depth", Types::OneOf(Types::Double(0),Types::Array(Types::ValueAtPoints(0., 2.))),
+                            "The depth in m from which the composition of this feature is present.");
+
+          prm.declare_entry("max depth", Types::OneOf(Types::Double(std::numeric_limits<double>::max()),Types::Array(Types::ValueAtPoints(std::numeric_limits<double>::max(), 2.))),
+                            "The depth in m to which the composition of this feature is present.");
+
+        }
+
+        void
+        Chapman::parse_entries(Parameters &prm, const std::vector<Point<2>> &coordinates)
+        {
+          min_depth_surface = Objects::Surface(prm.get("min depth",coordinates));
+          min_depth = min_depth_surface.minimum;
+          max_depth_surface = Objects::Surface(prm.get("max depth",coordinates));
+          max_depth = max_depth_surface.maximum;
+          WBAssert(max_depth >= min_depth, "max depth needs to be larger or equal to min depth.");
+          operation = string_operations_to_enum(prm.get<std::string>("operation"));
+          thermal_conductivity = prm.get<double>("thermal conductivity");
+          heat_production_per_unit_volume = prm.get<double>("heat generation per unit volume");
+          top_heat_flux = prm.get<double>("top heat flux");
+          top_temperature = prm.get<double>("top temperature");
+          WBAssert(!std::isnan(top_temperature), "Top surface temperature is not a number: " << top_temperature
+                   << ", based on a temperature model with the name " << this->name);
+        }
+
+
+        double
+        Chapman::get_temperature(const Point<3> & /*position_in_cartesian_coordinates*/,
+                                 const Objects::NaturalCoordinate &position_in_natural_coordinates,
+                                 const double depth,
+                                 const double gravity_norm,
+                                 double temperature_,
+                                 const double feature_min_depth,
+                                 const double feature_max_depth) const
+        {
+          if (depth <= max_depth && depth >= min_depth)
+            {
+              const double min_depth_local = min_depth_surface.constant_value ? min_depth : min_depth_surface.local_value(position_in_natural_coordinates.get_surface_point()).interpolated_value;
+              const double max_depth_local = max_depth_surface.constant_value ? max_depth : max_depth_surface.local_value(position_in_natural_coordinates.get_surface_point()).interpolated_value;
+              if (depth <= max_depth_local && depth >= min_depth_local)
+                {
+                  const double min_depth_local_local = std::max(feature_min_depth, min_depth_local);
+                  const double max_depth_local_local = std::min(feature_max_depth, max_depth_local);
+
+                  double top_temperature_local = top_temperature;
+
+                  // use adiabatic temperature if top temperature is below zero
+                  if (top_temperature_local < 0)
+                    {
+                      top_temperature_local =  this->world->potential_mantle_temperature *
+                                               std::exp(((this->world->thermal_expansion_coefficient * gravity_norm) /
+                                                         this->world->specific_heat) * min_depth_local_local);
+                    }
+
+                  // T(z) = t_top + (q_top / k) * (Δz) - (A / (2 * k)) * (Δz)^2
+                  const double new_temperature = top_temperature + (top_heat_flux / thermal_conductivity) * (depth-min_depth_local_local) - heat_production_per_unit_volume / (2. * thermal_conductivity) * (depth-min_depth_local_local)*(depth-min_depth_local_local);
+
+                  WBAssert(!std::isnan(new_temperature), "Temperature is not a number: " << new_temperature
+                           << ", based on a temperature model with the name " << this->name);
+                  WBAssert(std::isfinite(new_temperature), "Temperature is not a finite: " << new_temperature
+                           << ", based on a temperature model with the name " << this->name
+                           << ", top_temperature_local = " << top_temperature_local << ", depth = " << depth << ", min_depth_local = " << min_depth_local
+                           << ", top_heat_flux = " << top_heat_flux << ", thermal_conductivity=" << thermal_conductivity
+                           << ",max_depth_local_local=" << max_depth_local_local << ", min_depth_local_local =" << min_depth_local_local
+                           << ", feature_max_depth = " << feature_max_depth << ", feature_max_depth = " << feature_max_depth);
+
+                  return apply_operation(operation,temperature_,new_temperature);
+                }
+            }
+
+
+          return temperature_;
+        }
+
+        WB_REGISTER_FEATURE_CONTINENTAL_PLATE_TEMPERATURE_MODEL(Chapman, Chapman)
+      } // namespace Temperature
+    } // namespace ContinentalPlateModels
+  } // namespace Features
+} // namespace WorldBuilder
+