GitHub PullRequestReviewCommentEvent at 2014-05-19 22:41:28
Repository: aspect
Actor: bangerth
Payload: {"comment":{"url":"https://api.github.com/repos/geodynamics/aspect/pulls/comments/12820333","id":12820333,"diff_hunk":"@@ -0,0 +1,252 @@\n+/*\n+  Copyright (C) 2014 by the authors of the ASPECT code.\n+\n+  This file is part of ASPECT.\n+\n+  ASPECT is free software; you can redistribute it and/or modify\n+  it under the terms of the GNU General Public License as published by\n+  the Free Software Foundation; either version 2, or (at your option)\n+  any later version.\n+\n+  ASPECT is distributed in the hope that it will be useful,\n+  but WITHOUT ANY WARRANTY; without even the implied warranty of\n+  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the\n+  GNU General Public License for more details.\n+\n+  You should have received a copy of the GNU General Public License\n+  along with ASPECT; see the file doc/COPYING.  If not see\n+  <http://www.gnu.org/licenses/>.\n+*/\n+/*  $Id  */\n+\n+//#define _USE_MATH_DEFINES\n+#include <cmath>\n+#include <aspect/initial_conditions/solidus.h>\n+#include <aspect/geometry_model/spherical_shell.h>\n+#include <aspect/boundary_temperature/interface.h>\n+#include <boost/math/special_functions/spherical_harmonic.hpp>\n+\n+namespace aspect\n+{\n+    namespace melting\n+    {\n+        Melting_curve::Melting_curve(const std::string &filename)\n+        {\n+            read(filename);\n+        }\n+        void Melting_curve::read(const std::string &filename)\n+        {\n+            std::ifstream in(filename.c_str(), std::ios::in);\n+            char temp[256];\n+            std::string T_Unit,P_Unit;\n+          \t        Num_points=0;\n+            if(in.fail())return;\n+            in.getline(temp,256);\n+            in>>T_Unit>>P_Unit;\n+            in.getline(temp,256);\n+          \twhile(!in.eof())\n+          \t{\n+                double T,p;\n+                in>>T>>p;\n+                if(!in.fail())\n+                {\n+                    //Unit switching\n+                    if(T_Unit==\"C\")    T+=273.15; // Degree C to K\n+                    if(P_Unit==\"kbar\") p*=1.e8;   // kbar to Pa\n+                    if(P_Unit==\"GPa\")  p*=1.e9;   // GPa to Pa\n+                    is_radius=false;              // Second column is pressure\n+                    if(P_Unit==\"km\")\n+                    {\n+                        is_radius=true;           // Second column in radius instead of pressure\n+                        p*=1.e3;                  // km to meters\n+                    }\n+                    T_array.push_back(T);\n+                    P_array.push_back(p);\n+                    Num_points++;\n+              \t}\n+                in.getline(temp,256);\n+          \t}\n+      \t}\n+        \n+        double Melting_curve::T(const double p, const double radius) const\n+        {\n+            double T_value,P_value=is_radius?radius:p;\n+            if(T_array.size()==0)return(0.);\n+            for(unsigned i=1;i<Num_points;i++)\n+            {\n+                if(     (i==Num_points-1) || \n+                        (is_radius && P_value>P_array[i]) ||\n+                        (!is_radius && P_value<P_array[i]) )\n+                {\n+                    T_value=T_array[i-1]+(T_array[i]-T_array[i-1])/(P_array[i]-P_array[i-1])*(P_value-P_array[i-1]);\n+                    break;\n+                }\n+            }\n+            return(T_value);\n+        }\n+        \n+    }\n+\n+  namespace InitialConditions\n+  {\n+    template <int dim>\n+    Solidus<dim>::Solidus ()\n+    {}\n+\n+    template <int dim>\n+    double\n+    Solidus<dim>::\n+    initial_temperature (const Point<dim> &position) const\n+    {\n+        double T_min,T_litho;\n+        double T_solidus,T_perturbation;\n+        double litho_thick_theta;\n+\t\tdouble lateral_perturbation;\n+        double Depth=this->geometry_model->depth(position);\n+        static melting::Melting_curve Solidus_curve(solidus_filename);\n+\n+        AssertThrow(Solidus_curve.is_radius==true,ExcMessage(\"The solidus curve has to be depth dependent.\"));\n+        AssertThrow(Solidus_curve.Num_points!=0,ExcMessage(\"Error eading solidus file.\"));\n+        const GeometryModel::SphericalShell<dim> *spherical_geometry_model=\n+\t\t\tdynamic_cast< const GeometryModel::SphericalShell<dim> *>(this->geometry_model);\n+\t\tAssertThrow(spherical_geometry_model!=0,\n+                ExcMessage(\"This initial condition can only be work with spherical shell geometry model.\"));\n+        T_min=(this->get_boundary_temperature()).minimal_temperature();\n+\n+        // In case of spherical shell calculate spherical coordinates\n+\t\tconst Tensor<1,dim> scoord = spherical_surface_coordinates(position);\n+\t\tif(dim==2)\n+\t\t{\n+           // Use a sine as lateral perturbation that is scaled to the opening angle of the geometry.\n+           // This way the perturbation is alway 0 at the model boundaries.\n+           const double opening_angle = spherical_geometry_model->opening_angle()*numbers::PI/180.0;\n+           lateral_perturbation = std::sin(lateral_wave_number_1*scoord[1]*numbers::PI/opening_angle);\t\n+\t\t}\n+\t\telse if(dim==3)\n+\t\t{\n+\t\t\t// Spherical harmonics are only defined for order <= degree\n+            // and degree >= 0. Verify that it is indeed.\n+\t\t\tAssert ( std::abs(lateral_wave_number_2) <= lateral_wave_number_1,\n+                       ExcMessage (\"Spherical harmonics can only be computed for \"\n+                                   \"order <= degree.\"));\n+            Assert ( lateral_wave_number_1 >= 0,\n+                       ExcMessage (\"Spherical harmonics can only be computed for \"\n+                                   \"degree >= 0.\"));\n+            // use a spherical harmonic function as lateral perturbation\n+            lateral_perturbation = boost::math::spherical_harmonic_r(lateral_wave_number_1,lateral_wave_number_2,scoord[2],scoord[1]);\n+\t\t}\n+        litho_thick_theta=litho_thick-Magnitude_lith*lateral_perturbation;\n+        T_litho=Solidus_curve.T(0,spherical_geometry_model->R1-litho_thick_theta)+deltaT;\n+        \n+        if(litho_thick_theta>0 && Depth<litho_thick_theta)\n+            T_solidus=T_min+(T_litho-T_min)*(Depth/litho_thick_theta);\n+        else\n+            T_solidus=Solidus_curve.T(0,sqrt(position.square()))+deltaT;\n+\n+        T_perturbation=Depth/( this->geometry_model->maximal_depth() )*Magnitude_T*lateral_perturbation;\n+        return T_solidus+T_perturbation;\n+    }\n+    template <int dim>\n+    const Tensor<1,dim>\n+    Solidus<dim>::\n+    spherical_surface_coordinates(const Tensor<1,dim> &position) const\n+    {\n+      Tensor<1,dim> scoord;\n+\n+      scoord[0] = std::sqrt(position.norm_square()); // R\n+      scoord[1] = std::atan2(position[1],position[0]); // Phi\n+      if (scoord[1] < 0.0) scoord[1] = 2*numbers::PI + scoord[1]; // correct phi to [0,2*pi]\n+      if (dim==3)\n+        scoord[2] = std::acos(position[2]/std::sqrt(position.norm_square())); // Theta\n+\n+      return scoord;\n+    }\n+\n+    template <int dim>\n+    void\n+    Solidus<dim>::declare_parameters (ParameterHandler &prm)\n+    {\n+      prm.enter_subsection(\"Initial conditions\");\n+      {\n+        prm.declare_entry (\"Supersolidus\",\"0e0\",\n+                           Patterns::Double (),\n+                           \"The difference from solidus.\");\n+        prm.declare_entry (\"Lithosphere thickness\",\"0\",\n+                           Patterns::Double (0),\n+                           \"The thickness of lithosphere thickness. Unit: m\");        \n+        prm.enter_subsection(\"Perturbation\");\n+        {\n+            prm.declare_entry (\"Temperature amplitude\", \"0e0\",\n+                               Patterns::Double (0),\n+                               \"The amplitude of the initial spherical temperature perturbation in (K)\");\n+            prm.declare_entry (\"Lithosphere thickness amplitude\", \"0e0\",\n+                               Patterns::Double (),\n+                               \"The amplitude of the initial lithosphere thickness perturbation in (m)\");\n+            prm.declare_entry (\"Lateral wave number one\",\"3\",\n+                               Patterns::Integer(),\n+                               \"Doubled first lateral wave number of the harmonic perturbation. \"\n+                               \"Equals the spherical harmonic degree in 3D spherical shells. \"\n+                               \"In all other cases one equals half of a sine period over \"\n+                               \"the model domain. This allows for single up-/downswings. \"\n+                               \"Negative numbers reverse the sign of the perturbation but are \"\n+                               \"not allowed for the spherical harmonic case.\");\n+            prm.declare_entry (\"Lateral wave number two\", \"2\",\n+                               Patterns::Integer (),\n+                               \"Doubled second lateral wave number of the harmonic perturbation. \"\n+                               \"Equals the spherical harmonic order in 3D spherical shells. \"\n+                               \"In all other cases one equals half of a sine period over \"\n+                               \"the model domain. This allows for single up-/downswings. \"\n+                               \"Negative numbers reverse the sign of the perturbation.\");\n+        }\n+        prm.leave_subsection();\n+        prm.enter_subsection (\"Data\");\n+        {\n+            prm.declare_entry (\"Solidus filename\", \"\",\n+                               Patterns::Anything(),\n+                               \"The solidus filename.\");\n+      \t}\n+        prm.leave_subsection();\n+      }\n+      prm.leave_subsection();\n+    }\n+\n+\n+    template <int dim>\n+    void\n+    Solidus<dim>::parse_parameters (ParameterHandler &prm)\n+    {\n+      prm.enter_subsection(\"Initial conditions\");\n+      {\n+        deltaT=prm.get_double(\"Supersolidus\");\n+        litho_thick=prm.get_double(\"Lithosphere thickness\");\n+        prm.enter_subsection(\"Perturbation\");\n+        {\n+            Magnitude_T    = prm.get_double(\"Temperature amplitude\");\n+            Magnitude_lith = prm.get_double(\"Lithosphere thickness amplitude\");\n+            lateral_wave_number_1 = prm.get_integer (\"Lateral wave number one\");\n+            lateral_wave_number_2 = prm.get_integer (\"Lateral wave number two\");\n+        }\n+        prm.leave_subsection();\n+        prm.enter_subsection(\"Data\");\n+        {\n+            solidus_filename=prm.get (\"Solidus filename\");\n+        }\n+        prm.leave_subsection();\n+      }\n+      prm.leave_subsection();\n+    }\n+\n+  }\n+}\n+\n+// explicit instantiations\n+namespace aspect\n+{\n+  namespace InitialConditions\n+  {\n+    ASPECT_REGISTER_INITIAL_CONDITIONS(Solidus,\n+                                       \"solidus\",\n+                                       \"Temperature initial condition as solidus,\" \n+                                       \"with perturbation as sin() funciton.\");\n+  }","path":"source/initial_conditions/solidus.cc","position":251,"original_position":251,"commit_id":"5e237669b8542a072dafe51d93622465ae42acfd","original_commit_id":"5e237669b8542a072dafe51d93622465ae42acfd","user":{"login":"bangerth","id":7504421,"avatar_url":"https://avatars.githubusercontent.com/u/7504421?","gravatar_id":"bb0cb3b804901b1b3cb33cc459dd5cab","url":"https://api.github.com/users/bangerth","html_url":"https://github.com/bangerth","followers_url":"https://api.github.com/users/bangerth/followers","following_url":"https://api.github.com/users/bangerth/following{/other_user}","gists_url":"https://api.github.com/users/bangerth/gists{/gist_id}","starred_url":"https://api.github.com/users/bangerth/starred{/owner}{/repo}","subscriptions_url":"https://api.github.com/users/bangerth/subscriptions","organizations_url":"https://api.github.com/users/bangerth/orgs","repos_url":"https://api.github.com/users/bangerth/repos","events_url":"https://api.github.com/users/bangerth/events{/privacy}","received_events_url":"https://api.github.com/users/bangerth/received_events","type":"User","site_admin":false},"body":"This should have a significantly expanded description. Can you write a paragraph that explains what this class aims to do, where it gets its data from, etc.?","created_at":"2014-05-19T22:41:28Z","updated_at":"2014-05-19T22:41:28Z","html_url":"https://github.com/geodynamics/aspect/pull/19#discussion_r12820333","pull_request_url":"https://api.github.com/repos/geodynamics/aspect/pulls/19","_links":{"self":{"href":"https://api.github.com/repos/geodynamics/aspect/pulls/comments/12820333"},"html":{"href":"https://github.com/geodynamics/aspect/pull/19#discussion_r12820333"},"pull_request":{"href":"https://api.github.com/repos/geodynamics/aspect/pulls/19"}}}}