Adding tidal potential to gravity model

doc/modules/changes/20250716_hyunseong96

@@ -0,0 +1,4 @@
+Added: ASPECT now has a new gravity model plugin called 'Radial linear with tidal potential' for tidal forces.
+This plugin is useful for modeling long-term interior and surface evolution in moons orbiting a large planet.
+<br>
+(Hyunseong Kim, Antoniette Greta Grima, Wolfgang Bangerth 2025/07/16)

include/aspect/gravity_model/radial_with_tidal_potential.h

@@ -0,0 +1,89 @@
+/*
+  Copyright (C) 2014 - 2019 by the authors of the ASPECT code.
+
+  This file is part of ASPECT.
+
+  ASPECT is free software; you can redistribute it and/or modify
+  it under the terms of the GNU General Public License as published by
+  the Free Software Foundation; either version 2, or (at your option)
+  any later version.
+
+  ASPECT 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 General Public License for more details.
+
+  You should have received a copy of the GNU General Public License
+  along with ASPECT; see the file LICENSE.  If not see
+  <http://www.gnu.org/licenses/>.
+*/
+
+
+#ifndef _aspect_gravity_model_radial_with_tidal_potential_h
+#define _aspect_gravity_model_radial_with_tidal_potential_h
+
+#include <aspect/simulator_access.h>
+#include <aspect/gravity_model/interface.h>
+#include <aspect/gravity_model/radial.h>
+
+namespace aspect
+{
+  namespace GravityModel
+  {
+    /**
+     * A class that describes gravity as a radial vector of linearly
+     * changing magnitude, which is modified by a tidal potential from flattening and non-synchronous rotation.
+     *
+     * The equation implemented in this gravity model is from Tobie et al. (2025) (https://doi.org/10.1007/s11214-025-01136-y),
+     * which is defined as:
+     * g = -magnitude - gradient (-(tidal potential)).
+     * Tidal potential is positive because the formula follows conventions from geodesy research, where potential is taken as positive.
+     * (tidal potential) = (3 G M_p) / (2 a_s^3) * r^2 * (Tstar + T0)
+     * Tstar = 1/6 *(1-3*cos(theta)^2) and T0=1/2sin(theta)^2*cos(2*lambda + 2*b*t)
+     * where G = gravitational constant, M_p = mass of the perturbing body, a_s = semimajor axis of the orbit, b = angular rate of non-synchronous rotation.
+     * r, theta and lambda are radial distance, polar angle and azimuthal angle, respectively.
+     *
+     * @ingroup GravityModels
+     */
+    template <int dim>
+    class RadialWithTidalPotential : public Interface<dim>, public SimulatorAccess<dim>
+    {
+      public:
+        /**
+         * Return the gravity vector as a function of position.
+         */
+        Tensor<1,dim> gravity_vector (const Point<dim> &position) const override;
+
+        /**
+         * Declare the parameters this class takes through input files.
+         */
+        static
+        void
+        declare_parameters (ParameterHandler &prm);
+
+        /**
+         * Read the parameters this class declares from the parameter file.
+         */
+        void
+        parse_parameters (ParameterHandler &prm) override;
+
+      private:
+        /**
+         * Mass of the perturbing body
+         */
+        double M_p;
+
+        /**
+         * Semimajor axis of the orbit that causes the tidal perturbation
+         */
+        double a_s;
+
+        /**
+         * Angular rate of the non-synchronous rotation in degrees/year
+         */
+        double b;
+    };
+  }
+}
+
+#endif

source/gravity_model/radial_with_tidal_potential.cc

