From f6a4b3dfbcbc00ff2cb849a4b1c16281efb1be85 Mon Sep 17 00:00:00 2001 From: Hyunseong Kim Date: Fri, 13 Jun 2025 17:23:14 -0600 Subject: [PATCH 01/21] initiate gravity model for tidal effect --- include/aspect/gravity_model/radial.h | 46 ++++++ source/gravity_model/radial.cc | 170 ++++++++++++++++++++ tests/gravity_tidal_potential.prm | 100 ++++++++++++ tests/gravity_tidal_potential/screen-output | 139 ++++++++++++++++ tests/gravity_tidal_potential/statistics | 33 ++++ 5 files changed, 488 insertions(+) create mode 100644 tests/gravity_tidal_potential.prm create mode 100644 tests/gravity_tidal_potential/screen-output create mode 100644 tests/gravity_tidal_potential/statistics diff --git a/include/aspect/gravity_model/radial.h b/include/aspect/gravity_model/radial.h index 07059b5ecbc..3519f21b268 100644 --- a/include/aspect/gravity_model/radial.h +++ b/include/aspect/gravity_model/radial.h @@ -135,6 +135,52 @@ namespace aspect double magnitude_at_bottom; }; + + + + template + class RadialLinearWithTidalPotential : public RadialLinear + { + 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 hosting planet + */ + double M_p; + + /** + * Semimajor axis of satellite's orbit + */ + double a_s; + + /** + * Rate of non-synchronous rotation in m/year + */ + double b_NSR; + + /** + * Radius of modelling satellite + */ + double R1; + }; } } diff --git a/source/gravity_model/radial.cc b/source/gravity_model/radial.cc index 86b59a38141..5188b5210b0 100644 --- a/source/gravity_model/radial.cc +++ b/source/gravity_model/radial.cc @@ -25,6 +25,11 @@ #include +#include +#include +#include +#include + namespace aspect { namespace GravityModel @@ -178,6 +183,159 @@ namespace aspect "do not have either a spherical or ellipsoidal natural coordinate system.")); } + + + +// ----------------------------- RadialLinearWithTidalPotential ---------------------- + + + template + Tensor<1,dim> + RadialLinearWithTidalPotential::gravity_vector (const Point &p) const + { + /** + * This gravity model is a plugin for a case where the + * tidal potential by flattening and non-synchronnous rotation chages 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) / (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, 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. + * b [1/s] = b_NSR [m/yr] / (circumference of satellite [m]) / year_to_seconds [s/yr] * 2 * pi + * + * 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 Tensor<1,dim> e_x = Point::unit_vector(0); + const Tensor<1,dim> e_y = Point::unit_vector(1); + const Tensor<1,dim> e_z = Point::unit_vector(2); + + const double t = (this->simulator_is_past_initialization()) ? this->get_time() : 0.0; + const double x = p[0]; + const double y = p[1]; + const double z = p[2]; + + const double C1 = std::cos( 2. * b_NSR / R1 / year_in_seconds * t); + const double C2 = std::sin( 2. * b_NSR / R1 / year_in_seconds * t); + + const double dTstar_over_dx = 1. / 6. * ( 2. * x ); + const double dT0_over_dx = 1. / 2. * ( 2. * C1 * x - 2. * C2 * y ); + + const double dTstar_over_dy = 1. / 6. * ( 2. * y ); + const double dT0_over_dy = 1. / 2. * ( -2. * C1 * y - 2. * C2 * x ); + + const double dTstar_over_dz = 1. / 6. * ( -4. * z ); + const double dT0_over_dz = 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_over_dx + dT0_over_dx) * e_x + + (dTstar_over_dy + dT0_over_dy) * e_y + + (dTstar_over_dz + dT0_over_dz) * e_z); + + return RadialLinear::gravity_vector(p) + tidal_gravity; + } + + + template + void + RadialLinearWithTidalPotential::declare_parameters (ParameterHandler &prm) + { + RadialLinear::declare_parameters (prm); + + prm.enter_subsection("Gravity model"); + { + prm.enter_subsection("Radial linear with tidal potnetial"); + { + prm.declare_entry ("Mass of planet", "1.898e27", + Patterns::Double (), + "Mass of satellte's host planet. " + "Default value is for Europa, therefore, mass of Jupiter. " + "Units is $kg$."); + prm.declare_entry ("Semimajor axis of satellite", "670900000", + Patterns::Double (), + "Length of semimajor axis of satellite. " + "Default value is for Europa's semimajor axis" + "Units is $m$."); + prm.declare_entry ("Rate of nonsynchronous rotation", "1000", + Patterns::Double (), + "Rate of nonsynchronous rotation (NSR). " + "This will be converted to angular rate. " + "Units is $m/year$"); + } + prm.leave_subsection (); + } + prm.leave_subsection (); + } + + + template + void + RadialLinearWithTidalPotential::parse_parameters (ParameterHandler &prm) + { + AssertThrow (dim==3, ExcMessage ("The 'radial linear with tidal potential' gravity model " + "can only be used in 3D.")); + this->RadialLinear::parse_parameters (prm); + + prm.enter_subsection("Gravity model"); + { + prm.enter_subsection("Radial linear with tidal potnetial"); + { + M_p = prm.get_double ("Mass of planet"); + a_s = prm.get_double ("Semimajor axis of satellite"); + b_NSR = prm.get_double ("Rate of nonsynchronous rotation"); + } + prm.leave_subsection (); + } + prm.leave_subsection (); + // This effect of tidal potential only works if the geometry is derived from + // a spherical model (i.e. a sphere, spherical shell or chunk) + if (Plugins::plugin_type_matches>(this->get_geometry_model())) + { + R1 = Plugins::get_plugin_as_type> + (this->get_geometry_model()).radius(); + } + else if (Plugins::plugin_type_matches>(this->get_geometry_model())) + { + R1 = Plugins::get_plugin_as_type> + (this->get_geometry_model()).outer_radius(); + } + else if (Plugins::plugin_type_matches>(this->get_geometry_model())) + { + R1 = Plugins::get_plugin_as_type> + (this->get_geometry_model()).outer_radius(); + } + else if (Plugins::plugin_type_matches>(this->get_geometry_model())) + { + const auto &gm = Plugins::get_plugin_as_type> + (this->get_geometry_model()); + // TODO + // If the eccentricity of the EllipsoidalChunk is non-zero, the radius can vary along a boundary, + // but the maximal depth is the same everywhere and we could calculate a representative pressure + // profile. However, it requires some extra logic with ellipsoidal + // coordinates, so for now we only allow eccentricity zero. + // Using the EllipsoidalChunk with eccentricity zero can still be useful, + // because the domain can be non-coordinate parallel. + + AssertThrow(gm.get_eccentricity() == 0.0, + ExcNotImplemented("This plugin cannot be used with a non-zero eccentricity. ")); + + R1 = gm.get_semi_major_axis_a(); + } + else + { + Assert (false, ExcMessage ("This initial condition can only be used if the geometry " + "is a sphere, a spherical shell, a chunk or an " + "ellipsoidal chunk.")); + R1 = numbers::signaling_nan(); + } + } } } @@ -207,5 +365,17 @@ namespace aspect "is positive) and the magnitude changes linearly with depth. The " "magnitude of gravity at the surface and bottom is read from the " "input file in a section ``Gravity model/Radial linear''.") + + ASPECT_REGISTER_GRAVITY_MODEL(RadialLinearWithTidalPotential, + "radial linear with tidal potential", + "A gravity model that is the sum of the `radial linear' model " + "(which is radial, pointing inward if the gravity " + "is positive, and a magnitude that changes linearly with depth), " + "and a term that results from a tidal gravity potential and that " + "leads to a gravity field that varies with latitude and longitude." + "The magnitude of gravity for the linear radial part is read from the " + "input file in a section `Gravity model/Radial linear'; the " + "parameters that describe the tidal potential contribution are read " + "from a section ``Gravity model/Radial linear with tidal potential''.") } } diff --git a/tests/gravity_tidal_potential.prm b/tests/gravity_tidal_potential.prm new file mode 100644 index 00000000000..08f9bca6820 --- /dev/null +++ b/tests/gravity_tidal_potential.prm @@ -0,0 +1,100 @@ +# This parameter file tests the gravity model plugin for a case where the +# tidal potential by flattening and non-synchronnous rotation chages 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 of 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 in output instead of seconds = true +set End time = 5e5 + +set Output directory = gravity_tidal_potential + +set Maximum first time step = 5e4 +set CFL number = 0.8 +set Maximum time step = 5e4 + + +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 linear with tidal potential + subsection Radial linear + set Magnitude at surface = 1.3 + set Magnitude at bottom = 1.3 + end + subsection Radial linear with tidal potnetial + 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 = 5e4 + 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..a2c902a1fd1 --- /dev/null +++ b/tests/gravity_tidal_potential/screen-output @@ -0,0 +1,139 @@ +----------------------------------------------------------------------------- +----------------------------------------------------------------------------- +----------------------------------------------------------------------------- + +----------------------------------------------------------------------------- +----------------------------------------------------------------------------- +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)... 110+0 iterations. + + Postprocessing: + RMS, max velocity: 2.04e-05 m/year, 7.16e-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.228e+05 Pa, 5.829e+07 Pa, 1.201e+08 Pa + Average density / Average viscosity / Total mass: 917 kg/m^3, 1e+20 Pa s, 2.631e+21 kg + +*** Timestep 1: t=50000 years, dt=50000 years + Solving temperature system... 0 iterations. + Solving Stokes system (GMG)... 50+0 iterations. + + Postprocessing: + RMS, max velocity: 2.53e-05 m/year, 8.48e-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.228e+05 Pa, 5.829e+07 Pa, 1.201e+08 Pa + Average density / Average viscosity / Total mass: 917 kg/m^3, 1e+20 Pa s, 2.631e+21 kg + +*** Timestep 2: t=100000 years, dt=50000 years + Solving temperature system... 0 iterations. + Solving Stokes system (GMG)... 50+0 iterations. + + Postprocessing: + RMS, max velocity: 2.26e-05 m/year, 8.32e-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.233e+05 Pa, 5.829e+07 Pa, 1.2e+08 Pa + Average density / Average viscosity / Total mass: 917 kg/m^3, 1e+20 Pa s, 2.631e+21 kg + +*** Timestep 3: t=150000 years, dt=50000 years + Solving temperature system... 0 iterations. + Solving Stokes system (GMG)... 50+0 iterations. + + Postprocessing: + RMS, max velocity: 2.21e-05 m/year, 8.26e-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.228e+05 Pa, 5.829e+07 Pa, 1.2e+08 Pa + Average density / Average viscosity / Total mass: 917 kg/m^3, 1e+20 Pa s, 2.631e+21 kg + +*** Timestep 4: t=200000 years, dt=50000 years + Solving temperature system... 0 iterations. + Solving Stokes system (GMG)... 50+0 iterations. + + Postprocessing: + RMS, max velocity: 2.55e-05 m/year, 8.72e-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.259e+05 Pa, 5.829e+07 Pa, 1.201e+08 Pa + Average density / Average viscosity / Total mass: 917 kg/m^3, 1e+20 Pa s, 2.631e+21 kg + +*** Timestep 5: t=250000 years, dt=50000 years + Solving temperature system... 0 iterations. + Solving Stokes system (GMG)... 50+0 iterations. + + Postprocessing: + RMS, max velocity: 2.05e-05 m/year, 7.37e-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.228e+05 Pa, 5.829e+07 Pa, 1.201e+08 Pa + Average density / Average viscosity / Total mass: 917 kg/m^3, 1e+20 Pa s, 2.631e+21 kg + +*** Timestep 6: t=300000 years, dt=50000 years + Solving temperature system... 0 iterations. + Solving Stokes system (GMG)... 49+0 iterations. + + Postprocessing: + RMS, max velocity: 2.49e-05 m/year, 8.18e-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.228e+05 Pa, 5.829e+07 Pa, 1.201e+08 Pa + Average density / Average viscosity / Total mass: 917 kg/m^3, 1e+20 Pa s, 2.631e+21 kg + +*** Timestep 7: t=350000 years, dt=50000 years + Solving temperature system... 0 iterations. + Solving Stokes system (GMG)... 50+0 iterations. + + Postprocessing: + RMS, max velocity: 2.31e-05 m/year, 8.34e-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.229e+05 Pa, 5.829e+07 Pa, 1.2e+08 Pa + Average density / Average viscosity / Total mass: 917 kg/m^3, 1e+20 Pa s, 2.631e+21 kg + +*** Timestep 8: t=400000 years, dt=50000 years + Solving temperature system... 0 iterations. + Solving Stokes system (GMG)... 50+0 iterations. + + Postprocessing: + RMS, max velocity: 2.17e-05 m/year, 8.16e-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.228e+05 Pa, 5.829e+07 Pa, 1.2e+08 Pa + Average density / Average viscosity / Total mass: 917 kg/m^3, 1e+20 Pa s, 2.631e+21 kg + +*** Timestep 9: t=450000 years, dt=50000 years + Solving temperature system... 0 iterations. + Solving Stokes system (GMG)... 50+0 iterations. + + Postprocessing: + RMS, max velocity: 2.57e-05 m/year, 8.87e-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.29e+05 Pa, 5.829e+07 Pa, 1.201e+08 Pa + Average density / Average viscosity / Total mass: 917 kg/m^3, 1e+20 Pa s, 2.631e+21 kg + +*** Timestep 10: t=500000 years, dt=50000 years + Solving temperature system... 0 iterations. + Solving Stokes system (GMG)... 50+0 iterations. + + Postprocessing: + RMS, max velocity: 2.06e-05 m/year, 7.55e-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.228e+05 Pa, 5.829e+07 Pa, 1.201e+08 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..c5d50b9834b --- /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 110 112 112 2.04033429e-05 7.15757443e-05 1.00000000e+02 9.99999535e+01 1.00000000e+02 output-gravity_tidal_potential/solution/solution-00000 -5.22786493e+05 5.82904414e+07 1.20123467e+08 9.17000000e+02 1.00000000e+20 2.63118556e+21 + 1 5.000000000000e+04 5.000000000000e+04 96 3060 970 0 50 51 51 2.52626040e-05 8.48256428e-05 1.00000000e+02 9.99999535e+01 1.00000000e+02 output-gravity_tidal_potential/solution/solution-00001 -5.22786489e+05 5.82904414e+07 1.20095826e+08 9.17000000e+02 1.00000000e+20 2.63118556e+21 + 2 1.000000000000e+05 5.000000000000e+04 96 3060 970 0 50 51 51 2.26076512e-05 8.32470559e-05 1.00000000e+02 9.99999535e+01 1.00000000e+02 output-gravity_tidal_potential/solution/solution-00002 -5.23328149e+05 5.82904414e+07 1.19978315e+08 9.17000000e+02 1.00000000e+20 2.63118556e+21 + 3 1.500000000000e+05 5.000000000000e+04 96 3060 970 0 50 51 51 2.21153389e-05 8.26357877e-05 1.00000000e+02 9.99999535e+01 1.00000000e+02 output-gravity_tidal_potential/solution/solution-00003 -5.22786490e+05 5.82904414e+07 1.20014403e+08 9.17000000e+02 1.00000000e+20 2.63118556e+21 + 4 2.000000000000e+05 5.000000000000e+04 96 3060 970 0 50 51 51 2.55109712e-05 8.71612001e-05 1.00000000e+02 9.99999535e+01 1.00000000e+02 output-gravity_tidal_potential/solution/solution-00004 -5.25872516e+05 5.82904414e+07 1.20114229e+08 9.17000000e+02 1.00000000e+20 2.63118556e+21 + 5 2.500000000000e+05 5.000000000000e+04 96 3060 970 0 50 51 51 2.04568088e-05 7.37033618e-05 1.00000000e+02 9.99999535e+01 1.00000000e+02 output-gravity_tidal_potential/solution/solution-00005 -5.22786489e+05 5.82904414e+07 1.20120447e+08 9.17000000e+02 1.00000000e+20 2.63118556e+21 + 6 3.000000000000e+05 5.000000000000e+04 96 3060 970 0 49 50 50 2.49429919e-05 8.17556197e-05 1.00000000e+02 9.99999535e+01 1.00000000e+02 output-gravity_tidal_potential/solution/solution-00006 -5.22786489e+05 5.82904414e+07 1.20071635e+08 9.17000000e+02 1.00000000e+20 2.63118556e+21 + 7 3.500000000000e+05 5.000000000000e+04 96 3060 970 0 50 51 51 2.31118595e-05 8.33796208e-05 1.00000000e+02 9.99999535e+01 1.00000000e+02 output-gravity_tidal_potential/solution/solution-00007 -5.22887837e+05 5.82904414e+07 1.19977948e+08 9.17000000e+02 1.00000000e+20 2.63118556e+21 + 8 4.000000000000e+05 5.000000000000e+04 96 3060 970 0 50 51 51 2.16536241e-05 8.15533243e-05 1.00000000e+02 9.99999535e+01 1.00000000e+02 output-gravity_tidal_potential/solution/solution-00008 -5.22786490e+05 5.82904414e+07 1.20046645e+08 9.17000000e+02 1.00000000e+20 2.63118556e+21 + 9 4.500000000000e+05 5.000000000000e+04 96 3060 970 0 50 51 51 2.56815240e-05 8.87416779e-05 1.00000000e+02 9.99999535e+01 1.00000000e+02 output-gravity_tidal_potential/solution/solution-00009 -5.29005056e+05 5.82904414e+07 1.20126681e+08 9.17000000e+02 1.00000000e+20 2.63118556e+21 +10 5.000000000000e+05 5.000000000000e+04 96 3060 970 0 50 51 51 2.06145125e-05 7.55101235e-05 1.00000000e+02 9.99999535e+01 1.00000000e+02 output-gravity_tidal_potential/solution/solution-00010 -5.22786490e+05 5.82904414e+07 1.20111416e+08 9.17000000e+02 1.00000000e+20 2.63118556e+21 From 1ecc8f3357b0cb73fe7b4b3b3b8f408fdd1c3d71 Mon Sep 17 00:00:00 2001 From: Hyunseong Kim Date: Wed, 16 Jul 2025 12:39:34 +0100 Subject: [PATCH 02/21] modifying files --- doc/modules/changes/20250716_hyunseong96 | 4 ++ source/gravity_model/radial.cc | 4 +- tests/gravity_tidal_potential.prm | 10 ++-- tests/gravity_tidal_potential/screen-output | 54 ++++++++++----------- tests/gravity_tidal_potential/statistics | 20 ++++---- 5 files changed, 48 insertions(+), 44 deletions(-) create mode 100644 doc/modules/changes/20250716_hyunseong96 diff --git a/doc/modules/changes/20250716_hyunseong96 b/doc/modules/changes/20250716_hyunseong96 new file mode 100644 index 00000000000..c1a0f9bebce --- /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 potnetial' for tidal force. +This plugin is useful for modeling long-term interior and surface evolution in moons orbiting a large planet. +
+(Hyunseong Kim, 2025/07/16) diff --git a/source/gravity_model/radial.cc b/source/gravity_model/radial.cc index 5188b5210b0..c0a0b85c0b4 100644 --- a/source/gravity_model/radial.cc +++ b/source/gravity_model/radial.cc @@ -195,11 +195,11 @@ namespace aspect { /** * This gravity model is a plugin for a case where the - * tidal potential by flattening and non-synchronnous rotation chages gravity with position and time. + * 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)). + * g = -magnitude - gradient (-tidal potential)). * 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 planet, a_s = semimajor axis of satellite's orbit, b = angular rate of nonsynchronous rotation. diff --git a/tests/gravity_tidal_potential.prm b/tests/gravity_tidal_potential.prm index 08f9bca6820..d5c4fdc2ecc 100644 --- a/tests/gravity_tidal_potential.prm +++ b/tests/gravity_tidal_potential.prm @@ -3,7 +3,7 @@ # # 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 of tidal potential. +# 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. @@ -14,13 +14,13 @@ set Dimension = 3 set Use years in output instead of seconds = true -set End time = 5e5 +set End time = 1e4 set Output directory = gravity_tidal_potential -set Maximum first time step = 5e4 +set Maximum first time step = 1e3 set CFL number = 0.8 -set Maximum time step = 5e4 +set Maximum time step = 1e3 set Pressure normalization = surface @@ -93,7 +93,7 @@ subsection Postprocess pressure statistics, material statistics subsection Visualization - set Time between graphical output = 5e4 + 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 diff --git a/tests/gravity_tidal_potential/screen-output b/tests/gravity_tidal_potential/screen-output index a2c902a1fd1..2779dd879b0 100644 --- a/tests/gravity_tidal_potential/screen-output +++ b/tests/gravity_tidal_potential/screen-output @@ -18,111 +18,111 @@ Number of degrees of freedom: 4,030 (2,910+150+970) Pressure min/avg/max: -5.228e+05 Pa, 5.829e+07 Pa, 1.201e+08 Pa Average density / Average viscosity / Total mass: 917 kg/m^3, 1e+20 Pa s, 2.631e+21 kg -*** Timestep 1: t=50000 years, dt=50000 years +*** 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.48e-05 m/year + RMS, max velocity: 2.54e-05 m/year, 8.6e-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.228e+05 Pa, 5.829e+07 Pa, 1.201e+08 Pa + Pressure min/avg/max: -5.236e+05 Pa, 5.829e+07 Pa, 1.201e+08 Pa Average density / Average viscosity / Total mass: 917 kg/m^3, 1e+20 Pa s, 2.631e+21 kg -*** Timestep 2: t=100000 years, dt=50000 years +*** 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.26e-05 m/year, 8.32e-05 m/year + RMS, max velocity: 2.22e-05 m/year, 8.27e-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.233e+05 Pa, 5.829e+07 Pa, 1.2e+08 Pa + Pressure min/avg/max: -5.228e+05 Pa, 5.829e+07 Pa, 1.2e+08 Pa Average density / Average viscosity / Total mass: 917 kg/m^3, 1e+20 Pa s, 2.631e+21 kg -*** Timestep 3: t=150000 years, dt=50000 years +*** Timestep 3: t=3000 years, dt=1000 years Solving temperature system... 0 iterations. Solving Stokes system (GMG)... 50+0 iterations. Postprocessing: - RMS, max velocity: 2.21e-05 m/year, 8.26e-05 m/year + RMS, max velocity: 2.28e-05 m/year, 8.34e-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.228e+05 Pa, 5.829e+07 Pa, 1.2e+08 Pa + Pressure min/avg/max: -5.239e+05 Pa, 5.829e+07 Pa, 1.2e+08 Pa Average density / Average viscosity / Total mass: 917 kg/m^3, 1e+20 Pa s, 2.631e+21 kg -*** Timestep 4: t=200000 years, dt=50000 years +*** 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.55e-05 m/year, 8.72e-05 m/year + RMS, max velocity: 2.5e-05 m/year, 8.22e-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.259e+05 Pa, 5.829e+07 Pa, 1.201e+08 Pa + Pressure min/avg/max: -5.228e+05 Pa, 5.829e+07 Pa, 1.201e+08 Pa Average density / Average viscosity / Total mass: 917 kg/m^3, 1e+20 Pa s, 2.631e+21 kg -*** Timestep 5: t=250000 years, dt=50000 years +*** 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.05e-05 m/year, 7.37e-05 m/year + RMS, max velocity: 2.05e-05 m/year, 7.43e-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.228e+05 Pa, 5.829e+07 Pa, 1.201e+08 Pa Average density / Average viscosity / Total mass: 917 kg/m^3, 1e+20 Pa s, 2.631e+21 kg -*** Timestep 6: t=300000 years, dt=50000 years +*** Timestep 6: t=6000 years, dt=1000 years Solving temperature system... 0 iterations. Solving Stokes system (GMG)... 49+0 iterations. Postprocessing: - RMS, max velocity: 2.49e-05 m/year, 8.18e-05 m/year + RMS, max velocity: 2.57e-05 m/year, 8.85e-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.228e+05 Pa, 5.829e+07 Pa, 1.201e+08 Pa + Pressure min/avg/max: -5.285e+05 Pa, 5.829e+07 Pa, 1.201e+08 Pa Average density / Average viscosity / Total mass: 917 kg/m^3, 1e+20 Pa s, 2.631e+21 kg -*** Timestep 7: t=350000 years, dt=50000 years +*** 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.31e-05 m/year, 8.34e-05 m/year + RMS, max velocity: 2.15e-05 m/year, 8.12e-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.229e+05 Pa, 5.829e+07 Pa, 1.2e+08 Pa + Pressure min/avg/max: -5.228e+05 Pa, 5.829e+07 Pa, 1.201e+08 Pa Average density / Average viscosity / Total mass: 917 kg/m^3, 1e+20 Pa s, 2.631e+21 kg -*** Timestep 8: t=400000 years, dt=50000 years +*** 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.17e-05 m/year, 8.16e-05 m/year + RMS, max velocity: 2.35e-05 m/year, 8.32e-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.228e+05 Pa, 5.829e+07 Pa, 1.2e+08 Pa Average density / Average viscosity / Total mass: 917 kg/m^3, 1e+20 Pa s, 2.631e+21 kg -*** Timestep 9: t=450000 years, dt=50000 years +*** 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.57e-05 m/year, 8.87e-05 m/year + RMS, max velocity: 2.45e-05 m/year, 8.1e-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.29e+05 Pa, 5.829e+07 Pa, 1.201e+08 Pa + Pressure min/avg/max: -5.228e+05 Pa, 5.829e+07 Pa, 1.2e+08 Pa Average density / Average viscosity / Total mass: 917 kg/m^3, 1e+20 Pa s, 2.631e+21 kg -*** Timestep 10: t=500000 years, dt=50000 years +*** 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.06e-05 m/year, 7.55e-05 m/year + RMS, max velocity: 2.08e-05 m/year, 7.7e-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.228e+05 Pa, 5.829e+07 Pa, 1.201e+08 Pa diff --git a/tests/gravity_tidal_potential/statistics b/tests/gravity_tidal_potential/statistics index c5d50b9834b..b7210a2a72a 100644 --- a/tests/gravity_tidal_potential/statistics +++ b/tests/gravity_tidal_potential/statistics @@ -21,13 +21,13 @@ # 21: Average viscosity (Pa s) # 22: Total mass (kg) 0 0.000000000000e+00 0.000000000000e+00 96 3060 970 0 110 112 112 2.04033429e-05 7.15757443e-05 1.00000000e+02 9.99999535e+01 1.00000000e+02 output-gravity_tidal_potential/solution/solution-00000 -5.22786493e+05 5.82904414e+07 1.20123467e+08 9.17000000e+02 1.00000000e+20 2.63118556e+21 - 1 5.000000000000e+04 5.000000000000e+04 96 3060 970 0 50 51 51 2.52626040e-05 8.48256428e-05 1.00000000e+02 9.99999535e+01 1.00000000e+02 output-gravity_tidal_potential/solution/solution-00001 -5.22786489e+05 5.82904414e+07 1.20095826e+08 9.17000000e+02 1.00000000e+20 2.63118556e+21 - 2 1.000000000000e+05 5.000000000000e+04 96 3060 970 0 50 51 51 2.26076512e-05 8.32470559e-05 1.00000000e+02 9.99999535e+01 1.00000000e+02 output-gravity_tidal_potential/solution/solution-00002 -5.23328149e+05 5.82904414e+07 1.19978315e+08 9.17000000e+02 1.00000000e+20 2.63118556e+21 - 3 1.500000000000e+05 5.000000000000e+04 96 3060 970 0 50 51 51 2.21153389e-05 8.26357877e-05 1.00000000e+02 9.99999535e+01 1.00000000e+02 output-gravity_tidal_potential/solution/solution-00003 -5.22786490e+05 5.82904414e+07 1.20014403e+08 9.17000000e+02 1.00000000e+20 2.63118556e+21 - 4 2.000000000000e+05 5.000000000000e+04 96 3060 970 0 50 51 51 2.55109712e-05 8.71612001e-05 1.00000000e+02 9.99999535e+01 1.00000000e+02 output-gravity_tidal_potential/solution/solution-00004 -5.25872516e+05 5.82904414e+07 1.20114229e+08 9.17000000e+02 1.00000000e+20 2.63118556e+21 - 5 2.500000000000e+05 5.000000000000e+04 96 3060 970 0 50 51 51 2.04568088e-05 7.37033618e-05 1.00000000e+02 9.99999535e+01 1.00000000e+02 output-gravity_tidal_potential/solution/solution-00005 -5.22786489e+05 5.82904414e+07 1.20120447e+08 9.17000000e+02 1.00000000e+20 2.63118556e+21 - 6 3.000000000000e+05 5.000000000000e+04 96 3060 970 0 49 50 50 2.49429919e-05 8.17556197e-05 1.00000000e+02 9.99999535e+01 1.00000000e+02 output-gravity_tidal_potential/solution/solution-00006 -5.22786489e+05 5.82904414e+07 1.20071635e+08 9.17000000e+02 1.00000000e+20 2.63118556e+21 - 7 3.500000000000e+05 5.000000000000e+04 96 3060 970 0 50 51 51 2.31118595e-05 8.33796208e-05 1.00000000e+02 9.99999535e+01 1.00000000e+02 output-gravity_tidal_potential/solution/solution-00007 -5.22887837e+05 5.82904414e+07 1.19977948e+08 9.17000000e+02 1.00000000e+20 2.63118556e+21 - 8 4.000000000000e+05 5.000000000000e+04 96 3060 970 0 50 51 51 2.16536241e-05 8.15533243e-05 1.00000000e+02 9.99999535e+01 1.00000000e+02 output-gravity_tidal_potential/solution/solution-00008 -5.22786490e+05 5.82904414e+07 1.20046645e+08 9.17000000e+02 1.00000000e+20 2.63118556e+21 - 9 4.500000000000e+05 5.000000000000e+04 96 3060 970 0 50 51 51 2.56815240e-05 8.87416779e-05 1.00000000e+02 9.99999535e+01 1.00000000e+02 output-gravity_tidal_potential/solution/solution-00009 -5.29005056e+05 5.82904414e+07 1.20126681e+08 9.17000000e+02 1.00000000e+20 2.63118556e+21 -10 5.000000000000e+05 5.000000000000e+04 96 3060 970 0 50 51 51 2.06145125e-05 7.55101235e-05 1.00000000e+02 9.99999535e+01 1.00000000e+02 output-gravity_tidal_potential/solution/solution-00010 -5.22786490e+05 5.82904414e+07 1.20111416e+08 9.17000000e+02 1.00000000e+20 2.63118556e+21 + 1 1.000000000000e+03 1.000000000000e+03 96 3060 970 0 50 51 51 2.53870595e-05 8.60011863e-05 1.00000000e+02 9.99999535e+01 1.00000000e+02 output-gravity_tidal_potential/solution/solution-00001 -5.23573285e+05 5.82904414e+07 1.20105088e+08 9.17000000e+02 1.00000000e+20 2.63118556e+21 + 2 2.000000000000e+03 1.000000000000e+03 96 3060 970 0 50 51 51 2.21503532e-05 8.26962512e-05 1.00000000e+02 9.99999535e+01 1.00000000e+02 output-gravity_tidal_potential/solution/solution-00002 -5.22786489e+05 5.82904414e+07 1.20011855e+08 9.17000000e+02 1.00000000e+20 2.63118556e+21 + 3 3.000000000000e+03 1.000000000000e+03 96 3060 970 0 50 51 51 2.28043260e-05 8.33561483e-05 1.00000000e+02 9.99999535e+01 1.00000000e+02 output-gravity_tidal_potential/solution/solution-00003 -5.23892834e+05 5.82904414e+07 1.19978816e+08 9.17000000e+02 1.00000000e+20 2.63118556e+21 + 4 4.000000000000e+03 1.000000000000e+03 96 3060 970 0 50 51 51 2.49936170e-05 8.22470390e-05 1.00000000e+02 9.99999535e+01 1.00000000e+02 output-gravity_tidal_potential/solution/solution-00004 -5.22786489e+05 5.82904414e+07 1.20075508e+08 9.17000000e+02 1.00000000e+20 2.63118556e+21 + 5 5.000000000000e+03 1.000000000000e+03 96 3060 970 0 50 51 51 2.04959459e-05 7.43142701e-05 1.00000000e+02 9.99999535e+01 1.00000000e+02 output-gravity_tidal_potential/solution/solution-00005 -5.22786489e+05 5.82904414e+07 1.20118224e+08 9.17000000e+02 1.00000000e+20 2.63118556e+21 + 6 6.000000000000e+03 1.000000000000e+03 96 3060 970 0 49 50 50 2.56511938e-05 8.84619423e-05 1.00000000e+02 9.99999535e+01 1.00000000e+02 output-gravity_tidal_potential/solution/solution-00006 -5.28450605e+05 5.82904414e+07 1.20124477e+08 9.17000000e+02 1.00000000e+20 2.63118556e+21 + 7 7.000000000000e+03 1.000000000000e+03 96 3060 970 0 50 51 51 2.15471702e-05 8.12170166e-05 1.00000000e+02 9.99999535e+01 1.00000000e+02 output-gravity_tidal_potential/solution/solution-00007 -5.22786490e+05 5.82904414e+07 1.20053747e+08 9.17000000e+02 1.00000000e+20 2.63118556e+21 + 8 8.000000000000e+03 1.000000000000e+03 96 3060 970 0 50 51 51 2.34666825e-05 8.31817499e-05 1.00000000e+02 9.99999535e+01 1.00000000e+02 output-gravity_tidal_potential/solution/solution-00008 -5.22786490e+05 5.82904414e+07 1.19974416e+08 9.17000000e+02 1.00000000e+20 2.63118556e+21 + 9 9.000000000000e+03 1.000000000000e+03 96 3060 970 0 50 51 51 2.44893181e-05 8.10214829e-05 1.00000000e+02 9.99999535e+01 1.00000000e+02 output-gravity_tidal_potential/solution/solution-00009 -5.22786489e+05 5.82904414e+07 1.20036188e+08 9.17000000e+02 1.00000000e+20 2.63118556e+21 +10 1.000000000000e+04 1.000000000000e+03 96 3060 970 0 50 51 51 2.07657183e-05 7.69853998e-05 1.00000000e+02 9.99999535e+01 1.00000000e+02 output-gravity_tidal_potential/solution/solution-00010 -5.22786489e+05 5.82904414e+07 1.20102575e+08 9.17000000e+02 1.00000000e+20 2.63118556e+21 From c8e0290e0c9a27c1acade603b14ea8948defaac9 Mon Sep 17 00:00:00 2001 From: Hyunseong Kim Date: Wed, 16 Jul 2025 12:42:00 +0100 Subject: [PATCH 03/21] type correction --- tests/gravity_tidal_potential.prm | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/gravity_tidal_potential.prm b/tests/gravity_tidal_potential.prm index d5c4fdc2ecc..491312ef187 100644 --- a/tests/gravity_tidal_potential.prm +++ b/tests/gravity_tidal_potential.prm @@ -1,5 +1,5 @@ # This parameter file tests the gravity model plugin for a case where the -# tidal potential by flattening and non-synchronnous rotation chages gravity with position and time. +# 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: From 8ba55bac423cedc9bcd828b39a5aae597d238d83 Mon Sep 17 00:00:00 2001 From: Hyunseong Kim Date: Tue, 29 Jul 2025 21:32:42 +0100 Subject: [PATCH 04/21] test apply diff --- tests/gravity_tidal_potential/screen-output | 12 +---------- tests/gravity_tidal_potential/statistics | 22 ++++++++++----------- 2 files changed, 12 insertions(+), 22 deletions(-) diff --git a/tests/gravity_tidal_potential/screen-output b/tests/gravity_tidal_potential/screen-output index 2779dd879b0..2c5db84ada0 100644 --- a/tests/gravity_tidal_potential/screen-output +++ b/tests/gravity_tidal_potential/screen-output @@ -1,15 +1,10 @@ ------------------------------------------------------------------------------ ------------------------------------------------------------------------------ ------------------------------------------------------------------------------ ------------------------------------------------------------------------------ ------------------------------------------------------------------------------ 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)... 110+0 iterations. + Solving Stokes system (GMG)... 95+0 iterations. Postprocessing: RMS, max velocity: 2.04e-05 m/year, 7.16e-05 m/year @@ -131,9 +126,4 @@ Number of degrees of freedom: 4,030 (2,910+150+970) Termination requested by criterion: end time -+----------------------------------------------+------------+------------+ -+----------------------------------+-----------+------------+------------+ -+----------------------------------+-----------+------------+------------+ ------------------------------------------------------------------------------ ------------------------------------------------------------------------------ diff --git a/tests/gravity_tidal_potential/statistics b/tests/gravity_tidal_potential/statistics index b7210a2a72a..0e857453bdd 100644 --- a/tests/gravity_tidal_potential/statistics +++ b/tests/gravity_tidal_potential/statistics @@ -20,14 +20,14 @@ # 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 110 112 112 2.04033429e-05 7.15757443e-05 1.00000000e+02 9.99999535e+01 1.00000000e+02 output-gravity_tidal_potential/solution/solution-00000 -5.22786493e+05 5.82904414e+07 1.20123467e+08 9.17000000e+02 1.00000000e+20 2.63118556e+21 - 1 1.000000000000e+03 1.000000000000e+03 96 3060 970 0 50 51 51 2.53870595e-05 8.60011863e-05 1.00000000e+02 9.99999535e+01 1.00000000e+02 output-gravity_tidal_potential/solution/solution-00001 -5.23573285e+05 5.82904414e+07 1.20105088e+08 9.17000000e+02 1.00000000e+20 2.63118556e+21 - 2 2.000000000000e+03 1.000000000000e+03 96 3060 970 0 50 51 51 2.21503532e-05 8.26962512e-05 1.00000000e+02 9.99999535e+01 1.00000000e+02 output-gravity_tidal_potential/solution/solution-00002 -5.22786489e+05 5.82904414e+07 1.20011855e+08 9.17000000e+02 1.00000000e+20 2.63118556e+21 - 3 3.000000000000e+03 1.000000000000e+03 96 3060 970 0 50 51 51 2.28043260e-05 8.33561483e-05 1.00000000e+02 9.99999535e+01 1.00000000e+02 output-gravity_tidal_potential/solution/solution-00003 -5.23892834e+05 5.82904414e+07 1.19978816e+08 9.17000000e+02 1.00000000e+20 2.63118556e+21 - 4 4.000000000000e+03 1.000000000000e+03 96 3060 970 0 50 51 51 2.49936170e-05 8.22470390e-05 1.00000000e+02 9.99999535e+01 1.00000000e+02 output-gravity_tidal_potential/solution/solution-00004 -5.22786489e+05 5.82904414e+07 1.20075508e+08 9.17000000e+02 1.00000000e+20 2.63118556e+21 - 5 5.000000000000e+03 1.000000000000e+03 96 3060 970 0 50 51 51 2.04959459e-05 7.43142701e-05 1.00000000e+02 9.99999535e+01 1.00000000e+02 output-gravity_tidal_potential/solution/solution-00005 -5.22786489e+05 5.82904414e+07 1.20118224e+08 9.17000000e+02 1.00000000e+20 2.63118556e+21 - 6 6.000000000000e+03 1.000000000000e+03 96 3060 970 0 49 50 50 2.56511938e-05 8.84619423e-05 1.00000000e+02 9.99999535e+01 1.00000000e+02 output-gravity_tidal_potential/solution/solution-00006 -5.28450605e+05 5.82904414e+07 1.20124477e+08 9.17000000e+02 1.00000000e+20 2.63118556e+21 - 7 7.000000000000e+03 1.000000000000e+03 96 3060 970 0 50 51 51 2.15471702e-05 8.12170166e-05 1.00000000e+02 9.99999535e+01 1.00000000e+02 output-gravity_tidal_potential/solution/solution-00007 -5.22786490e+05 5.82904414e+07 1.20053747e+08 9.17000000e+02 1.00000000e+20 2.63118556e+21 - 8 8.000000000000e+03 1.000000000000e+03 96 3060 970 0 50 51 51 2.34666825e-05 8.31817499e-05 1.00000000e+02 9.99999535e+01 1.00000000e+02 output-gravity_tidal_potential/solution/solution-00008 -5.22786490e+05 5.82904414e+07 1.19974416e+08 9.17000000e+02 1.00000000e+20 2.63118556e+21 - 9 9.000000000000e+03 1.000000000000e+03 96 3060 970 0 50 51 51 2.44893181e-05 8.10214829e-05 1.00000000e+02 9.99999535e+01 1.00000000e+02 output-gravity_tidal_potential/solution/solution-00009 -5.22786489e+05 5.82904414e+07 1.20036188e+08 9.17000000e+02 1.00000000e+20 2.63118556e+21 -10 1.000000000000e+04 1.000000000000e+03 96 3060 970 0 50 51 51 2.07657183e-05 7.69853998e-05 1.00000000e+02 9.99999535e+01 1.00000000e+02 output-gravity_tidal_potential/solution/solution-00010 -5.22786489e+05 5.82904414e+07 1.20102575e+08 9.17000000e+02 1.00000000e+20 2.63118556e+21 + 0 0.000000000000e+00 0.000000000000e+00 96 3060 970 0 95 96 96 2.04033436e-05 7.15759315e-05 1.00000000e+02 9.99999535e+01 1.00000000e+02 output-gravity_tidal_potential/solution/solution-00000 -5.22786377e+05 5.82904414e+07 1.20123467e+08 9.17000000e+02 1.00000000e+20 2.63118556e+21 + 1 1.000000000000e+03 1.000000000000e+03 96 3060 970 0 50 51 51 2.53870602e-05 8.60009790e-05 1.00000000e+02 9.99999535e+01 1.00000000e+02 output-gravity_tidal_potential/solution/solution-00001 -5.23573278e+05 5.82904414e+07 1.20105088e+08 9.17000000e+02 1.00000000e+20 2.63118556e+21 + 2 2.000000000000e+03 1.000000000000e+03 96 3060 970 0 50 51 51 2.21503531e-05 8.26962710e-05 1.00000000e+02 9.99999535e+01 1.00000000e+02 output-gravity_tidal_potential/solution/solution-00002 -5.22786367e+05 5.82904414e+07 1.20011855e+08 9.17000000e+02 1.00000000e+20 2.63118556e+21 + 3 3.000000000000e+03 1.000000000000e+03 96 3060 970 0 50 51 51 2.28043259e-05 8.33561663e-05 1.00000000e+02 9.99999535e+01 1.00000000e+02 output-gravity_tidal_potential/solution/solution-00003 -5.23892832e+05 5.82904414e+07 1.19978816e+08 9.17000000e+02 1.00000000e+20 2.63118556e+21 + 4 4.000000000000e+03 1.000000000000e+03 96 3060 970 0 50 51 51 2.49936171e-05 8.22468386e-05 1.00000000e+02 9.99999535e+01 1.00000000e+02 output-gravity_tidal_potential/solution/solution-00004 -5.22786367e+05 5.82904414e+07 1.20075507e+08 9.17000000e+02 1.00000000e+20 2.63118556e+21 + 5 5.000000000000e+03 1.000000000000e+03 96 3060 970 0 50 51 51 2.04959460e-05 7.43144417e-05 1.00000000e+02 9.99999535e+01 1.00000000e+02 output-gravity_tidal_potential/solution/solution-00005 -5.22786369e+05 5.82904414e+07 1.20118224e+08 9.17000000e+02 1.00000000e+20 2.63118556e+21 + 6 6.000000000000e+03 1.000000000000e+03 96 3060 970 0 49 50 50 2.56511932e-05 8.84617474e-05 1.00000000e+02 9.99999535e+01 1.00000000e+02 output-gravity_tidal_potential/solution/solution-00006 -5.28450608e+05 5.82904414e+07 1.20124477e+08 9.17000000e+02 1.00000000e+20 2.63118556e+21 + 7 7.000000000000e+03 1.000000000000e+03 96 3060 970 0 50 51 51 2.15471700e-05 8.12170330e-05 1.00000000e+02 9.99999535e+01 1.00000000e+02 output-gravity_tidal_potential/solution/solution-00007 -5.22786369e+05 5.82904414e+07 1.20053747e+08 9.17000000e+02 1.00000000e+20 2.63118556e+21 + 8 8.000000000000e+03 1.000000000000e+03 96 3060 970 0 50 51 51 2.34666825e-05 8.31817714e-05 1.00000000e+02 9.99999535e+01 1.00000000e+02 output-gravity_tidal_potential/solution/solution-00008 -5.22786369e+05 5.82904414e+07 1.19974416e+08 9.17000000e+02 1.00000000e+20 2.63118556e+21 + 9 9.000000000000e+03 1.000000000000e+03 96 3060 970 0 50 51 51 2.44893184e-05 8.10214991e-05 1.00000000e+02 9.99999535e+01 1.00000000e+02 output-gravity_tidal_potential/solution/solution-00009 -5.22786370e+05 5.82904414e+07 1.20036187e+08 9.17000000e+02 1.00000000e+20 2.63118556e+21 +10 1.000000000000e+04 1.000000000000e+03 96 3060 970 0 50 51 51 2.07657185e-05 7.69854185e-05 1.00000000e+02 9.99999535e+01 1.00000000e+02 output-gravity_tidal_potential/solution/solution-00010 -5.22786370e+05 5.82904414e+07 1.20102575e+08 9.17000000e+02 1.00000000e+20 2.63118556e+21 From 98784e42763bd67bc6fd045d2c681547a9e60379 Mon Sep 17 00:00:00 2001 From: Hyunseong Kim Date: Tue, 16 Sep 2025 16:21:57 +0100 Subject: [PATCH 05/21] adding contributors to change log --- doc/modules/changes/20250716_hyunseong96 | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/doc/modules/changes/20250716_hyunseong96 b/doc/modules/changes/20250716_hyunseong96 index c1a0f9bebce..f60eb1ba058 100644 --- a/doc/modules/changes/20250716_hyunseong96 +++ b/doc/modules/changes/20250716_hyunseong96 @@ -1,4 +1,4 @@ -Added: ASPECT now has a new gravity model plugin called 'Radial linear with tidal potnetial' for tidal force. +Added: ASPECT now has a new gravity model plugin called 'Radial linear with tidal potnetial' for tidal forces. This plugin is useful for modeling long-term interior and surface evolution in moons orbiting a large planet.
-(Hyunseong Kim, 2025/07/16) +(Hyunseong Kim, Antoniette Greta Grima, Wolfgang Bangerth 2025/07/16) From d814c23f88a8def3e08d18475a3d87107b0aec7a Mon Sep 17 00:00:00 2001 From: Hyunseong Kim Date: Tue, 16 Sep 2025 23:34:55 +0100 Subject: [PATCH 06/21] changing time parameter name --- tests/gravity_tidal_potential.prm | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/gravity_tidal_potential.prm b/tests/gravity_tidal_potential.prm index 491312ef187..e3e9de82f94 100644 --- a/tests/gravity_tidal_potential.prm +++ b/tests/gravity_tidal_potential.prm @@ -13,7 +13,7 @@ # Visualization 'gravity' shows gravity distribution. set Dimension = 3 -set Use years in output instead of seconds = true +set Use years instead of seconds = true set End time = 1e4 set Output directory = gravity_tidal_potential From 8a8e39f7f43f96c26371b5a166cd6e8291ae36ea Mon Sep 17 00:00:00 2001 From: Hyunseong Kim Date: Wed, 17 Sep 2025 15:58:26 +0100 Subject: [PATCH 07/21] applying test diff --- tests/gravity_tidal_potential/statistics | 22 +++++++++++----------- 1 file changed, 11 insertions(+), 11 deletions(-) diff --git a/tests/gravity_tidal_potential/statistics b/tests/gravity_tidal_potential/statistics index 0e857453bdd..02ef831ab03 100644 --- a/tests/gravity_tidal_potential/statistics +++ b/tests/gravity_tidal_potential/statistics @@ -20,14 +20,14 @@ # 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 95 96 96 2.04033436e-05 7.15759315e-05 1.00000000e+02 9.99999535e+01 1.00000000e+02 output-gravity_tidal_potential/solution/solution-00000 -5.22786377e+05 5.82904414e+07 1.20123467e+08 9.17000000e+02 1.00000000e+20 2.63118556e+21 - 1 1.000000000000e+03 1.000000000000e+03 96 3060 970 0 50 51 51 2.53870602e-05 8.60009790e-05 1.00000000e+02 9.99999535e+01 1.00000000e+02 output-gravity_tidal_potential/solution/solution-00001 -5.23573278e+05 5.82904414e+07 1.20105088e+08 9.17000000e+02 1.00000000e+20 2.63118556e+21 - 2 2.000000000000e+03 1.000000000000e+03 96 3060 970 0 50 51 51 2.21503531e-05 8.26962710e-05 1.00000000e+02 9.99999535e+01 1.00000000e+02 output-gravity_tidal_potential/solution/solution-00002 -5.22786367e+05 5.82904414e+07 1.20011855e+08 9.17000000e+02 1.00000000e+20 2.63118556e+21 - 3 3.000000000000e+03 1.000000000000e+03 96 3060 970 0 50 51 51 2.28043259e-05 8.33561663e-05 1.00000000e+02 9.99999535e+01 1.00000000e+02 output-gravity_tidal_potential/solution/solution-00003 -5.23892832e+05 5.82904414e+07 1.19978816e+08 9.17000000e+02 1.00000000e+20 2.63118556e+21 - 4 4.000000000000e+03 1.000000000000e+03 96 3060 970 0 50 51 51 2.49936171e-05 8.22468386e-05 1.00000000e+02 9.99999535e+01 1.00000000e+02 output-gravity_tidal_potential/solution/solution-00004 -5.22786367e+05 5.82904414e+07 1.20075507e+08 9.17000000e+02 1.00000000e+20 2.63118556e+21 - 5 5.000000000000e+03 1.000000000000e+03 96 3060 970 0 50 51 51 2.04959460e-05 7.43144417e-05 1.00000000e+02 9.99999535e+01 1.00000000e+02 output-gravity_tidal_potential/solution/solution-00005 -5.22786369e+05 5.82904414e+07 1.20118224e+08 9.17000000e+02 1.00000000e+20 2.63118556e+21 - 6 6.000000000000e+03 1.000000000000e+03 96 3060 970 0 49 50 50 2.56511932e-05 8.84617474e-05 1.00000000e+02 9.99999535e+01 1.00000000e+02 output-gravity_tidal_potential/solution/solution-00006 -5.28450608e+05 5.82904414e+07 1.20124477e+08 9.17000000e+02 1.00000000e+20 2.63118556e+21 - 7 7.000000000000e+03 1.000000000000e+03 96 3060 970 0 50 51 51 2.15471700e-05 8.12170330e-05 1.00000000e+02 9.99999535e+01 1.00000000e+02 output-gravity_tidal_potential/solution/solution-00007 -5.22786369e+05 5.82904414e+07 1.20053747e+08 9.17000000e+02 1.00000000e+20 2.63118556e+21 - 8 8.000000000000e+03 1.000000000000e+03 96 3060 970 0 50 51 51 2.34666825e-05 8.31817714e-05 1.00000000e+02 9.99999535e+01 1.00000000e+02 output-gravity_tidal_potential/solution/solution-00008 -5.22786369e+05 5.82904414e+07 1.19974416e+08 9.17000000e+02 1.00000000e+20 2.63118556e+21 - 9 9.000000000000e+03 1.000000000000e+03 96 3060 970 0 50 51 51 2.44893184e-05 8.10214991e-05 1.00000000e+02 9.99999535e+01 1.00000000e+02 output-gravity_tidal_potential/solution/solution-00009 -5.22786370e+05 5.82904414e+07 1.20036187e+08 9.17000000e+02 1.00000000e+20 2.63118556e+21 -10 1.000000000000e+04 1.000000000000e+03 96 3060 970 0 50 51 51 2.07657185e-05 7.69854185e-05 1.00000000e+02 9.99999535e+01 1.00000000e+02 output-gravity_tidal_potential/solution/solution-00010 -5.22786370e+05 5.82904414e+07 1.20102575e+08 9.17000000e+02 1.00000000e+20 2.63118556e+21 + 0 0.000000000000e+00 0.000000000000e+00 96 3060 970 0 95 96 96 2.04033436e-05 7.15759491e-05 1.00000000e+02 9.99999535e+01 1.00000000e+02 output-gravity_tidal_potential/solution/solution-00000 -5.22786368e+05 5.82904414e+07 1.20123467e+08 9.17000000e+02 1.00000000e+20 2.63118556e+21 + 1 1.000000000000e+03 1.000000000000e+03 96 3060 970 0 50 51 51 2.53870602e-05 8.60009600e-05 1.00000000e+02 9.99999535e+01 1.00000000e+02 output-gravity_tidal_potential/solution/solution-00001 -5.23573277e+05 5.82904414e+07 1.20105088e+08 9.17000000e+02 1.00000000e+20 2.63118556e+21 + 2 2.000000000000e+03 1.000000000000e+03 96 3060 970 0 50 51 51 2.21503531e-05 8.26962723e-05 1.00000000e+02 9.99999535e+01 1.00000000e+02 output-gravity_tidal_potential/solution/solution-00002 -5.22786356e+05 5.82904414e+07 1.20011855e+08 9.17000000e+02 1.00000000e+20 2.63118556e+21 + 3 3.000000000000e+03 1.000000000000e+03 96 3060 970 0 50 51 51 2.28043259e-05 8.33561683e-05 1.00000000e+02 9.99999535e+01 1.00000000e+02 output-gravity_tidal_potential/solution/solution-00003 -5.23892832e+05 5.82904414e+07 1.19978816e+08 9.17000000e+02 1.00000000e+20 2.63118556e+21 + 4 4.000000000000e+03 1.000000000000e+03 96 3060 970 0 50 51 51 2.49936171e-05 8.22468188e-05 1.00000000e+02 9.99999535e+01 1.00000000e+02 output-gravity_tidal_potential/solution/solution-00004 -5.22786356e+05 5.82904414e+07 1.20075507e+08 9.17000000e+02 1.00000000e+20 2.63118556e+21 + 5 5.000000000000e+03 1.000000000000e+03 96 3060 970 0 50 51 51 2.04959462e-05 7.43144594e-05 1.00000000e+02 9.99999535e+01 1.00000000e+02 output-gravity_tidal_potential/solution/solution-00005 -5.22786357e+05 5.82904414e+07 1.20118224e+08 9.17000000e+02 1.00000000e+20 2.63118556e+21 + 6 6.000000000000e+03 1.000000000000e+03 96 3060 970 0 49 50 50 2.56511932e-05 8.84617285e-05 1.00000000e+02 9.99999535e+01 1.00000000e+02 output-gravity_tidal_potential/solution/solution-00006 -5.28450608e+05 5.82904414e+07 1.20124477e+08 9.17000000e+02 1.00000000e+20 2.63118556e+21 + 7 7.000000000000e+03 1.000000000000e+03 96 3060 970 0 50 51 51 2.15471701e-05 8.12170364e-05 1.00000000e+02 9.99999535e+01 1.00000000e+02 output-gravity_tidal_potential/solution/solution-00007 -5.22786357e+05 5.82904414e+07 1.20053747e+08 9.17000000e+02 1.00000000e+20 2.63118556e+21 + 8 8.000000000000e+03 1.000000000000e+03 96 3060 970 0 50 51 51 2.34666825e-05 8.31817721e-05 1.00000000e+02 9.99999535e+01 1.00000000e+02 output-gravity_tidal_potential/solution/solution-00008 -5.22786357e+05 5.82904414e+07 1.19974416e+08 9.17000000e+02 1.00000000e+20 2.63118556e+21 + 9 9.000000000000e+03 1.000000000000e+03 96 3060 970 0 50 51 51 2.44893181e-05 8.10215036e-05 1.00000000e+02 9.99999535e+01 1.00000000e+02 output-gravity_tidal_potential/solution/solution-00009 -5.22786359e+05 5.82904414e+07 1.20036187e+08 9.17000000e+02 1.00000000e+20 2.63118556e+21 +10 1.000000000000e+04 1.000000000000e+03 96 3060 970 0 50 51 51 2.07657184e-05 7.69854194e-05 1.00000000e+02 9.99999535e+01 1.00000000e+02 output-gravity_tidal_potential/solution/solution-00010 -5.22786358e+05 5.82904414e+07 1.20102575e+08 9.17000000e+02 1.00000000e+20 2.63118556e+21 From 66a6297c9b8e73bc176b14fef27954227d768285 Mon Sep 17 00:00:00 2001 From: Hyunseong Kim Date: Thu, 25 Sep 2025 15:31:52 +0100 Subject: [PATCH 08/21] modifying by suggestions --- include/aspect/gravity_model/radial.h | 48 +---- .../radial_with_tidal_potential.h | 101 +++++++++ source/gravity_model/radial.cc | 165 --------------- .../radial_with_tidal_potential.cc | 199 ++++++++++++++++++ tests/gravity_tidal_potential.prm | 5 +- 5 files changed, 304 insertions(+), 214 deletions(-) create mode 100644 include/aspect/gravity_model/radial_with_tidal_potential.h create mode 100644 source/gravity_model/radial_with_tidal_potential.cc diff --git a/include/aspect/gravity_model/radial.h b/include/aspect/gravity_model/radial.h index 3519f21b268..0a94d0eb54e 100644 --- a/include/aspect/gravity_model/radial.h +++ b/include/aspect/gravity_model/radial.h @@ -62,6 +62,8 @@ namespace aspect * Magnitude of the gravity vector. */ double magnitude; + + template friend class RadialWithTidalPotential; }; @@ -135,52 +137,6 @@ namespace aspect double magnitude_at_bottom; }; - - - - template - class RadialLinearWithTidalPotential : public RadialLinear - { - 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 hosting planet - */ - double M_p; - - /** - * Semimajor axis of satellite's orbit - */ - double a_s; - - /** - * Rate of non-synchronous rotation in m/year - */ - double b_NSR; - - /** - * Radius of modelling satellite - */ - double R1; - }; } } 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..cebbce42da9 --- /dev/null +++ b/include/aspect/gravity_model/radial_with_tidal_potential.h @@ -0,0 +1,101 @@ +/* + 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 by tidal potential from flattening and non-synchronnous rotation with respect to 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)). + * Tidal potential is positive because it is 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 planet, 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. + * b [1/s] = b_NSR [m/yr] / (circumference of satellite [m]) / year_to_seconds [s/yr] * 2 * pi + * + * @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: + /** + * Declare a member variable of type RadialConstant that allows us to call + * functions from radial_constant.cc. + */ + RadialConstant radialconstant; + + /** + * Mass of the perturbing body + */ + double M_p; + + /** + * Semimajor axis of the modeled body + */ + double a_s; + + /** + * Rate of the non-synchronous rotation in m/year + */ + double b_NSR; + + /** + * Radius of the modeled body + */ + double R1; + }; + } +} + +#endif diff --git a/source/gravity_model/radial.cc b/source/gravity_model/radial.cc index c0a0b85c0b4..c57dc1d9328 100644 --- a/source/gravity_model/radial.cc +++ b/source/gravity_model/radial.cc @@ -183,159 +183,6 @@ namespace aspect "do not have either a spherical or ellipsoidal natural coordinate system.")); } - - - -// ----------------------------- RadialLinearWithTidalPotential ---------------------- - - - template - Tensor<1,dim> - RadialLinearWithTidalPotential::gravity_vector (const Point &p) const - { - /** - * This gravity model is a 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) / (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, 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. - * b [1/s] = b_NSR [m/yr] / (circumference of satellite [m]) / year_to_seconds [s/yr] * 2 * pi - * - * 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 Tensor<1,dim> e_x = Point::unit_vector(0); - const Tensor<1,dim> e_y = Point::unit_vector(1); - const Tensor<1,dim> e_z = Point::unit_vector(2); - - const double t = (this->simulator_is_past_initialization()) ? this->get_time() : 0.0; - const double x = p[0]; - const double y = p[1]; - const double z = p[2]; - - const double C1 = std::cos( 2. * b_NSR / R1 / year_in_seconds * t); - const double C2 = std::sin( 2. * b_NSR / R1 / year_in_seconds * t); - - const double dTstar_over_dx = 1. / 6. * ( 2. * x ); - const double dT0_over_dx = 1. / 2. * ( 2. * C1 * x - 2. * C2 * y ); - - const double dTstar_over_dy = 1. / 6. * ( 2. * y ); - const double dT0_over_dy = 1. / 2. * ( -2. * C1 * y - 2. * C2 * x ); - - const double dTstar_over_dz = 1. / 6. * ( -4. * z ); - const double dT0_over_dz = 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_over_dx + dT0_over_dx) * e_x - + (dTstar_over_dy + dT0_over_dy) * e_y - + (dTstar_over_dz + dT0_over_dz) * e_z); - - return RadialLinear::gravity_vector(p) + tidal_gravity; - } - - - template - void - RadialLinearWithTidalPotential::declare_parameters (ParameterHandler &prm) - { - RadialLinear::declare_parameters (prm); - - prm.enter_subsection("Gravity model"); - { - prm.enter_subsection("Radial linear with tidal potnetial"); - { - prm.declare_entry ("Mass of planet", "1.898e27", - Patterns::Double (), - "Mass of satellte's host planet. " - "Default value is for Europa, therefore, mass of Jupiter. " - "Units is $kg$."); - prm.declare_entry ("Semimajor axis of satellite", "670900000", - Patterns::Double (), - "Length of semimajor axis of satellite. " - "Default value is for Europa's semimajor axis" - "Units is $m$."); - prm.declare_entry ("Rate of nonsynchronous rotation", "1000", - Patterns::Double (), - "Rate of nonsynchronous rotation (NSR). " - "This will be converted to angular rate. " - "Units is $m/year$"); - } - prm.leave_subsection (); - } - prm.leave_subsection (); - } - - - template - void - RadialLinearWithTidalPotential::parse_parameters (ParameterHandler &prm) - { - AssertThrow (dim==3, ExcMessage ("The 'radial linear with tidal potential' gravity model " - "can only be used in 3D.")); - this->RadialLinear::parse_parameters (prm); - - prm.enter_subsection("Gravity model"); - { - prm.enter_subsection("Radial linear with tidal potnetial"); - { - M_p = prm.get_double ("Mass of planet"); - a_s = prm.get_double ("Semimajor axis of satellite"); - b_NSR = prm.get_double ("Rate of nonsynchronous rotation"); - } - prm.leave_subsection (); - } - prm.leave_subsection (); - // This effect of tidal potential only works if the geometry is derived from - // a spherical model (i.e. a sphere, spherical shell or chunk) - if (Plugins::plugin_type_matches>(this->get_geometry_model())) - { - R1 = Plugins::get_plugin_as_type> - (this->get_geometry_model()).radius(); - } - else if (Plugins::plugin_type_matches>(this->get_geometry_model())) - { - R1 = Plugins::get_plugin_as_type> - (this->get_geometry_model()).outer_radius(); - } - else if (Plugins::plugin_type_matches>(this->get_geometry_model())) - { - R1 = Plugins::get_plugin_as_type> - (this->get_geometry_model()).outer_radius(); - } - else if (Plugins::plugin_type_matches>(this->get_geometry_model())) - { - const auto &gm = Plugins::get_plugin_as_type> - (this->get_geometry_model()); - // TODO - // If the eccentricity of the EllipsoidalChunk is non-zero, the radius can vary along a boundary, - // but the maximal depth is the same everywhere and we could calculate a representative pressure - // profile. However, it requires some extra logic with ellipsoidal - // coordinates, so for now we only allow eccentricity zero. - // Using the EllipsoidalChunk with eccentricity zero can still be useful, - // because the domain can be non-coordinate parallel. - - AssertThrow(gm.get_eccentricity() == 0.0, - ExcNotImplemented("This plugin cannot be used with a non-zero eccentricity. ")); - - R1 = gm.get_semi_major_axis_a(); - } - else - { - Assert (false, ExcMessage ("This initial condition can only be used if the geometry " - "is a sphere, a spherical shell, a chunk or an " - "ellipsoidal chunk.")); - R1 = numbers::signaling_nan(); - } - } } } @@ -365,17 +212,5 @@ namespace aspect "is positive) and the magnitude changes linearly with depth. The " "magnitude of gravity at the surface and bottom is read from the " "input file in a section ``Gravity model/Radial linear''.") - - ASPECT_REGISTER_GRAVITY_MODEL(RadialLinearWithTidalPotential, - "radial linear with tidal potential", - "A gravity model that is the sum of the `radial linear' model " - "(which is radial, pointing inward if the gravity " - "is positive, and a magnitude that changes linearly with depth), " - "and a term that results from a tidal gravity potential and that " - "leads to a gravity field that varies with latitude and longitude." - "The magnitude of gravity for the linear radial part is read from the " - "input file in a section `Gravity model/Radial linear'; the " - "parameters that describe the tidal potential contribution are read " - "from a section ``Gravity model/Radial linear with tidal potential''.") } } 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..dd1a309e401 --- /dev/null +++ b/source/gravity_model/radial_with_tidal_potential.cc @@ -0,0 +1,199 @@ +/* + 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 +#include +#include +#include +#include + +namespace aspect +{ + namespace GravityModel + { + template + Tensor<1,dim> + RadialWithTidalPotential::gravity_vector (const Point &p) const + { + /** + * 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 Tensor<1,dim> e_x = Point::unit_vector(0); + const Tensor<1,dim> e_y = Point::unit_vector(1); + const Tensor<1,dim> e_z = Point::unit_vector(2); + + const double t = (this->simulator_is_past_initialization()) ? this->get_time() : 0.0; + const double x = p[0]; + const double y = p[1]; + const double z = p[2]; + + const double C1 = std::cos( 2. * b_NSR / R1 / year_in_seconds * t); + const double C2 = std::sin( 2. * b_NSR / R1 / year_in_seconds * t); + + const double dTstar_over_dx = 1. / 6. * ( 2. * x ); + const double dT0_over_dx = 1. / 2. * ( 2. * C1 * x - 2. * C2 * y ); + + const double dTstar_over_dy = 1. / 6. * ( 2. * y ); + const double dT0_over_dy = 1. / 2. * ( -2. * C1 * y - 2. * C2 * x ); + + const double dTstar_over_dz = 1. / 6. * ( -4. * z ); + const double dT0_over_dz = 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_over_dx + dT0_over_dx) * e_x + + (dTstar_over_dy + dT0_over_dy) * e_y + + (dTstar_over_dz + dT0_over_dz) * e_z); + + if (p.norm() == 0.0) + return Tensor<1,dim>(); + + const double r = p.norm(); + return -radialconstant.magnitude * p/r + tidal_gravity; + } + + + template + void + RadialWithTidalPotential::declare_parameters (ParameterHandler &prm) + { + prm.enter_subsection("Gravity model"); + { + prm.enter_subsection("Radial with tidal potential"); + { + prm.declare_entry ("Mass of planet", "1.898e27", + Patterns::Double (), + "Mass of satellte's host planet. " + "Default value is for Europa, therefore, mass of Jupiter. " + "Units is $kg$."); + prm.declare_entry ("Semimajor axis of satellite", "670900000", + Patterns::Double (), + "Length of semimajor axis of satellite. " + "Default value is for Europa's semimajor axis" + "Units is $m$."); + prm.declare_entry ("Rate of nonsynchronous rotation", "1000", + Patterns::Double (), + "Rate of nonsynchronous rotation (NSR). " + "This will be converted to angular rate. " + "Units is $m/year$"); + } + 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 planet"); + a_s = prm.get_double ("Semimajor axis of satellite"); + b_NSR = prm.get_double ("Rate of nonsynchronous rotation"); + } + prm.leave_subsection (); + } + prm.leave_subsection (); + + radialconstant.initialize_simulator (this->get_simulator()); + radialconstant.parse_parameters(prm); + radialconstant.initialize(); + + // This effect of tidal potential only works if the geometry is derived from + // a spherical model (i.e. a sphere, spherical shell or chunk) + if (Plugins::plugin_type_matches>(this->get_geometry_model())) + { + R1 = Plugins::get_plugin_as_type> + (this->get_geometry_model()).radius(); + } + else if (Plugins::plugin_type_matches>(this->get_geometry_model())) + { + R1 = Plugins::get_plugin_as_type> + (this->get_geometry_model()).outer_radius(); + } + else if (Plugins::plugin_type_matches>(this->get_geometry_model())) + { + R1 = Plugins::get_plugin_as_type> + (this->get_geometry_model()).outer_radius(); + } + else if (Plugins::plugin_type_matches>(this->get_geometry_model())) + { + const auto &gm = Plugins::get_plugin_as_type> + (this->get_geometry_model()); + // TODO + // If the eccentricity of the EllipsoidalChunk is non-zero, the radius can vary along a boundary, + // but the maximal depth is the same everywhere and we could calculate a representative pressure + // profile. However, it requires some extra logic with ellipsoidal + // coordinates, so for now we only allow eccentricity zero. + // Using the EllipsoidalChunk with eccentricity zero can still be useful, + // because the domain can be non-coordinate parallel. + + AssertThrow(gm.get_eccentricity() == 0.0, + ExcNotImplemented("This plugin cannot be used with a non-zero eccentricity. ")); + + R1 = gm.get_semi_major_axis_a(); + } + else + { + Assert (false, ExcMessage ("This initial condition can only be used if the geometry " + "is a sphere, a spherical shell, a chunk or an " + "ellipsoidal chunk.")); + R1 = numbers::signaling_nan(); + } + } + } +} + +// 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 linear with tidal potential''.") + } +} diff --git a/tests/gravity_tidal_potential.prm b/tests/gravity_tidal_potential.prm index e3e9de82f94..565d49d470f 100644 --- a/tests/gravity_tidal_potential.prm +++ b/tests/gravity_tidal_potential.prm @@ -54,9 +54,8 @@ end subsection Gravity model set Model name = radial linear with tidal potential - subsection Radial linear - set Magnitude at surface = 1.3 - set Magnitude at bottom = 1.3 + subsection Radial constant + set Magnitude = 1.3 end subsection Radial linear with tidal potnetial end From 545332f95e3b9e195305cc29876a6c5697b1f34f Mon Sep 17 00:00:00 2001 From: Hyunseong Kim Date: Thu, 25 Sep 2025 18:14:22 +0100 Subject: [PATCH 09/21] minor corrections --- .../gravity_model/radial_with_tidal_potential.h | 2 +- .../gravity_model/radial_with_tidal_potential.cc | 16 +++++++++------- 2 files changed, 10 insertions(+), 8 deletions(-) diff --git a/include/aspect/gravity_model/radial_with_tidal_potential.h b/include/aspect/gravity_model/radial_with_tidal_potential.h index cebbce42da9..b794ed6b8a8 100644 --- a/include/aspect/gravity_model/radial_with_tidal_potential.h +++ b/include/aspect/gravity_model/radial_with_tidal_potential.h @@ -34,7 +34,7 @@ namespace aspect * A class that describes gravity as a radial vector of linearly * changing magnitude by tidal potential from flattening and non-synchronnous rotation with respect to 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), + * 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 it is from geodesy research, where potential is taken as positive. diff --git a/source/gravity_model/radial_with_tidal_potential.cc b/source/gravity_model/radial_with_tidal_potential.cc index dd1a309e401..33d5337ea15 100644 --- a/source/gravity_model/radial_with_tidal_potential.cc +++ b/source/gravity_model/radial_with_tidal_potential.cc @@ -89,19 +89,21 @@ namespace aspect { prm.enter_subsection("Radial with tidal potential"); { - prm.declare_entry ("Mass of planet", "1.898e27", + prm.declare_entry ("Mass of perterbing body", "1.898e27", Patterns::Double (), - "Mass of satellte's host planet. " - "Default value is for Europa, therefore, mass of Jupiter. " + "Mass of body that perturbs gravity of modeled body. " + "Default value is for modeling Europa, therefore, mass of Jupiter. " "Units is $kg$."); - prm.declare_entry ("Semimajor axis of satellite", "670900000", + prm.declare_entry ("Semimajor axis of orbit", "670900000", Patterns::Double (), - "Length of semimajor axis of satellite. " + "Length of semimajor axis of orbit between modeled body and perturbing body. " "Default value is for Europa's semimajor axis" "Units is $m$."); prm.declare_entry ("Rate of nonsynchronous rotation", "1000", Patterns::Double (), "Rate of nonsynchronous rotation (NSR). " + "This works for the modeled body having decoupled rotation between interior layers. " + "Default value is for Europa's icy shell. " "This will be converted to angular rate. " "Units is $m/year$"); } @@ -122,8 +124,8 @@ namespace aspect { prm.enter_subsection("Radial with tidal potential"); { - M_p = prm.get_double ("Mass of planet"); - a_s = prm.get_double ("Semimajor axis of satellite"); + M_p = prm.get_double ("Mass of perturbing body"); + a_s = prm.get_double ("Semimajor axis of orbit"); b_NSR = prm.get_double ("Rate of nonsynchronous rotation"); } prm.leave_subsection (); From 58b8278860d138b2549d23ddc070bf4d85a8ca86 Mon Sep 17 00:00:00 2001 From: Hyunseong Kim Date: Fri, 26 Sep 2025 00:21:54 +0100 Subject: [PATCH 10/21] minor corrections about plugin --- source/gravity_model/radial_with_tidal_potential.cc | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/source/gravity_model/radial_with_tidal_potential.cc b/source/gravity_model/radial_with_tidal_potential.cc index 33d5337ea15..ebb124b1480 100644 --- a/source/gravity_model/radial_with_tidal_potential.cc +++ b/source/gravity_model/radial_with_tidal_potential.cc @@ -196,6 +196,6 @@ namespace aspect "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 linear with tidal potential''.") + "from a section ``Gravity model/Radial with tidal potential''.") } } From 19db0c757e50ad18fb19025f1165d2bf668196e0 Mon Sep 17 00:00:00 2001 From: Hyunseong Kim Date: Wed, 1 Oct 2025 18:41:40 +0100 Subject: [PATCH 11/21] updating angular rate --- .../gravity_model/radial_with_tidal_potential.h | 7 +------ .../gravity_model/radial_with_tidal_potential.cc | 16 +++++++--------- 2 files changed, 8 insertions(+), 15 deletions(-) diff --git a/include/aspect/gravity_model/radial_with_tidal_potential.h b/include/aspect/gravity_model/radial_with_tidal_potential.h index b794ed6b8a8..cd82ef43f79 100644 --- a/include/aspect/gravity_model/radial_with_tidal_potential.h +++ b/include/aspect/gravity_model/radial_with_tidal_potential.h @@ -86,14 +86,9 @@ namespace aspect double a_s; /** - * Rate of the non-synchronous rotation in m/year + * Angular rate of the non-synchronous rotation in m/year */ double b_NSR; - - /** - * Radius of the modeled body - */ - double R1; }; } } diff --git a/source/gravity_model/radial_with_tidal_potential.cc b/source/gravity_model/radial_with_tidal_potential.cc index ebb124b1480..d4779f00d56 100644 --- a/source/gravity_model/radial_with_tidal_potential.cc +++ b/source/gravity_model/radial_with_tidal_potential.cc @@ -53,8 +53,8 @@ namespace aspect const double y = p[1]; const double z = p[2]; - const double C1 = std::cos( 2. * b_NSR / R1 / year_in_seconds * t); - const double C2 = std::sin( 2. * b_NSR / R1 / year_in_seconds * t); + const double C1 = std::cos( 2. * b_NSR * t); + const double C2 = std::sin( 2. * b_NSR * t); const double dTstar_over_dx = 1. / 6. * ( 2. * x ); const double dT0_over_dx = 1. / 2. * ( 2. * C1 * x - 2. * C2 * y ); @@ -99,13 +99,12 @@ namespace aspect "Length of semimajor axis of orbit between modeled body and perturbing body. " "Default value is for Europa's semimajor axis" "Units is $m$."); - prm.declare_entry ("Rate of nonsynchronous rotation", "1000", + prm.declare_entry ("Angular rate of nonsynchronous rotation", "1000", Patterns::Double (), - "Rate of nonsynchronous rotation (NSR). " + "Angular rate of nonsynchronous rotation (NSR). " "This works for the modeled body having decoupled rotation between interior layers. " "Default value is for Europa's icy shell. " - "This will be converted to angular rate. " - "Units is $m/year$"); + "Units is $degrees/year$"); } prm.leave_subsection (); } @@ -126,15 +125,14 @@ namespace aspect { M_p = prm.get_double ("Mass of perturbing body"); a_s = prm.get_double ("Semimajor axis of orbit"); - b_NSR = prm.get_double ("Rate of nonsynchronous rotation"); + const double time_scale = this->get_parameters().convert_to_years ? constants::year_in_seconds : 1.0; + b_NSR = prm.get_double ("Angular rate of nonsynchronous rotation") * constants::degree_to_radians / time_scale; } prm.leave_subsection (); } prm.leave_subsection (); - radialconstant.initialize_simulator (this->get_simulator()); radialconstant.parse_parameters(prm); - radialconstant.initialize(); // This effect of tidal potential only works if the geometry is derived from // a spherical model (i.e. a sphere, spherical shell or chunk) From 4478a2868171f7f94d13e27964e7c632dcb3a9ae Mon Sep 17 00:00:00 2001 From: Hyunseong Kim Date: Mon, 6 Oct 2025 15:02:45 +0100 Subject: [PATCH 12/21] default value change and the status working with friend class --- doc/modules/changes/20250716_hyunseong96 | 2 +- .../gravity_model/radial_with_tidal_potential.cc | 16 +++------------- 2 files changed, 4 insertions(+), 14 deletions(-) diff --git a/doc/modules/changes/20250716_hyunseong96 b/doc/modules/changes/20250716_hyunseong96 index f60eb1ba058..360e963c29f 100644 --- a/doc/modules/changes/20250716_hyunseong96 +++ b/doc/modules/changes/20250716_hyunseong96 @@ -1,4 +1,4 @@ -Added: ASPECT now has a new gravity model plugin called 'Radial linear with tidal potnetial' for tidal forces. +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/source/gravity_model/radial_with_tidal_potential.cc b/source/gravity_model/radial_with_tidal_potential.cc index d4779f00d56..e4e44cf9765 100644 --- a/source/gravity_model/radial_with_tidal_potential.cc +++ b/source/gravity_model/radial_with_tidal_potential.cc @@ -85,11 +85,12 @@ namespace aspect 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 perterbing body", "1.898e27", + prm.declare_entry ("Mass of perturbing body", "1.898e27", Patterns::Double (), "Mass of body that perturbs gravity of modeled body. " "Default value is for modeling Europa, therefore, mass of Jupiter. " @@ -99,7 +100,7 @@ namespace aspect "Length of semimajor axis of orbit between modeled body and perturbing body. " "Default value is for Europa's semimajor axis" "Units is $m$."); - prm.declare_entry ("Angular rate of nonsynchronous rotation", "1000", + 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. " @@ -131,25 +132,17 @@ namespace aspect prm.leave_subsection (); } prm.leave_subsection (); - - radialconstant.parse_parameters(prm); // This effect of tidal potential only works if the geometry is derived from // a spherical model (i.e. a sphere, spherical shell or chunk) if (Plugins::plugin_type_matches>(this->get_geometry_model())) { - R1 = Plugins::get_plugin_as_type> - (this->get_geometry_model()).radius(); } else if (Plugins::plugin_type_matches>(this->get_geometry_model())) { - R1 = Plugins::get_plugin_as_type> - (this->get_geometry_model()).outer_radius(); } else if (Plugins::plugin_type_matches>(this->get_geometry_model())) { - R1 = Plugins::get_plugin_as_type> - (this->get_geometry_model()).outer_radius(); } else if (Plugins::plugin_type_matches>(this->get_geometry_model())) { @@ -165,15 +158,12 @@ namespace aspect AssertThrow(gm.get_eccentricity() == 0.0, ExcNotImplemented("This plugin cannot be used with a non-zero eccentricity. ")); - - R1 = gm.get_semi_major_axis_a(); } else { Assert (false, ExcMessage ("This initial condition can only be used if the geometry " "is a sphere, a spherical shell, a chunk or an " "ellipsoidal chunk.")); - R1 = numbers::signaling_nan(); } } } From 2a0042b5598f9a6fd8c6f637122cc6a0ddfd1a84 Mon Sep 17 00:00:00 2001 From: Hyunseong Kim Date: Mon, 6 Oct 2025 21:16:53 +0100 Subject: [PATCH 13/21] working with public --- include/aspect/gravity_model/radial.h | 2 -- .../radial_with_tidal_potential.h | 19 ++++--------- .../radial_with_tidal_potential.cc | 28 +++++++++---------- 3 files changed, 19 insertions(+), 30 deletions(-) diff --git a/include/aspect/gravity_model/radial.h b/include/aspect/gravity_model/radial.h index 0a94d0eb54e..07059b5ecbc 100644 --- a/include/aspect/gravity_model/radial.h +++ b/include/aspect/gravity_model/radial.h @@ -62,8 +62,6 @@ namespace aspect * Magnitude of the gravity vector. */ double magnitude; - - template friend class RadialWithTidalPotential; }; diff --git a/include/aspect/gravity_model/radial_with_tidal_potential.h b/include/aspect/gravity_model/radial_with_tidal_potential.h index cd82ef43f79..ea60bf09fb6 100644 --- a/include/aspect/gravity_model/radial_with_tidal_potential.h +++ b/include/aspect/gravity_model/radial_with_tidal_potential.h @@ -32,17 +32,16 @@ namespace aspect { /** * A class that describes gravity as a radial vector of linearly - * changing magnitude by tidal potential from flattening and non-synchronnous rotation with respect to position and time. + * 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 it is from geodesy research, where potential is taken as positive. + * 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 planet, a_s = semimajor axis of satellite's orbit, b = angular rate of nonsynchronous rotation. + * 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. - * b [1/s] = b_NSR [m/yr] / (circumference of satellite [m]) / year_to_seconds [s/yr] * 2 * pi * * @ingroup GravityModels */ @@ -69,26 +68,20 @@ namespace aspect parse_parameters (ParameterHandler &prm) override; private: - /** - * Declare a member variable of type RadialConstant that allows us to call - * functions from radial_constant.cc. - */ - RadialConstant radialconstant; - /** * Mass of the perturbing body */ double M_p; /** - * Semimajor axis of the modeled body + * Semimajor axis of the orbit */ double a_s; /** - * Angular rate of the non-synchronous rotation in m/year + * Angular rate of the non-synchronous rotation in degrees/year */ - double b_NSR; + double b; }; } } diff --git a/source/gravity_model/radial_with_tidal_potential.cc b/source/gravity_model/radial_with_tidal_potential.cc index e4e44cf9765..899acc8aa0a 100644 --- a/source/gravity_model/radial_with_tidal_potential.cc +++ b/source/gravity_model/radial_with_tidal_potential.cc @@ -53,8 +53,8 @@ namespace aspect const double y = p[1]; const double z = p[2]; - const double C1 = std::cos( 2. * b_NSR * t); - const double C2 = std::sin( 2. * b_NSR * t); + const double C1 = std::cos( 2. * b * t); + const double C2 = std::sin( 2. * b * t); const double dTstar_over_dx = 1. / 6. * ( 2. * x ); const double dT0_over_dx = 1. / 2. * ( 2. * C1 * x - 2. * C2 * y ); @@ -66,18 +66,15 @@ namespace aspect const double dT0_over_dz = 0; const double G = aspect::constants::big_g; - const double T_factor = 3. * G * M_p / 2. / a_s / a_s / a_s; + const double T_factor = 3. * G * M_p / ( 2. * a_s * a_s * a_s ); const Tensor<1,dim> tidal_gravity = T_factor * ( (dTstar_over_dx + dT0_over_dx) * e_x + (dTstar_over_dy + dT0_over_dy) * e_y + (dTstar_over_dz + dT0_over_dz) * e_z); - if (p.norm() == 0.0) - return Tensor<1,dim>(); - - const double r = p.norm(); - return -radialconstant.magnitude * p/r + tidal_gravity; + RadialConstant radialconstant; + return radialconstant.gravity_vector(p) + tidal_gravity; } @@ -93,19 +90,20 @@ namespace aspect prm.declare_entry ("Mass of perturbing body", "1.898e27", Patterns::Double (), "Mass of body that perturbs gravity of modeled body. " - "Default value is for modeling Europa, therefore, mass of Jupiter. " + "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 (), - "Length of semimajor axis of orbit between modeled body and perturbing body. " - "Default value is for Europa's semimajor axis" + "The length of the semimajor axis of the orbit between the modeled body and the perturbing body. " + "The default value is for Europa's semimajor axis. " "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. " - "Default value is for Europa's icy shell. " - "Units is $degrees/year$"); + "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 (); } @@ -127,7 +125,7 @@ namespace aspect 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_NSR = prm.get_double ("Angular rate of nonsynchronous rotation") * constants::degree_to_radians / time_scale; + b = prm.get_double ("Angular rate of nonsynchronous rotation") * constants::degree_to_radians / time_scale; } prm.leave_subsection (); } @@ -161,7 +159,7 @@ namespace aspect } else { - Assert (false, ExcMessage ("This initial condition can only be used if the geometry " + Assert (false, ExcMessage ("This gravity model can only be used if the geometry " "is a sphere, a spherical shell, a chunk or an " "ellipsoidal chunk.")); } From a6daae7fb08b1301388cf465ba2df4125ffa9734 Mon Sep 17 00:00:00 2001 From: Hyunseong Kim Date: Tue, 7 Oct 2025 15:24:18 +0100 Subject: [PATCH 14/21] modifying template and tensor expressions --- .../radial_with_tidal_potential.h | 2 +- .../radial_with_tidal_potential.cc | 71 ++++--------------- tests/gravity_tidal_potential.prm | 4 +- tests/gravity_tidal_potential/screen-output | 62 +++++++++------- tests/gravity_tidal_potential/statistics | 22 +++--- 5 files changed, 65 insertions(+), 96 deletions(-) diff --git a/include/aspect/gravity_model/radial_with_tidal_potential.h b/include/aspect/gravity_model/radial_with_tidal_potential.h index ea60bf09fb6..fb7e1ab8ef9 100644 --- a/include/aspect/gravity_model/radial_with_tidal_potential.h +++ b/include/aspect/gravity_model/radial_with_tidal_potential.h @@ -74,7 +74,7 @@ namespace aspect double M_p; /** - * Semimajor axis of the orbit + * Semimajor axis of the orbit in consideration for tidal perturbation */ double a_s; diff --git a/source/gravity_model/radial_with_tidal_potential.cc b/source/gravity_model/radial_with_tidal_potential.cc index 899acc8aa0a..27730d4d918 100644 --- a/source/gravity_model/radial_with_tidal_potential.cc +++ b/source/gravity_model/radial_with_tidal_potential.cc @@ -26,10 +26,6 @@ #include #include -#include -#include -#include -#include namespace aspect { @@ -39,39 +35,35 @@ namespace aspect 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) */ - const Tensor<1,dim> e_x = Point::unit_vector(0); - const Tensor<1,dim> e_y = Point::unit_vector(1); - const Tensor<1,dim> e_z = Point::unit_vector(2); - const double t = (this->simulator_is_past_initialization()) ? this->get_time() : 0.0; - const double x = p[0]; - const double y = p[1]; - const double z = p[2]; const double C1 = std::cos( 2. * b * t); const double C2 = std::sin( 2. * b * t); - const double dTstar_over_dx = 1. / 6. * ( 2. * x ); - const double dT0_over_dx = 1. / 2. * ( 2. * C1 * x - 2. * C2 * y ); + const Tensor<1,dim> dTstar_gradient ({1./3. * p[0], 1./3. * p[1], -2./3. * p[2]}); - const double dTstar_over_dy = 1. / 6. * ( 2. * y ); - const double dT0_over_dy = 1. / 2. * ( -2. * C1 * y - 2. * C2 * x ); - - const double dTstar_over_dz = 1. / 6. * ( -4. * z ); - const double dT0_over_dz = 0; + 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_over_dx + dT0_over_dx) * e_x - + (dTstar_over_dy + dT0_over_dy) * e_y - + (dTstar_over_dz + dT0_over_dz) * e_z); + (dTstar_gradient + dT0_gradient); RadialConstant radialconstant; return radialconstant.gravity_vector(p) + tidal_gravity; @@ -95,13 +87,13 @@ namespace aspect prm.declare_entry ("Semimajor axis of orbit", "670900000", Patterns::Double (), "The length of the semimajor axis of the orbit between the modeled body and the perturbing body. " - "The default value is for Europa's semimajor axis. " + "The default value is for 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. " - "Default value is the angular rate of Europa's icy shell. " + "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. "); } @@ -130,39 +122,6 @@ namespace aspect prm.leave_subsection (); } prm.leave_subsection (); - - // This effect of tidal potential only works if the geometry is derived from - // a spherical model (i.e. a sphere, spherical shell or chunk) - if (Plugins::plugin_type_matches>(this->get_geometry_model())) - { - } - else if (Plugins::plugin_type_matches>(this->get_geometry_model())) - { - } - else if (Plugins::plugin_type_matches>(this->get_geometry_model())) - { - } - else if (Plugins::plugin_type_matches>(this->get_geometry_model())) - { - const auto &gm = Plugins::get_plugin_as_type> - (this->get_geometry_model()); - // TODO - // If the eccentricity of the EllipsoidalChunk is non-zero, the radius can vary along a boundary, - // but the maximal depth is the same everywhere and we could calculate a representative pressure - // profile. However, it requires some extra logic with ellipsoidal - // coordinates, so for now we only allow eccentricity zero. - // Using the EllipsoidalChunk with eccentricity zero can still be useful, - // because the domain can be non-coordinate parallel. - - AssertThrow(gm.get_eccentricity() == 0.0, - ExcNotImplemented("This plugin cannot be used with a non-zero eccentricity. ")); - } - else - { - Assert (false, ExcMessage ("This gravity model can only be used if the geometry " - "is a sphere, a spherical shell, a chunk or an " - "ellipsoidal chunk.")); - } } } } diff --git a/tests/gravity_tidal_potential.prm b/tests/gravity_tidal_potential.prm index 565d49d470f..0e3f47b936c 100644 --- a/tests/gravity_tidal_potential.prm +++ b/tests/gravity_tidal_potential.prm @@ -53,11 +53,11 @@ end subsection Gravity model - set Model name = radial linear with tidal potential + set Model name = radial with tidal potential subsection Radial constant set Magnitude = 1.3 end - subsection Radial linear with tidal potnetial + subsection Radial with tidal potential end end diff --git a/tests/gravity_tidal_potential/screen-output b/tests/gravity_tidal_potential/screen-output index 2c5db84ada0..1c22ad2f081 100644 --- a/tests/gravity_tidal_potential/screen-output +++ b/tests/gravity_tidal_potential/screen-output @@ -1,16 +1,21 @@ +----------------------------------------------------------------------------- +----------------------------------------------------------------------------- +----------------------------------------------------------------------------- +----------------------------------------------------------------------------- +----------------------------------------------------------------------------- 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)... 95+0 iterations. + Solving Stokes system (GMG)... 61+0 iterations. Postprocessing: - RMS, max velocity: 2.04e-05 m/year, 7.16e-05 m/year + 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.228e+05 Pa, 5.829e+07 Pa, 1.201e+08 Pa + Pressure min/avg/max: -5.226e+05 Pa, 0.002525 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 @@ -18,10 +23,10 @@ Number of degrees of freedom: 4,030 (2,910+150+970) Solving Stokes system (GMG)... 50+0 iterations. Postprocessing: - RMS, max velocity: 2.54e-05 m/year, 8.6e-05 m/year + 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.236e+05 Pa, 5.829e+07 Pa, 1.201e+08 Pa + Pressure min/avg/max: -5.226e+05 Pa, 0.002249 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 @@ -29,21 +34,21 @@ Number of degrees of freedom: 4,030 (2,910+150+970) Solving Stokes system (GMG)... 50+0 iterations. Postprocessing: - RMS, max velocity: 2.22e-05 m/year, 8.27e-05 m/year + 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.228e+05 Pa, 5.829e+07 Pa, 1.2e+08 Pa + Pressure min/avg/max: -5.226e+05 Pa, 0.002251 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)... 50+0 iterations. + Solving Stokes system (GMG)... 47+0 iterations. Postprocessing: - RMS, max velocity: 2.28e-05 m/year, 8.34e-05 m/year + 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.239e+05 Pa, 5.829e+07 Pa, 1.2e+08 Pa + Pressure min/avg/max: -5.226e+05 Pa, 0.002249 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 @@ -51,10 +56,10 @@ Number of degrees of freedom: 4,030 (2,910+150+970) Solving Stokes system (GMG)... 50+0 iterations. Postprocessing: - RMS, max velocity: 2.5e-05 m/year, 8.22e-05 m/year + 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.228e+05 Pa, 5.829e+07 Pa, 1.201e+08 Pa + Pressure min/avg/max: -5.226e+05 Pa, 0.002248 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 @@ -62,32 +67,32 @@ Number of degrees of freedom: 4,030 (2,910+150+970) Solving Stokes system (GMG)... 50+0 iterations. Postprocessing: - RMS, max velocity: 2.05e-05 m/year, 7.43e-05 m/year + 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.228e+05 Pa, 5.829e+07 Pa, 1.201e+08 Pa + Pressure min/avg/max: -5.226e+05 Pa, 0.002247 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)... 49+0 iterations. + Solving Stokes system (GMG)... 29+0 iterations. Postprocessing: - RMS, max velocity: 2.57e-05 m/year, 8.85e-05 m/year + 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.285e+05 Pa, 5.829e+07 Pa, 1.201e+08 Pa + Pressure min/avg/max: -5.226e+05 Pa, 0.002249 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. + Solving Stokes system (GMG)... 49+0 iterations. Postprocessing: - RMS, max velocity: 2.15e-05 m/year, 8.12e-05 m/year + 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.228e+05 Pa, 5.829e+07 Pa, 1.201e+08 Pa + 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 @@ -95,10 +100,10 @@ Number of degrees of freedom: 4,030 (2,910+150+970) Solving Stokes system (GMG)... 50+0 iterations. Postprocessing: - RMS, max velocity: 2.35e-05 m/year, 8.32e-05 m/year + 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.228e+05 Pa, 5.829e+07 Pa, 1.2e+08 Pa + Pressure min/avg/max: -5.226e+05 Pa, 0.002251 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 @@ -106,10 +111,10 @@ Number of degrees of freedom: 4,030 (2,910+150+970) Solving Stokes system (GMG)... 50+0 iterations. Postprocessing: - RMS, max velocity: 2.45e-05 m/year, 8.1e-05 m/year + 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.228e+05 Pa, 5.829e+07 Pa, 1.2e+08 Pa + Pressure min/avg/max: -5.226e+05 Pa, 0.002251 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 @@ -117,13 +122,18 @@ Number of degrees of freedom: 4,030 (2,910+150+970) Solving Stokes system (GMG)... 50+0 iterations. Postprocessing: - RMS, max velocity: 2.08e-05 m/year, 7.7e-05 m/year + 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.228e+05 Pa, 5.829e+07 Pa, 1.201e+08 Pa + Pressure min/avg/max: -5.226e+05 Pa, 0.002252 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 index 02ef831ab03..9533e172b3f 100644 --- a/tests/gravity_tidal_potential/statistics +++ b/tests/gravity_tidal_potential/statistics @@ -20,14 +20,14 @@ # 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 95 96 96 2.04033436e-05 7.15759491e-05 1.00000000e+02 9.99999535e+01 1.00000000e+02 output-gravity_tidal_potential/solution/solution-00000 -5.22786368e+05 5.82904414e+07 1.20123467e+08 9.17000000e+02 1.00000000e+20 2.63118556e+21 - 1 1.000000000000e+03 1.000000000000e+03 96 3060 970 0 50 51 51 2.53870602e-05 8.60009600e-05 1.00000000e+02 9.99999535e+01 1.00000000e+02 output-gravity_tidal_potential/solution/solution-00001 -5.23573277e+05 5.82904414e+07 1.20105088e+08 9.17000000e+02 1.00000000e+20 2.63118556e+21 - 2 2.000000000000e+03 1.000000000000e+03 96 3060 970 0 50 51 51 2.21503531e-05 8.26962723e-05 1.00000000e+02 9.99999535e+01 1.00000000e+02 output-gravity_tidal_potential/solution/solution-00002 -5.22786356e+05 5.82904414e+07 1.20011855e+08 9.17000000e+02 1.00000000e+20 2.63118556e+21 - 3 3.000000000000e+03 1.000000000000e+03 96 3060 970 0 50 51 51 2.28043259e-05 8.33561683e-05 1.00000000e+02 9.99999535e+01 1.00000000e+02 output-gravity_tidal_potential/solution/solution-00003 -5.23892832e+05 5.82904414e+07 1.19978816e+08 9.17000000e+02 1.00000000e+20 2.63118556e+21 - 4 4.000000000000e+03 1.000000000000e+03 96 3060 970 0 50 51 51 2.49936171e-05 8.22468188e-05 1.00000000e+02 9.99999535e+01 1.00000000e+02 output-gravity_tidal_potential/solution/solution-00004 -5.22786356e+05 5.82904414e+07 1.20075507e+08 9.17000000e+02 1.00000000e+20 2.63118556e+21 - 5 5.000000000000e+03 1.000000000000e+03 96 3060 970 0 50 51 51 2.04959462e-05 7.43144594e-05 1.00000000e+02 9.99999535e+01 1.00000000e+02 output-gravity_tidal_potential/solution/solution-00005 -5.22786357e+05 5.82904414e+07 1.20118224e+08 9.17000000e+02 1.00000000e+20 2.63118556e+21 - 6 6.000000000000e+03 1.000000000000e+03 96 3060 970 0 49 50 50 2.56511932e-05 8.84617285e-05 1.00000000e+02 9.99999535e+01 1.00000000e+02 output-gravity_tidal_potential/solution/solution-00006 -5.28450608e+05 5.82904414e+07 1.20124477e+08 9.17000000e+02 1.00000000e+20 2.63118556e+21 - 7 7.000000000000e+03 1.000000000000e+03 96 3060 970 0 50 51 51 2.15471701e-05 8.12170364e-05 1.00000000e+02 9.99999535e+01 1.00000000e+02 output-gravity_tidal_potential/solution/solution-00007 -5.22786357e+05 5.82904414e+07 1.20053747e+08 9.17000000e+02 1.00000000e+20 2.63118556e+21 - 8 8.000000000000e+03 1.000000000000e+03 96 3060 970 0 50 51 51 2.34666825e-05 8.31817721e-05 1.00000000e+02 9.99999535e+01 1.00000000e+02 output-gravity_tidal_potential/solution/solution-00008 -5.22786357e+05 5.82904414e+07 1.19974416e+08 9.17000000e+02 1.00000000e+20 2.63118556e+21 - 9 9.000000000000e+03 1.000000000000e+03 96 3060 970 0 50 51 51 2.44893181e-05 8.10215036e-05 1.00000000e+02 9.99999535e+01 1.00000000e+02 output-gravity_tidal_potential/solution/solution-00009 -5.22786359e+05 5.82904414e+07 1.20036187e+08 9.17000000e+02 1.00000000e+20 2.63118556e+21 -10 1.000000000000e+04 1.000000000000e+03 96 3060 970 0 50 51 51 2.07657184e-05 7.69854194e-05 1.00000000e+02 9.99999535e+01 1.00000000e+02 output-gravity_tidal_potential/solution/solution-00010 -5.22786358e+05 5.82904414e+07 1.20102575e+08 9.17000000e+02 1.00000000e+20 2.63118556e+21 + 0 0.000000000000e+00 0.000000000000e+00 96 3060 970 0 61 62 62 2.04002876e-05 7.18187431e-05 1.00000000e+02 9.99999535e+01 1.00000000e+02 output-gravity_tidal_potential/solution/solution-00000 -5.22625769e+05 2.52537853e-03 1.04524984e+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.53157546e-05 8.53402095e-05 1.00000000e+02 9.99999535e+01 1.00000000e+02 output-gravity_tidal_potential/solution/solution-00001 -5.22625789e+05 2.24906057e-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.30686075e-05 1.00000000e+02 9.99999535e+01 1.00000000e+02 output-gravity_tidal_potential/solution/solution-00002 -5.22625788e+05 2.25137919e-03 8.95530998e+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.30686169e-05 1.00000000e+02 9.99999535e+01 1.00000000e+02 output-gravity_tidal_potential/solution/solution-00003 -5.22625788e+05 2.24898011e-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.53157550e-05 8.53402007e-05 1.00000000e+02 9.99999535e+01 1.00000000e+02 output-gravity_tidal_potential/solution/solution-00004 -5.22625789e+05 2.24831020e-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.04002860e-05 7.18187334e-05 1.00000000e+02 9.99999535e+01 1.00000000e+02 output-gravity_tidal_potential/solution/solution-00005 -5.22625790e+05 2.24697792e-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.53402040e-05 1.00000000e+02 9.99999535e+01 1.00000000e+02 output-gravity_tidal_potential/solution/solution-00006 -5.22625787e+05 2.24930624e-03 1.01845607e+06 9.17000000e+02 1.00000000e+20 2.63118556e+21 + 7 7.000000000000e+03 1.000000000000e+03 96 3060 970 0 49 50 50 2.24054782e-05 8.30686193e-05 1.00000000e+02 9.99999535e+01 1.00000000e+02 output-gravity_tidal_potential/solution/solution-00007 -5.22625794e+05 2.24988606e-03 8.95530994e+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.24054785e-05 8.30686283e-05 1.00000000e+02 9.99999535e+01 1.00000000e+02 output-gravity_tidal_potential/solution/solution-00008 -5.22625791e+05 2.25079936e-03 8.95531002e+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.53157549e-05 8.53402019e-05 1.00000000e+02 9.99999535e+01 1.00000000e+02 output-gravity_tidal_potential/solution/solution-00009 -5.22625790e+05 2.25142332e-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.04002857e-05 7.18187306e-05 1.00000000e+02 9.99999535e+01 1.00000000e+02 output-gravity_tidal_potential/solution/solution-00010 -5.22625794e+05 2.25245757e-03 1.04524981e+06 9.17000000e+02 1.00000000e+20 2.63118556e+21 From 59b5abb4a26457a1a7617d22d32fa09e01b47637 Mon Sep 17 00:00:00 2001 From: Hyunseong Kim Date: Tue, 7 Oct 2025 16:47:18 +0100 Subject: [PATCH 15/21] adding further description on semimajor axis --- include/aspect/gravity_model/radial_with_tidal_potential.h | 2 +- source/gravity_model/radial_with_tidal_potential.cc | 6 ++++-- 2 files changed, 5 insertions(+), 3 deletions(-) diff --git a/include/aspect/gravity_model/radial_with_tidal_potential.h b/include/aspect/gravity_model/radial_with_tidal_potential.h index fb7e1ab8ef9..9213b7b0267 100644 --- a/include/aspect/gravity_model/radial_with_tidal_potential.h +++ b/include/aspect/gravity_model/radial_with_tidal_potential.h @@ -74,7 +74,7 @@ namespace aspect double M_p; /** - * Semimajor axis of the orbit in consideration for tidal perturbation + * Semimajor axis of the orbit that causes the tidal perturbation */ double a_s; diff --git a/source/gravity_model/radial_with_tidal_potential.cc b/source/gravity_model/radial_with_tidal_potential.cc index 27730d4d918..10b123430b3 100644 --- a/source/gravity_model/radial_with_tidal_potential.cc +++ b/source/gravity_model/radial_with_tidal_potential.cc @@ -86,8 +86,10 @@ namespace aspect "Units is $kg$."); prm.declare_entry ("Semimajor axis of orbit", "670900000", Patterns::Double (), - "The length of the semimajor axis of the orbit between the modeled body and the perturbing body. " - "The default value is for semimajor axis of Europa's orbit. " + "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 (), From ed5d7b5d7615f7ddc2c5dbb3d5bc343ce93212c4 Mon Sep 17 00:00:00 2001 From: Hyunseong Kim Date: Tue, 7 Oct 2025 16:55:32 +0100 Subject: [PATCH 16/21] removing unused lines on radial --- source/gravity_model/radial.cc | 5 ----- 1 file changed, 5 deletions(-) diff --git a/source/gravity_model/radial.cc b/source/gravity_model/radial.cc index c57dc1d9328..86b59a38141 100644 --- a/source/gravity_model/radial.cc +++ b/source/gravity_model/radial.cc @@ -25,11 +25,6 @@ #include -#include -#include -#include -#include - namespace aspect { namespace GravityModel From a268a78a3e8bb71c9ecb3a8f4a4c9c2a2e6b3eb0 Mon Sep 17 00:00:00 2001 From: Hyunseong Kim Date: Thu, 9 Oct 2025 13:15:04 +0100 Subject: [PATCH 17/21] correction on unused parameter error --- source/gravity_model/radial_with_tidal_potential.cc | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/source/gravity_model/radial_with_tidal_potential.cc b/source/gravity_model/radial_with_tidal_potential.cc index 10b123430b3..81f14fdf315 100644 --- a/source/gravity_model/radial_with_tidal_potential.cc +++ b/source/gravity_model/radial_with_tidal_potential.cc @@ -33,11 +33,11 @@ namespace aspect { template Tensor<1,dim> - RadialWithTidalPotential::gravity_vector (const Point &p) const + RadialWithTidalPotential::gravity_vector (const Point &/*p*/) const { - // This plugin is not implemented for 2D models - AssertThrow(false, ExcNotImplemented()); - return Tensor<1,dim>(); + // This plugin is not implemented for 2D models + AssertThrow(false, ExcNotImplemented()); + return Tensor<1,dim>(); } template <> From bf76c59ad26c821176f8a7ec6ff4eaf6ad859266 Mon Sep 17 00:00:00 2001 From: Hyunseong Kim Date: Thu, 9 Oct 2025 13:30:55 +0100 Subject: [PATCH 18/21] indent apply --- source/gravity_model/radial_with_tidal_potential.cc | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/source/gravity_model/radial_with_tidal_potential.cc b/source/gravity_model/radial_with_tidal_potential.cc index 81f14fdf315..912da3e6d7b 100644 --- a/source/gravity_model/radial_with_tidal_potential.cc +++ b/source/gravity_model/radial_with_tidal_potential.cc @@ -57,7 +57,7 @@ namespace aspect 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 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 ); @@ -111,7 +111,7 @@ namespace aspect { 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"); From 053571da9a5625ba6f884c1f2e46d4a9827963d1 Mon Sep 17 00:00:00 2001 From: Hyunseong Kim Date: Thu, 9 Oct 2025 15:11:24 +0100 Subject: [PATCH 19/21] test diff correction --- tests/gravity_tidal_potential/screen-output | 34 ++++++++------------- tests/gravity_tidal_potential/statistics | 22 ++++++------- 2 files changed, 23 insertions(+), 33 deletions(-) diff --git a/tests/gravity_tidal_potential/screen-output b/tests/gravity_tidal_potential/screen-output index 1c22ad2f081..cc8652dedb8 100644 --- a/tests/gravity_tidal_potential/screen-output +++ b/tests/gravity_tidal_potential/screen-output @@ -1,21 +1,16 @@ ------------------------------------------------------------------------------ ------------------------------------------------------------------------------ ------------------------------------------------------------------------------ ------------------------------------------------------------------------------ ------------------------------------------------------------------------------ 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)... 61+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.002525 Pa, 1.045e+06 Pa + 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 @@ -26,7 +21,7 @@ Number of degrees of freedom: 4,030 (2,910+150+970) 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.002249 Pa, 1.018e+06 Pa + 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 @@ -37,7 +32,7 @@ Number of degrees of freedom: 4,030 (2,910+150+970) 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.002251 Pa, 8.955e+05 Pa + 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 @@ -48,7 +43,7 @@ Number of degrees of freedom: 4,030 (2,910+150+970) 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.002249 Pa, 8.955e+05 Pa + 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 @@ -59,7 +54,7 @@ Number of degrees of freedom: 4,030 (2,910+150+970) 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.002248 Pa, 1.018e+06 Pa + 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 @@ -70,7 +65,7 @@ Number of degrees of freedom: 4,030 (2,910+150+970) 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.002247 Pa, 1.045e+06 Pa + 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 @@ -81,12 +76,12 @@ Number of degrees of freedom: 4,030 (2,910+150+970) 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.002249 Pa, 1.018e+06 Pa + 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)... 49+0 iterations. + Solving Stokes system (GMG)... 50+0 iterations. Postprocessing: RMS, max velocity: 2.24e-05 m/year, 8.31e-05 m/year @@ -103,7 +98,7 @@ Number of degrees of freedom: 4,030 (2,910+150+970) 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.002251 Pa, 8.955e+05 Pa + 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 @@ -114,7 +109,7 @@ Number of degrees of freedom: 4,030 (2,910+150+970) 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.002251 Pa, 1.018e+06 Pa + 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 @@ -125,15 +120,10 @@ Number of degrees of freedom: 4,030 (2,910+150+970) 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.002252 Pa, 1.045e+06 Pa + 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 index 9533e172b3f..1340acc34a5 100644 --- a/tests/gravity_tidal_potential/statistics +++ b/tests/gravity_tidal_potential/statistics @@ -20,14 +20,14 @@ # 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 61 62 62 2.04002876e-05 7.18187431e-05 1.00000000e+02 9.99999535e+01 1.00000000e+02 output-gravity_tidal_potential/solution/solution-00000 -5.22625769e+05 2.52537853e-03 1.04524984e+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.53157546e-05 8.53402095e-05 1.00000000e+02 9.99999535e+01 1.00000000e+02 output-gravity_tidal_potential/solution/solution-00001 -5.22625789e+05 2.24906057e-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.30686075e-05 1.00000000e+02 9.99999535e+01 1.00000000e+02 output-gravity_tidal_potential/solution/solution-00002 -5.22625788e+05 2.25137919e-03 8.95530998e+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.30686169e-05 1.00000000e+02 9.99999535e+01 1.00000000e+02 output-gravity_tidal_potential/solution/solution-00003 -5.22625788e+05 2.24898011e-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.53157550e-05 8.53402007e-05 1.00000000e+02 9.99999535e+01 1.00000000e+02 output-gravity_tidal_potential/solution/solution-00004 -5.22625789e+05 2.24831020e-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.04002860e-05 7.18187334e-05 1.00000000e+02 9.99999535e+01 1.00000000e+02 output-gravity_tidal_potential/solution/solution-00005 -5.22625790e+05 2.24697792e-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.53402040e-05 1.00000000e+02 9.99999535e+01 1.00000000e+02 output-gravity_tidal_potential/solution/solution-00006 -5.22625787e+05 2.24930624e-03 1.01845607e+06 9.17000000e+02 1.00000000e+20 2.63118556e+21 - 7 7.000000000000e+03 1.000000000000e+03 96 3060 970 0 49 50 50 2.24054782e-05 8.30686193e-05 1.00000000e+02 9.99999535e+01 1.00000000e+02 output-gravity_tidal_potential/solution/solution-00007 -5.22625794e+05 2.24988606e-03 8.95530994e+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.24054785e-05 8.30686283e-05 1.00000000e+02 9.99999535e+01 1.00000000e+02 output-gravity_tidal_potential/solution/solution-00008 -5.22625791e+05 2.25079936e-03 8.95531002e+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.53157549e-05 8.53402019e-05 1.00000000e+02 9.99999535e+01 1.00000000e+02 output-gravity_tidal_potential/solution/solution-00009 -5.22625790e+05 2.25142332e-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.04002857e-05 7.18187306e-05 1.00000000e+02 9.99999535e+01 1.00000000e+02 output-gravity_tidal_potential/solution/solution-00010 -5.22625794e+05 2.25245757e-03 1.04524981e+06 9.17000000e+02 1.00000000e+20 2.63118556e+21 + 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 From dc77c7161d38b8da9908256d0ba9d18e0eeffbdd Mon Sep 17 00:00:00 2001 From: Hyunseong Kim Date: Fri, 17 Oct 2025 16:33:04 +0100 Subject: [PATCH 20/21] changing parameter from angular rate to period --- .../radial_with_tidal_potential.h | 5 +++-- .../radial_with_tidal_potential.cc | 20 ++++++++++--------- 2 files changed, 14 insertions(+), 11 deletions(-) diff --git a/include/aspect/gravity_model/radial_with_tidal_potential.h b/include/aspect/gravity_model/radial_with_tidal_potential.h index 9213b7b0267..f32dee19e78 100644 --- a/include/aspect/gravity_model/radial_with_tidal_potential.h +++ b/include/aspect/gravity_model/radial_with_tidal_potential.h @@ -41,6 +41,7 @@ namespace aspect * (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 @@ -79,9 +80,9 @@ namespace aspect double a_s; /** - * Angular rate of the non-synchronous rotation in degrees/year + * Period of the non-synchronous rotation in year or second */ - double b; + double P; }; } } diff --git a/source/gravity_model/radial_with_tidal_potential.cc b/source/gravity_model/radial_with_tidal_potential.cc index 912da3e6d7b..0cc840acf3a 100644 --- a/source/gravity_model/radial_with_tidal_potential.cc +++ b/source/gravity_model/radial_with_tidal_potential.cc @@ -49,11 +49,13 @@ namespace aspect * 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 C1 = std::cos( 2. * b * t); - const double C2 = std::sin( 2. * b * t); + 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]}); @@ -91,13 +93,13 @@ namespace aspect "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", + prm.declare_entry ("Period of nonsynchronous rotation", "10000", Patterns::Double (), - "Angular rate of nonsynchronous rotation (NSR). " + "Period 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. "); + "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 (); } @@ -119,7 +121,7 @@ namespace aspect 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; + P = prm.get_double ("Period of nonsynchronous rotation") * time_scale; } prm.leave_subsection (); } @@ -143,6 +145,6 @@ namespace aspect "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''.") + "from a section `Gravity model/Radial with tidal potential'.") } } From 4986161a7b05b4119264bf4263cfe0130605c154 Mon Sep 17 00:00:00 2001 From: Hyunseong Kim Date: Fri, 17 Oct 2025 21:31:58 +0100 Subject: [PATCH 21/21] indent correction --- include/aspect/gravity_model/radial_with_tidal_potential.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/include/aspect/gravity_model/radial_with_tidal_potential.h b/include/aspect/gravity_model/radial_with_tidal_potential.h index f32dee19e78..6107417f9e0 100644 --- a/include/aspect/gravity_model/radial_with_tidal_potential.h +++ b/include/aspect/gravity_model/radial_with_tidal_potential.h @@ -41,7 +41,7 @@ namespace aspect * (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. + * 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