diff --git a/doc/modules/changes/20250716_hyunseong96 b/doc/modules/changes/20250716_hyunseong96 new file mode 100644 index 00000000000..360e963c29f --- /dev/null +++ b/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. +
+(Hyunseong Kim, Antoniette Greta Grima, Wolfgang Bangerth 2025/07/16) diff --git a/include/aspect/gravity_model/radial_with_tidal_potential.h b/include/aspect/gravity_model/radial_with_tidal_potential.h new file mode 100644 index 00000000000..6107417f9e0 --- /dev/null +++ b/include/aspect/gravity_model/radial_with_tidal_potential.h @@ -0,0 +1,90 @@ +/* + 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 + . +*/ + + +#ifndef _aspect_gravity_model_radial_with_tidal_potential_h +#define _aspect_gravity_model_radial_with_tidal_potential_h + +#include +#include +#include + +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. + * b = 2 * pi / P where P is period of NSR. + * r, theta and lambda are radial distance, polar angle and azimuthal angle, respectively. + * + * @ingroup GravityModels + */ + template + class RadialWithTidalPotential : public Interface, public SimulatorAccess + { + public: + /** + * Return the gravity vector as a function of position. + */ + Tensor<1,dim> gravity_vector (const Point &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; + + /** + * Period of the non-synchronous rotation in year or second + */ + double P; + }; + } +} + +#endif diff --git a/source/gravity_model/radial_with_tidal_potential.cc b/source/gravity_model/radial_with_tidal_potential.cc new file mode 100644 index 00000000000..0cc840acf3a --- /dev/null +++ b/source/gravity_model/radial_with_tidal_potential.cc @@ -0,0 +1,150 @@ +/* + 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 + . +*/ + + +#include +#include +#include + +#include + +#include + +namespace aspect +{ + namespace GravityModel + { + template + Tensor<1,dim> + RadialWithTidalPotential::gravity_vector (const Point &/*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) + * b = 2 * pi / P + */ + const double t = (this->simulator_is_past_initialization()) ? this->get_time() : 0.0; + + const double angular_frequency = 2. * dealii::numbers::PI / P; + const double C1 = std::cos( 2. * angular_frequency * t); + const double C2 = std::sin( 2. * angular_frequency * 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 radialconstant; + return radialconstant.gravity_vector(p) + tidal_gravity; + } + + + template + void + RadialWithTidalPotential::declare_parameters (ParameterHandler &prm) + { + RadialLinear::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 ("Period of nonsynchronous rotation", "10000", + Patterns::Double (), + "Period of nonsynchronous rotation (NSR). " + "This works for the modeled body having decoupled rotation between interior layers. " + "The default value is the period of NSR on Europa's icy shell. " + "Units is $year$ when 'Use years instead of seconds' is true, " + "and $second$ when 'Use years instead of seconds' is false. "); + } + prm.leave_subsection (); + } + prm.leave_subsection (); + } + + + template + void + RadialWithTidalPotential::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; + P = prm.get_double ("Period of nonsynchronous rotation") * 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'.") + } +} diff --git a/tests/gravity_tidal_potential.prm b/tests/gravity_tidal_potential.prm new file mode 100644 index 00000000000..0e3f47b936c --- /dev/null +++ b/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 diff --git a/tests/gravity_tidal_potential/screen-output b/tests/gravity_tidal_potential/screen-output new file mode 100644 index 00000000000..cc8652dedb8 --- /dev/null +++ b/tests/gravity_tidal_potential/screen-output @@ -0,0 +1,129 @@ + +Number of active cells: 96 (on 1 levels) +Number of degrees of freedom: 4,030 (2,910+150+970) + +*** Timestep 0: t=0 years, dt=0 years + Solving temperature system... 0 iterations. + Solving Stokes system (GMG)... 38+0 iterations. + + Postprocessing: + RMS, max velocity: 2.04e-05 m/year, 7.18e-05 m/year + Temperature min/avg/max: 100 K, 100 K, 100 K + Writing graphical output: output-gravity_tidal_potential/solution/solution-00000 + Pressure min/avg/max: -5.226e+05 Pa, 0.00209 Pa, 1.045e+06 Pa + Average density / Average viscosity / Total mass: 917 kg/m^3, 1e+20 Pa s, 2.631e+21 kg + +*** Timestep 1: t=1000 years, dt=1000 years + Solving temperature system... 0 iterations. + Solving Stokes system (GMG)... 50+0 iterations. + + Postprocessing: + RMS, max velocity: 2.53e-05 m/year, 8.53e-05 m/year + Temperature min/avg/max: 100 K, 100 K, 100 K + Writing graphical output: output-gravity_tidal_potential/solution/solution-00001 + Pressure min/avg/max: -5.226e+05 Pa, 0.00225 Pa, 1.018e+06 Pa + Average density / Average viscosity / Total mass: 917 kg/m^3, 1e+20 Pa s, 2.631e+21 kg + +*** Timestep 2: t=2000 years, dt=1000 years + Solving temperature system... 0 iterations. + Solving Stokes system (GMG)... 50+0 iterations. + + Postprocessing: + RMS, max velocity: 2.24e-05 m/year, 8.31e-05 m/year + Temperature min/avg/max: 100 K, 100 K, 100 K + Writing graphical output: output-gravity_tidal_potential/solution/solution-00002 + Pressure min/avg/max: -5.226e+05 Pa, 0.00225 Pa, 8.955e+05 Pa + Average density / Average viscosity / Total mass: 917 kg/m^3, 1e+20 Pa s, 2.631e+21 kg + +*** Timestep 3: t=3000 years, dt=1000 years + Solving temperature system... 0 iterations. + Solving Stokes system (GMG)... 47+0 iterations. + + Postprocessing: + RMS, max velocity: 2.24e-05 m/year, 8.31e-05 m/year + Temperature min/avg/max: 100 K, 100 K, 100 K + Writing graphical output: output-gravity_tidal_potential/solution/solution-00003 + Pressure min/avg/max: -5.226e+05 Pa, 0.00225 Pa, 8.955e+05 Pa + Average density / Average viscosity / Total mass: 917 kg/m^3, 1e+20 Pa s, 2.631e+21 kg + +*** Timestep 4: t=4000 years, dt=1000 years + Solving temperature system... 0 iterations. + Solving Stokes system (GMG)... 50+0 iterations. + + Postprocessing: + RMS, max velocity: 2.53e-05 m/year, 8.53e-05 m/year + Temperature min/avg/max: 100 K, 100 K, 100 K + Writing graphical output: output-gravity_tidal_potential/solution/solution-00004 + Pressure min/avg/max: -5.226e+05 Pa, 0.00225 Pa, 1.018e+06 Pa + Average density / Average viscosity / Total mass: 917 kg/m^3, 1e+20 Pa s, 2.631e+21 kg + +*** Timestep 5: t=5000 years, dt=1000 years + Solving temperature system... 0 iterations. + Solving Stokes system (GMG)... 50+0 iterations. + + Postprocessing: + RMS, max velocity: 2.04e-05 m/year, 7.18e-05 m/year + Temperature min/avg/max: 100 K, 100 K, 100 K + Writing graphical output: output-gravity_tidal_potential/solution/solution-00005 + Pressure min/avg/max: -5.226e+05 Pa, 0.00225 Pa, 1.045e+06 Pa + Average density / Average viscosity / Total mass: 917 kg/m^3, 1e+20 Pa s, 2.631e+21 kg + +*** Timestep 6: t=6000 years, dt=1000 years + Solving temperature system... 0 iterations. + Solving Stokes system (GMG)... 29+0 iterations. + + Postprocessing: + RMS, max velocity: 2.53e-05 m/year, 8.53e-05 m/year + Temperature min/avg/max: 100 K, 100 K, 100 K + Writing graphical output: output-gravity_tidal_potential/solution/solution-00006 + Pressure min/avg/max: -5.226e+05 Pa, 0.00225 Pa, 1.018e+06 Pa + Average density / Average viscosity / Total mass: 917 kg/m^3, 1e+20 Pa s, 2.631e+21 kg + +*** Timestep 7: t=7000 years, dt=1000 years + Solving temperature system... 0 iterations. + Solving Stokes system (GMG)... 50+0 iterations. + + Postprocessing: + RMS, max velocity: 2.24e-05 m/year, 8.31e-05 m/year + Temperature min/avg/max: 100 K, 100 K, 100 K + Writing graphical output: output-gravity_tidal_potential/solution/solution-00007 + Pressure min/avg/max: -5.226e+05 Pa, 0.00225 Pa, 8.955e+05 Pa + Average density / Average viscosity / Total mass: 917 kg/m^3, 1e+20 Pa s, 2.631e+21 kg + +*** Timestep 8: t=8000 years, dt=1000 years + Solving temperature system... 0 iterations. + Solving Stokes system (GMG)... 50+0 iterations. + + Postprocessing: + RMS, max velocity: 2.24e-05 m/year, 8.31e-05 m/year + Temperature min/avg/max: 100 K, 100 K, 100 K + Writing graphical output: output-gravity_tidal_potential/solution/solution-00008 + Pressure min/avg/max: -5.226e+05 Pa, 0.00225 Pa, 8.955e+05 Pa + Average density / Average viscosity / Total mass: 917 kg/m^3, 1e+20 Pa s, 2.631e+21 kg + +*** Timestep 9: t=9000 years, dt=1000 years + Solving temperature system... 0 iterations. + Solving Stokes system (GMG)... 50+0 iterations. + + Postprocessing: + RMS, max velocity: 2.53e-05 m/year, 8.53e-05 m/year + Temperature min/avg/max: 100 K, 100 K, 100 K + Writing graphical output: output-gravity_tidal_potential/solution/solution-00009 + Pressure min/avg/max: -5.226e+05 Pa, 0.00225 Pa, 1.018e+06 Pa + Average density / Average viscosity / Total mass: 917 kg/m^3, 1e+20 Pa s, 2.631e+21 kg + +*** Timestep 10: t=10000 years, dt=1000 years + Solving temperature system... 0 iterations. + Solving Stokes system (GMG)... 50+0 iterations. + + Postprocessing: + RMS, max velocity: 2.04e-05 m/year, 7.18e-05 m/year + Temperature min/avg/max: 100 K, 100 K, 100 K + Writing graphical output: output-gravity_tidal_potential/solution/solution-00010 + Pressure min/avg/max: -5.226e+05 Pa, 0.00225 Pa, 1.045e+06 Pa + Average density / Average viscosity / Total mass: 917 kg/m^3, 1e+20 Pa s, 2.631e+21 kg + +Termination requested by criterion: end time + + + diff --git a/tests/gravity_tidal_potential/statistics b/tests/gravity_tidal_potential/statistics new file mode 100644 index 00000000000..1340acc34a5 --- /dev/null +++ b/tests/gravity_tidal_potential/statistics @@ -0,0 +1,33 @@ +# 1: Time step number +# 2: Time (years) +# 3: Time step size (years) +# 4: Number of mesh cells +# 5: Number of Stokes degrees of freedom +# 6: Number of temperature degrees of freedom +# 7: Iterations for temperature solver +# 8: Iterations for Stokes solver +# 9: Velocity iterations in Stokes preconditioner +# 10: Schur complement iterations in Stokes preconditioner +# 11: RMS velocity (m/year) +# 12: Max. velocity (m/year) +# 13: Minimal temperature (K) +# 14: Average temperature (K) +# 15: Maximal temperature (K) +# 16: Visualization file name +# 17: Minimal pressure (Pa) +# 18: Average pressure (Pa) +# 19: Maximal pressure (Pa) +# 20: Average density (kg/m^3) +# 21: Average viscosity (Pa s) +# 22: Total mass (kg) + 0 0.000000000000e+00 0.000000000000e+00 96 3060 970 0 38 39 39 2.04002849e-05 7.18187163e-05 1.00000000e+02 9.99999535e+01 1.00000000e+02 output-gravity_tidal_potential/solution/solution-00000 -5.22625857e+05 2.09015736e-03 1.04524968e+06 9.17000000e+02 1.00000000e+20 2.63118556e+21 + 1 1.000000000000e+03 1.000000000000e+03 96 3060 970 0 50 51 51 2.53157545e-05 8.53402100e-05 1.00000000e+02 9.99999535e+01 1.00000000e+02 output-gravity_tidal_potential/solution/solution-00001 -5.22625788e+05 2.24999227e-03 1.01845606e+06 9.17000000e+02 1.00000000e+20 2.63118556e+21 + 2 2.000000000000e+03 1.000000000000e+03 96 3060 970 0 50 51 51 2.24054782e-05 8.30686057e-05 1.00000000e+02 9.99999535e+01 1.00000000e+02 output-gravity_tidal_potential/solution/solution-00002 -5.22625788e+05 2.24991582e-03 8.95531001e+05 9.17000000e+02 1.00000000e+20 2.63118556e+21 + 3 3.000000000000e+03 1.000000000000e+03 96 3060 970 0 47 48 48 2.24054780e-05 8.30686165e-05 1.00000000e+02 9.99999535e+01 1.00000000e+02 output-gravity_tidal_potential/solution/solution-00003 -5.22625787e+05 2.24969716e-03 8.95530999e+05 9.17000000e+02 1.00000000e+20 2.63118556e+21 + 4 4.000000000000e+03 1.000000000000e+03 96 3060 970 0 50 51 51 2.53157551e-05 8.53401992e-05 1.00000000e+02 9.99999535e+01 1.00000000e+02 output-gravity_tidal_potential/solution/solution-00004 -5.22625788e+05 2.25002031e-03 1.01845608e+06 9.17000000e+02 1.00000000e+20 2.63118556e+21 + 5 5.000000000000e+03 1.000000000000e+03 96 3060 970 0 50 51 51 2.04002862e-05 7.18187350e-05 1.00000000e+02 9.99999535e+01 1.00000000e+02 output-gravity_tidal_potential/solution/solution-00005 -5.22625787e+05 2.24982528e-03 1.04524980e+06 9.17000000e+02 1.00000000e+20 2.63118556e+21 + 6 6.000000000000e+03 1.000000000000e+03 96 3060 970 0 29 30 30 2.53157554e-05 8.53402043e-05 1.00000000e+02 9.99999535e+01 1.00000000e+02 output-gravity_tidal_potential/solution/solution-00006 -5.22625787e+05 2.25015373e-03 1.01845607e+06 9.17000000e+02 1.00000000e+20 2.63118556e+21 + 7 7.000000000000e+03 1.000000000000e+03 96 3060 970 0 50 51 51 2.24054782e-05 8.30686059e-05 1.00000000e+02 9.99999535e+01 1.00000000e+02 output-gravity_tidal_potential/solution/solution-00007 -5.22625788e+05 2.25002696e-03 8.95531001e+05 9.17000000e+02 1.00000000e+20 2.63118556e+21 + 8 8.000000000000e+03 1.000000000000e+03 96 3060 970 0 50 51 51 2.24054786e-05 8.30686277e-05 1.00000000e+02 9.99999535e+01 1.00000000e+02 output-gravity_tidal_potential/solution/solution-00008 -5.22625788e+05 2.24991928e-03 8.95530999e+05 9.17000000e+02 1.00000000e+20 2.63118556e+21 + 9 9.000000000000e+03 1.000000000000e+03 96 3060 970 0 50 51 51 2.53157551e-05 8.53401992e-05 1.00000000e+02 9.99999535e+01 1.00000000e+02 output-gravity_tidal_potential/solution/solution-00009 -5.22625788e+05 2.25004559e-03 1.01845608e+06 9.17000000e+02 1.00000000e+20 2.63118556e+21 +10 1.000000000000e+04 1.000000000000e+03 96 3060 970 0 50 51 51 2.04002862e-05 7.18187350e-05 1.00000000e+02 9.99999535e+01 1.00000000e+02 output-gravity_tidal_potential/solution/solution-00010 -5.22625787e+05 2.25020422e-03 1.04524980e+06 9.17000000e+02 1.00000000e+20 2.63118556e+21