@@ -0,0 +1,148 @@
+/*
+  Copyright (C) 2011 - 2020 by the authors of the ASPECT code.
+
+  This file is part of ASPECT.
+
+  ASPECT is free software; you can redistribute it and/or modify
+  it under the terms of the GNU General Public License as published by
+  the Free Software Foundation; either version 2, or (at your option)
+  any later version.
+
+  ASPECT 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 General Public License for more details.
+
+  You should have received a copy of the GNU General Public License
+  along with ASPECT; see the file LICENSE.  If not see
+  <http://www.gnu.org/licenses/>.
+*/
+
+
+#include <aspect/gravity_model/radial_with_tidal_potential.h>
+#include <aspect/geometry_model/interface.h>
+#include <aspect/utilities.h>
+
+#include <deal.II/base/tensor.h>
+
+#include <aspect/gravity_model/radial.h>
+
+namespace aspect
+{
+  namespace GravityModel
+  {
+    template <int dim>
+    Tensor<1,dim>
+    RadialWithTidalPotential<dim>::gravity_vector (const Point<dim> &/*p*/) const
+    {
+      // This plugin is not implemented for 2D models
+      AssertThrow(false, ExcNotImplemented());
+      return Tensor<1,dim>();
+    }
+
+    template <>
+    Tensor<1,3>
+    RadialWithTidalPotential<3>::gravity_vector (const Point<3> &p) const
+    {
+      const unsigned int dim = 3;
+      /**
+       * Notation of this potential equation is converted from spherical coordinates to cartesian coordinates.
+       * Therefore, gradient of potential is (3 G M_p) / (2 a_s^3) * ( 1 / 6 * ( x^2 + y^2 - 2 * z^2) + 1 / 2 * (C1*(x^2 + y^2) - 2 * C2 * x * y)))
+       * where C1 = cos(2*b*t) and C2 = sin(2*b*t)
+       */
+      const double t = (this->simulator_is_past_initialization()) ? this->get_time() : 0.0;
+
+      const double C1 = std::cos( 2. * b * t);
+      const double C2 = std::sin( 2. * b * t);
+
+      const Tensor<1,dim> dTstar_gradient ({1./3. * p[0], 1./3. * p[1], -2./3. * p[2]});
+
+      const Tensor<1,dim> dT0_gradient ({C1 *p[0] - C2 *p[1], -C1 *p[1] - C2 *p[0], 0});
+
+      const double G = aspect::constants::big_g;
+      const double T_factor = 3. * G * M_p / ( 2. * a_s * a_s * a_s );
+
+      const Tensor<1,dim> tidal_gravity = T_factor *
+                                          (dTstar_gradient + dT0_gradient);
+
+      RadialConstant<dim> radialconstant;
+      return radialconstant.gravity_vector(p) + tidal_gravity;
+    }
+
+
+    template <int dim>
+    void
+    RadialWithTidalPotential<dim>::declare_parameters (ParameterHandler &prm)
+    {
+      RadialLinear<dim>::declare_parameters(prm);
+      prm.enter_subsection("Gravity model");
+      {
+        prm.enter_subsection("Radial with tidal potential");
+        {
+          prm.declare_entry ("Mass of perturbing body", "1.898e27",
+                             Patterns::Double (),
+                             "Mass of body that perturbs gravity of modeled body. "
+                             "The default value is chosen for modeling Europa, therefore, it is the mass of Jupiter. "
+                             "Units is $kg$.");
+          prm.declare_entry ("Semimajor axis of orbit", "670900000",
+                             Patterns::Double (),
+                             "The length of the semimajor axis of the orbit that cause the tidal perturbation. "
+                             "For example, tidal perturbation on Europa happens by Europa orbiting Jupiter, "
+                             "and that on Earth, if Moon is in consideration, happens by Moon orbiting Earth. "
+                             "The default value is for the semimajor axis of Europa's orbit. "
+                             "Units is $m$.");
+          prm.declare_entry ("Angular rate of nonsynchronous rotation", "0.036",
+                             Patterns::Double (),
+                             "Angular rate of nonsynchronous rotation (NSR). "
+                             "This works for the modeled body having decoupled rotation between interior layers. "
+                             "The default value is the angular rate of Europa's icy shell. "
+                             "Units is $degrees/year$ when 'Use years instead of seconds' is true, "
+                             "and $degress/second$ when 'Use years instead of seconds' is false. ");
+        }
+        prm.leave_subsection ();
+      }
+      prm.leave_subsection ();
+    }
+
+
+    template <int dim>
+    void
+    RadialWithTidalPotential<dim>::parse_parameters (ParameterHandler &prm)
+    {
+      AssertThrow (dim==3, ExcMessage ("The 'radial with tidal potential' gravity model "
+                                       "can only be used in 3D."));
+
+      prm.enter_subsection("Gravity model");
+      {
+        prm.enter_subsection("Radial with tidal potential");
+        {
+          M_p = prm.get_double ("Mass of perturbing body");
+          a_s = prm.get_double ("Semimajor axis of orbit");
+          const double time_scale = this->get_parameters().convert_to_years ? constants::year_in_seconds : 1.0;
+          b = prm.get_double ("Angular rate of nonsynchronous rotation") * constants::degree_to_radians / time_scale;
+        }
+        prm.leave_subsection ();
+      }
+      prm.leave_subsection ();
+    }
+  }
+}
+
+// explicit instantiations
+namespace aspect
+{
+  namespace GravityModel
+  {
+    ASPECT_REGISTER_GRAVITY_MODEL(RadialWithTidalPotential,
+                                  "radial with tidal potential",
+                                  "A gravity model that is the sum of the `radial constant' model "
+                                  "(which is radial, pointing inward if the gravity "
+                                  "is positive), "
+                                  "and a term that results from a tidal potential and that "
+                                  "leads to a gravity field that varies with latitude and longitude."
+                                  "The magnitude of gravity for the radial constant part is read from the "
+                                  "input file in a section `Gravity model/Radial constant'; the "
+                                  "parameters that describe the tidal potential contribution are read "
+                                  "from a section ``Gravity model/Radial with tidal potential''.")
+  }
+}

