add temperature plate model constant age to oceanic plate models

CHANGELOG.md

@@ -15,6 +15,7 @@ This project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.htm
 - The option to use a half space cooling model for the temperature of an oceanic plate has been added to the oceanic plate feature. \[Magali Billen; 2021-07-15; [#317](github.com/GeodynamicWorldBuilder/WorldBuilder/pull/317)\] 
 - Information about the depth of the reference surface in the PointDistanceFromCurvedPlanes structure. \[Menno Fraters; 2021-07-20; [#320](github.com/GeodynamicWorldBuilder/WorldBuilder/pull/320)\]
 -The option for a smooth slab temperature structure that conserves mass with distance along the slab surface. \[Magali Billen; 2021-07-28; [#323] (github.com/GeodynamicWorldBuilder/WorldBuilder/pull/323)\]
+-The option for a temperature model for an oceanic plate with a constant age using the plate model. \[Haoyuan Li; 2021-10-29; [#344] (https://github.com/GeodynamicWorldBuilder/WorldBuilder/pull/344)\]
 ### Changed
 - The World Builder Visualizer will now use zlib compression for vtu files by default. If zlib is not available binary output will be used. \[Menno Fraters; 2021-06-26; [#282](github.com/GeodynamicWorldBuilder/WorldBuilder/pull/282)\]
 - The return argument type of the distance_point_from_curved_planes function has been converted from a map to a struct, requiring a change in the plugin interfaces for 'fault_models' and 'subducting_plate_models', but significantly increasing the speed of the function and all functions that access its returned values. \[Rene Gassmoeller; 2021-06-27; [#289](github.com/GeodynamicWorldBuilder/WorldBuilder/issues/289)\]

include/world_builder/features/oceanic_plate_models/temperature/plate_model_constant_age.h

@@ -0,0 +1,100 @@
+/*
+  Copyright (C) 2018 - 2021 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/>.
+*/
+
+#ifndef WORLD_BUILDER_FEATURES_OCEANIC_PLATE_MODELS_TEMPERATURE_PLATE_MODEL_CONSTANT_AGE_H
+#define WORLD_BUILDER_FEATURES_OCEANIC_PLATE_MODELS_TEMPERATURE_PLATE_MODEL_CONSTANT_AGE_H
+
+
+#include "world_builder/features/oceanic_plate_models/temperature/interface.h"
+#include "world_builder/features/utilities.h"
+
+
+namespace WorldBuilder
+{
+  class Parameters;
+  class World;
+
+  namespace Features
+  {
+    namespace OceanicPlateModels
+    {
+      namespace Temperature
+      {
+        /**
+         * This class represents a oceanic plate of constant age 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.
+         * In this plugin, the temperature of the plate is derived from the plate model,
+         * featuring cooling of a plate from a ridge by vertical heat conduction,
+         * while the effects of horizontal heat conduction are ignored.
+         */
+        class PlateModelConstantAge final: public Interface
+        {
+          public:
+            /**
+             * constructor
+             */
+            PlateModelConstantAge(WorldBuilder::World *world);
+
+            /**
+             * Destructor
+             */
+            ~PlateModelConstantAge() 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) 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 double depth,
+                                   const double gravity,
+                                   double temperature,
+                                   const double feature_min_depth,
+                                   const double feature_max_depth) const override final;
+
+
+          private:
+            // plate model constant age temperature submodule parameters
+            double min_depth;
+            double max_depth;
+            double top_temperature;
+            double bottom_temperature;
+            double plate_age;
+            Utilities::Operations operation;
+
+        };
+      } // namespace Temperature
+    } // namespace OceanicPlateModels
+  } // namespace Features
+} // namespace WorldBuilder
+
+#endif

source/features/oceanic_plate_models/temperature/plate_model_constant_age.cc

@@ -0,0 +1,151 @@
+/*
+  Copyright (C) 2018 - 2021 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/oceanic_plate_models/temperature/plate_model_constant_age.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/point.h"
+#include "world_builder/utilities.h"
+#include "world_builder/world.h"
+
+
+namespace WorldBuilder
+{
+  using namespace Utilities;
+
+  namespace Features
+  {
+    namespace OceanicPlateModels
+    {
+      namespace Temperature
+      {
+        PlateModelConstantAge::PlateModelConstantAge(WorldBuilder::World *world_)
+          :
+          min_depth(NaN::DSNAN),
+          max_depth(NaN::DSNAN),
+          top_temperature(NaN::DSNAN),
+          bottom_temperature(NaN::DSNAN),
+          operation(Utilities::Operations::REPLACE)
+        {
+          this->world = world_;
+          this->name = "plate model constant age";
+        }
+
+        PlateModelConstantAge::~PlateModelConstantAge()
+          = default;
+
+        void
+        PlateModelConstantAge::declare_entries(Parameters &prm, const std::string & /*unused*/)
+        {
+
+          // Add temperature to the required parameters.
+          prm.declare_entry("", Types::Object({"max depth"}), "Temperature model object");
+
+          prm.declare_entry("min depth", Types::Double(0),
+                            "The depth in meters from which the temperature of this feature is present.");
+
+          prm.declare_entry("max depth", Types::Double(std::numeric_limits<double>::max()),
+                            "The depth in meters to which the temperature of this feature is present.");
+
+          prm.declare_entry("top temperature", Types::Double(293.15),
+                            "The temperature in degree Kelvin which this feature should have");
+
+          prm.declare_entry("bottom temperature", Types::Double(-1),
+                            "The temperature in degree Kelvin which this feature should have");
+
+          prm.declare_entry("plate age", Types::Double(80e3),
+                            "The age of the plate in year. "
+                            "This age is assigned to the whole plate. ");
+        }
+
+        void
+        PlateModelConstantAge::parse_entries(Parameters &prm)
+        {
+
+          min_depth = prm.get<double>("min depth");
+          max_depth = prm.get<double>("max depth");
+          operation = Utilities::string_operations_to_enum(prm.get<std::string>("operation"));
+          top_temperature = prm.get<double>("top temperature");
+          bottom_temperature = prm.get<double>("bottom temperature");
+          plate_age = prm.get<double>("plate age")*31557600;
+        }
+
+
+        double
+        PlateModelConstantAge::get_temperature(const Point<3> & /*position*/,
+                                               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)
+            {
+              double bottom_temperature_local = bottom_temperature;
+
+              if (bottom_temperature_local < 0)
+                {
+                  bottom_temperature_local =  this->world->potential_mantle_temperature *
+                                              std::exp(((this->world->thermal_expansion_coefficient* gravity_norm) /
+                                                        this->world->specific_heat) * depth);
+                }
+              const int sommation_number = 100;
+
+              // some aliases
+              //const double top_temperature = top_temperature;
+              const double thermal_diffusivity = this->world->thermal_diffusivity;
+              double temperature = top_temperature + (bottom_temperature_local - top_temperature) * (depth / max_depth);
+
+              // This formula ignores the horizontal heat transfer by just having the age of the plate in it.
+              // (Chapter 7 Heat, Fowler M. The solid earth: an introduction to global geophysics[M]. Cambridge University Press, 1990)
+              for (int i = 1; i<sommation_number+1; ++i)
+                {
+                  temperature = temperature + (bottom_temperature_local - top_temperature) *
+                                ((2 / (double(i) * const_pi)) * std::sin((double(i) * const_pi * depth) / max_depth) *
+                                 std::exp(-1.0 * i * i * const_pi * const_pi * thermal_diffusivity * plate_age / (max_depth * max_depth)));
+                }
+
+              WBAssert(!std::isnan(temperature), "Temparture inside plate model constant age is not a number: " << temperature
+                       << ". Relevant variables: bottom_temperature_local = " << bottom_temperature_local
+                       << ", top_temperature = " << top_temperature
+                       << ", max_depth = " << max_depth
+                       << ", thermal_diffusivity = " << thermal_diffusivity
+                       << ", age = " << plate_age << ".");
+              WBAssert(std::isfinite(temperature), "Temparture inside plate model constant age is not a finite: " << temperature                           << ". Relevant variables: bottom_temperature_local = " << bottom_temperature_local
+                       << ", top_temperature = " << top_temperature
+                       << ", thermal_diffusivity = " << thermal_diffusivity
+                       << ", age = " << plate_age << ".");
+
+
+              return Utilities::apply_operation(operation,temperature_,temperature);
+
+            }
+          return temperature_;
+        }
+
+        WB_REGISTER_FEATURE_OCEANIC_PLATE_TEMPERATURE_MODEL(PlateModelConstantAge, plate model constant age)
+      } // namespace Temperature
+    } // namespace OceanicPlateModels
+  } // namespace Features
+} // namespace WorldBuilder
+

tests/app/app_oceanic_plate_constant_age.dat

@@ -0,0 +1,16 @@
+# This is a comment in the data
+# file. 
+# Now define parameters:
+# dim = 2
+# compositions = 0
+# x z d g T 
+3000e3 0 10e3 10
+3000e3 0 50e3 10
+3000e3 0 80e3 10
+3000e3 0 100e3 10
+3000e3 0 150e3 10
+3000e3 0 250e3 10
+3000e3 0 300e3 10
+4500e3 0 10e3 10
+4500e3 0 50e3 10
+4500e3 0 80e3 10

tests/app/app_oceanic_plate_constant_age.wb

@@ -0,0 +1,15 @@
+{
+  "version":"0.5",
+  "coordinate system":{"model":"cartesian"},
+  "cross section":[[0,0],[1000e3,0]],"surface temperature":273, "force surface temperature":true,
+  "potential mantle temperature":1673, "thermal expansion coefficient":3.1e-5, 
+  "specific heat":1250, "thermal diffusivity":1.0e-6,
+  "features":
+  [
+    { "model":"oceanic plate", "name":"Subducting", "max depth":400e3,"min depth":0,
+      "coordinates" :[[1500e3,100e3],[1500e3,-100e3],[4500e3,-100e3],[4500e3,100e3]],
+      "temperature models":[
+        {"model":"plate model constant age", "min depth":-10e3,  "max depth":250e3, "plate age":40e6,
+        "top temperature":273}]}      
+  ]
+}

tests/app/app_oceanic_plate_constant_age/screen-output.log

@@ -0,0 +1,11 @@
+# x z d g T 
+3000e3 0 10e3 10 494.513 
+3000e3 0 50e3 10 1239.65 
+3000e3 0 80e3 10 1546.91 
+3000e3 0 100e3 10 1647.86 
+3000e3 0 150e3 10 1732.26 
+3000e3 0 250e3 10 1780.01 
+3000e3 0 300e3 10 1802.22 
+4500e3 0 10e3 10 494.513 
+4500e3 0 50e3 10 1239.65 
+4500e3 0 80e3 10 1546.91