tests/gravity_tidal_potential.prm

@@ -0,0 +1,99 @@
+# This parameter file tests the gravity model plugin for a case where the
+# tidal potential by flattening and non-synchronnous rotation changes gravity with position and time.
+# 
+# The equation implemented in this heating model is from Tobie et al. (2025) (https://doi.org/10.1007/s11214-025-01136-y),
+# which is defined as:
+# g = - magnitude - gradient ( - tidal potential ).
+# potential = 3 G M_p R_s^2 / 2 a_s^3 r^2 (Tstar + T0)
+# Tstar = 1/6 *(1-3*cos(theta)^2) and T0=1/2sin(theta)^2*cos(2*lambda + 2*b*t),
+# where G: gravitational constant, M_p: mass of planet, R_s: radius of satellite, a_s: semimajor axis of satellite's orbit, b = angular rate of nonsynchronous rotation.
+# r, theta and lambda are radial distance, polar angle and azimuthal angle, respectively.
+#
+# Model shows the Europa's icy shell without conduction in simpler model.
+# Visualization 'gravity' shows gravity distribution.
+
+set Dimension                              = 3
+set Use years instead of seconds           = true
+set End time                               = 1e4
+
+set Output directory                       = gravity_tidal_potential
+
+set Maximum first time step                    = 1e3
+set CFL number                                 = 0.8
+set Maximum time step                          = 1e3
+
+
+set Pressure normalization                 = surface
+set Surface pressure                       = 0
+
+
+subsection Geometry model
+  set Model name = spherical shell
+  subsection Spherical shell
+    set Outer radius = 1560800
+    set Inner radius = 1460800
+    set Opening angle = 360
+  end
+end
+
+
+subsection Initial temperature model
+  set Model name = function
+  subsection Function
+    set Coordinate system = spherical
+    set Variable names = r, phi,theta
+    set Function expression = 100
+  end
+end
+
+
+subsection Boundary velocity model
+  set Zero velocity boundary indicators = top, bottom
+end
+
+
+subsection Gravity model
+  set Model name = radial with tidal potential
+  subsection Radial constant
+    set Magnitude = 1.3
+  end
+  subsection Radial with tidal potential
+  end
+end
+
+
+subsection Material model
+  set Model name = simpler
+  subsection Simpler model
+    set Reference density = 917
+    set Reference specific heat = 2110
+    set Reference temperature = 100
+    set Thermal conductivity = 0 #1.93
+    set Thermal expansion coefficient = 0 #1.6e-4
+    set Viscosity = 1e20
+  end
+end
+
+
+subsection Formulation
+  set Formulation = Boussinesq approximation
+end
+
+
+subsection Mesh refinement
+  set Initial global refinement                = 0
+  set Initial adaptive refinement              = 0
+  set Time steps between mesh refinement       = 0
+end
+
+
+subsection Postprocess
+  set List of postprocessors = velocity statistics, temperature statistics, visualization, basic statistics, \
+			       pressure statistics, material statistics
+
+  subsection Visualization
+    set Time between graphical output = 1e3
+    set Output format                 = vtu
+    set List of output variables      = material properties, strain rate, shear stress, stress, nonadiabatic pressure, gravity
+  end
+end