Adding tidal potential to gravity model #6614

hyunseong96 · 2025-07-16T10:48:40Z

Apologies for opening this PR so late. Antoniette and I took some time to check sign convention for the equation.
I open this PR to add gravity model plugin where I can add tidal potential exerted as satellite orbits around the host planet.

The equation of tidal potential is from Tobie et al. (2025) (https://doi.org/10.1007/s11214-025-01136-y).

I only use T* and T0 because ep 4, 5 and 6 are dependent to short time scale (~3.55 days) and they are replaced by tidal heating.
T* is flattening and T0 is nonsynchronous rotation. T0 happens because ice surface is decoupled from rocky interior by the subsurface ocean and have difference rotational velocity.

To implement in cartesian coordinates, these equations are converted to ones using x, y and z instead of spherical variables.
which ends up as:
dV_T/dx = (3* G M_p)/(2 a^3) (1/6(2x)+1/2(2C1x - 2C2y))
dV_T/dy = (3* G M_p)/(2 a^3) (1/6(2y)+1/2(2C1y - 2C2x))
dV_T/dz = (3* G M_p)/(2 a^3) (1/6(-4z)).
where G is gravitational constant, M_p is mass of perturbing body, a is semi-major axis, C1 = cos(2bt) and C2 = sin(2b*t).

This $V_T$ is from geodesy research, which takes potential and acceleration always positive as unit vector is inward to center of body. Therefore, for us, where unit vector is outward from the center of body, the return of plugin is $-g\vec{r}-\nabla(-V_T)\vec{r}$

test model is just to see how gravity changes with time. temperature does not change with time. viscosity is high throughout the model so there is no convection. Colorbar is for gravity magnitude, and arrow denotes direction of tidal acceleration.

I will edit some parts quickly in test .prm, make more clarification in .cc and add change file!

hyunseong96 · 2025-07-16T13:36:02Z

Corrections are done but I am not sure why I have failure from dealii 9.5 test.
Is it okay if I just use git apply diff file?

bangerth · 2025-07-28T03:18:52Z

@hyunseong96 Thank you for this patch! I will look at it later this week, but in the meantime, can you take a look at the test failure? It currently reports this:

******* Check /__w/aspect/aspect/build/tests/output-gravity_tidal_potential/statistics /__w/aspect/aspect/tests/gravity_tidal_potential/statistics
******* First 50 of 93 lines of diffs are:
----------------
##23      <==  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 
##23      ==>  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 

##23      #:8   <== 110
##23      #:8   ==> 95
@ Absolute error = 1.5000000000e+1, Relative error = 1.5789473684e-1
##23      #:9   <== 112
##23      #:9   ==> 96
@ Absolute error = 1.6000000000e+1, Relative error = 1.6666666667e-1
##23      #:10  <== 112
##23      #:10  ==> 96
@ Absolute error = 1.6000000000e+1, Relative error = 1.6666666667e-1
##23      #:12  <== 7.15757443e-05
##23      #:12  ==> 7.15759315e-05
@ Absolute error = 1.8720000000e-10, Relative error = 2.6154111540e-6
##23      #:17  <== -5.22786493e+05
##23      #:17  ==> -5.22786377e+05
@ Absolute error = 1.1600000000e-1, Relative error = 2.2188795482e-7
----------------
##24      <==  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 
##24      ==>  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 

##24      #:12  <== 8.60011863e-05
##24      #:12  ==> 8.60009790e-05
@ Absolute error = 2.0730000000e-10, Relative error = 2.4104376765e-6
##24      #:17  <== -5.23573285e+05
##24      #:17  ==> -5.23573278e+05
@ Absolute error = 7.0000000000e-3, Relative error = 1.3369666280e-8
----------------
##25      <==  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 
##25      ==>  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 

##25      #:17  <== -5.22786489e+05
##25      #:17  ==> -5.22786367e+05
@ Absolute error = 1.2200000000e-1, Relative error = 2.3336492246e-7
----------------
##27      <==  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 
##27      ==>  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 

##27      #:12  <== 8.22470390e-05
##27      #:12  ==> 8.22468386e-05
@ Absolute error = 2.0040000000e-10, Relative error = 2.4365678172e-6
##27      #:17  <== -5.22786489e+05
##27      #:17  ==> -5.22786367e+05
@ Absolute error = 1.2200000000e-1, Relative error = 2.3336492246e-7
----------------
##28      <==  5 5.000000000000e+03 1.000000000000e+03 96 3060 970 0  50  51  51 2.0495[945](https://github.com/geodynamics/aspect/actions/runs/16318503026/job/46090091115?pr=6614#step:6:946)9e-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 
##28      ==>  5 5.000000000000e+03 1.000000000000e+03 96 3060 970 0 50 51 51 2.0495[946](https://github.com/geodynamics/aspect/actions/runs/16318503026/job/46090091115?pr=6614#step:6:947)0e-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 

ninja: build stopped: subcommand failed.

CMake Error at /__w/aspect/aspect/tests/run_test.cmake:14 (message):
  *** Test aborted.

hyunseong96 · 2025-07-30T13:32:35Z

@bangerth Thanks in advance! :) The issue with test is solved now.

hyunseong96 · 2025-09-16T15:46:49Z

Hi @bangerth , hope you are doing well! just a gentle reminder to ask for reviewing.
Meanwhile, I was advised to add you and @Geoniette to contributors on change log, so I modified change log in my latest push.
If there are points I should modify in the plugin before your review, I don't mind working on that!
Thanks for your help :)

gassmoeller

I will go ahead and leave some comments, Wolfgang can of course chime in on the mathematics (I havent compared the implemented formula with the paper). Some comments on the code and structure below. I think this is a valuable addition, but I would change the structure of the code a bit. Let me know in case you have questions or want some guidance on how to do something.

gassmoeller · 2025-09-22T12:48:07Z

include/aspect/gravity_model/radial.h

+
+
+
+    template <int dim>


As I explain in #6681 we usually associate one class/plugin with one file. The radial gravity models were an exception from the earliest days of ASPECT's life. Can you introduce your gravity profile in a new file instead? For reasons I will explain below I would name the new file (and class and everything) RadialWithTidalPotential and radial_with_tidal_potential.h/cc respectively.

I will make it into a different plugin file then!

gassmoeller · 2025-09-22T13:01:33Z

include/aspect/gravity_model/radial.h

+
+
+    template <int dim>
+    class RadialLinearWithTidalPotential : public RadialLinear<dim>


While it works to derive a plugin from another plugin if you want to use the behavior of the parent class in the child class, I have noted over the years that it is rarely a good idea in ASPECT. It creates complicated class structures, and sometimes it straight out prevents making full use of the parent class (in particular if the parent has protected member variables or is itself derived from another class like here). I would recommend the following approach instead:

Suggested change

class RadialLinearWithTidalPotential : public RadialLinear<dim>

class RadialLinearWithTidalPotential : public Interface<dim>, public SimulatorAccess<dim>

Then give this class a private member variable of type RadialLinear<dim> as a "base_gravity_model" and use its functionality by calling functions of this member variable. You can see how this works in the class initial_temperature/patch_on_S40RTS.h and its member variable S40RTSPerturbation<dim> s40rts;. Note in particular the lines 166-168 in initial_temperature/patch_on_S40RTS.cc where your class would initialize the member variable and let it read its own input parameters and line 76 where it uses its function get_Vs.

For now this is just a way to achieve the same functionality as your current approach (just with a bit more code, but easier to understand concept). But in the future this allows to expand the plugin to work with more than just the radial linear gravity model. For example the material model prescribed_viscosity can use an arbitrary material model as base model, and then modify it. In the same way your plugin could use any of the radial gravity profiles and apply a tidal potential on top.

gassmoeller · 2025-09-22T13:04:14Z

include/aspect/gravity_model/radial.h

+
+
+
+    template <int dim>


Please add doxygen documentation above the class name about what this class does. You can use the text you wrote at the bottom of the source file, but bring it into the doxygen format (/** ... */). This will be used in the online documentation (here the example for RadialLinear: https://aspect.geodynamics.org/doc/doxygen/classaspect_1_1GravityModel_1_1RadialLinear.html.

gassmoeller · 2025-09-22T13:05:38Z

include/aspect/gravity_model/radial.h

+        double b_NSR;
+
+        /**
+         * Radius of modelling satellite


By satellite you mean the moon/planet this gravity is computed for? This would be worth spelling out, since it may not be obvious for geoscientists focused on Earth.

Honestly, I haven't thought about the planet's perspective. My idea was about effect of gravitation by planet on satellite, because I wanted to implement gravity variation on Europa. That's why I used host planet, which affects the gravity of satellite.
Would it be better to differentiate as modeling object and gravitationally interacting object?
I also think I should add few more comments on parameters for cases other than icy moons.

gassmoeller · 2025-09-22T13:06:33Z

include/aspect/gravity_model/radial.h

+        double a_s;
+
+        /**
+         * Rate of non-synchronous rotation in m/year


Please use SI units wherever possible. You can still allow m/yr as input parameter if the global parameter Use years instead of seconds is set, and then convert into m/s when reading and storing the parameter.

And why m as length unit? Wouldnt it be much easier to ask the user to compute the angular velocity themselves beforehand and use that as input? Then you dont need to input the radius of the satellite either, or do you need that for anything else?

gassmoeller · 2025-09-22T13:19:29Z

source/gravity_model/radial.cc

+        {
+          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");


here you can do the following to allow for different time units in the input file, but still use SI units in the code:

Suggested change

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 ("Rate of nonsynchronous rotation") / time_scale;

gassmoeller · 2025-09-22T13:32:52Z

source/gravity_model/radial.cc

+        {
+          prm.declare_entry ("Mass of planet", "1.898e27",
+                             Patterns::Double (),
+                             "Mass of satellte's host planet. "


Suggested change

"Mass of satellte's host planet. "

"Mass of satellite's host planet. "

gassmoeller · 2025-09-22T13:34:46Z

source/gravity_model/radial.cc

+                             "Units is $kg$.");
+          prm.declare_entry ("Semimajor axis of satellite", "670900000",
+                             Patterns::Double (),
+                             "Length of semimajor axis of satellite. "


Better to be more verbose:

Suggested change

"Length of semimajor axis of satellite. "

"Length of the semimajor axis of the orbit of the modeled planet. "

gassmoeller · 2025-09-22T13:36:27Z

source/gravity_model/radial.cc

+    void
+    RadialLinearWithTidalPotential<dim>::parse_parameters (ParameterHandler &prm)
+    {
+      AssertThrow (dim==3, ExcMessage ("The 'radial linear with tidal potential' gravity model "


update name

gassmoeller · 2025-09-22T13:36:52Z

source/gravity_model/radial.cc

+        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");


The current spelling is easy to confuse, if you model an elliptical body on an elliptical orbit.

Suggested change

a_s = prm.get_double ("Semimajor axis of satellite");

a_s = prm.get_double ("Semimajor axis of orbit");

hyunseong96 · 2025-09-25T17:09:01Z

@gassmoeller Hi Rene, thanks for reviewing the PR. I have few questions.
First, when you mentioned using angular rate instead of m/year. I thought about just using angular rate as 1/s, but seeing the suggestion below confused me a bit.

const double time_scale = this->get_parameters().convert_to_years ? constants::year_in_seconds : 1.0;
b_NSR = prm.get_double ("Rate of nonsynchronous rotation") / time_scale;

Did you mean the parameter, b_NSR, as angle/year? or is there any perspective I am not aware of inside the code while solving equations?

Second, separating into different files seems a bit risky with my poor skill of git commands sadly. I found out that I need to add this line to radial constant plugin.
template <int dim2> friend class RadialWithTidalPotential;

After the incident of having hundred of others' commits on my PR at Hackathon, I am still not sure how I can update the branch so I can adjust my plugin to the update you will do about radial plugins. I searched for this and tried on branch I am not using anymore.

Sync fork of my repository (main) at github - git remote update - git checkout - git merge origin/main

Would this four steps be valid to keep up with updates? Otherwise I think it can be easier for me to close the PR and open a new one after your radial separation PR is merged.

gassmoeller · 2025-09-26T11:26:29Z

First, when you mentioned using angular rate instead of m/year. I thought about just using angular rate as 1/s, but seeing the suggestion below confused me a bit.
Did you mean the parameter, b_NSR, as angle/year? or is there any perspective I am not aware of inside the code while solving equations?

Yes I could have been more clear in my comment there. In general there are two separate issues:

I think expressing the rotation velocity in angular velocity instead of velocity makes more sense, because then you do not need the input of the satellite radius and the internal conversion to angular velocity.
In many places in ASPECT we expect that the input of velocities is given as unit per second (if the parameter Use years instead of seconds is false), or as unit per years (if the parameter is true). But internally we always want to use consistent SI units, so if the velocity is given by the input file as 1/yr then after reading it we need to convert to 1/s. This is all I wanted to express with my suggestions. The code I mentioned should do exactly that: Read the input, then if the input is given in years instead of seconds scale it with a time scale and convert to the proper SI unit 1/s:

const double time_scale = this->get_parameters().convert_to_years ? constants::year_in_seconds : 1.0;
b_NSR = prm.get_double ("Rate of nonsynchronous rotation") / time_scale;

I found out that I need to add this line to radial constant plugin.
template <int dim2> friend class RadialWithTidalPotential;

I do not see why you would need to do this. You can declare a class (in this RadialConstant) as a friend of another class (in this case RadialWithTidalPotential), which allows the other class to access private member variables of the first class. But I do not see where you do this, the only operation your new class needs of RadialConstant is to call its public member functions declare_parameters, parse_parameters and gravity_vector. No access to private variables is needed. Can you upload your full new code following the instructions I give below?

After the incident of having hundred of others' commits on my PR at Hackathon, I am still not sure how I can update the branch so I can adjust my plugin to the update you will do about radial plugins. I searched for this and tried on branch I am not using anymore.

Sync fork of my repository (main) at github - git remote update - git checkout - git merge origin/main

Would this four steps be valid to keep up with updates? Otherwise I think it can be easier for me to close the PR and open a new one after your radial separation PR is merged.

Please do not use this workflow, it is what caused all the commits to appear in your branch. It is a valid git workflow and some projects use it, but we do not. I could explain our alternate workflow and the reasons for it (using an operation called git rebase), but it is not necessary for what you want to do here. No update from the main branch is necessary, just continue on your branch and move all the new code you added into a new file and remove it from radial.h/.cc. Afterwards, there should not be any changes in radial.h/.cc left and therefore it will not conflict with my changes on the main branch. We will have to squash all your commits before merging this PR in the end, but I will write you instructions for it when we are at that stage. For now just keep making changes on this branch without updates from main.

hyunseong96 · 2025-10-01T17:57:12Z

@gassmoeller Thanks for clarifying! I changed parameters accordingly. I will modify the tests later once codes are sorted out!
In the this commit, I separated plugin to a new files. However, this still brings RadialConstant through private variable. Because plugin was not working with public variables.
I was trying to have plugin work with parse_parameter, declare_parameter and gravity_vector from RadialConstant but I kept getting errors as:

In file included from /Users/hyunseong/aspect-kim/Hackathons/tidal_potential/CMakeFiles/aspect.exe.release.dir/Unity/unity_8_cxx.cxx:19:
/Users/hyunseong/aspect/source/gravity_model/radial_with_tidal_potential.cc:76:73: error: call to non-static member function without an object argument
const Tensor<1,dim> Constant = GravityModel::RadialConstant::gravity_vector(p);
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~^~~~~~~~~~~~~~
/Users/hyunseong/aspect/source/gravity_model/radial_with_tidal_potential.cc:186:35: note: in instantiation of member function 'aspect::GravityModel::RadialWithTidalPotential<2>::gravity_vector' requested here
ASPECT_REGISTER_GRAVITY_MODEL(RadialWithTidalPotential,
^
/Users/hyunseong/aspect/source/gravity_model/radial_with_tidal_potential.cc:134:42: error: call to non-static member function without an object argument
GravityModel::RadialConstant::parse_parameters (prm);
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~^~~~~~~~~~~~~~~~
/Users/hyunseong/aspect/source/gravity_model/radial_with_tidal_potential.cc:186:35: note: in instantiation of member function 'aspect::GravityModel::RadialWithTidalPotential<2>::parse_parameters' requested here
ASPECT_REGISTER_GRAVITY_MODEL(RadialWithTidalPotential,
^
/Users/hyunseong/aspect/source/gravity_model/radial_with_tidal_potential.cc:76:73: error: call to non-static member function without an object argument
const Tensor<1,dim> Constant = GravityModel::RadialConstant::gravity_vector(p);
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~^~~~~~~~~~~~~~
/Users/hyunseong/aspect/source/gravity_model/radial_with_tidal_potential.cc:186:35: note: in instantiation of member function 'aspect::GravityModel::RadialWithTidalPotential<3>::gravity_vector' requested here
ASPECT_REGISTER_GRAVITY_MODEL(RadialWithTidalPotential,
^
/Users/hyunseong/aspect/source/gravity_model/radial_with_tidal_potential.cc:134:42: error: call to non-static member function without an object argument
GravityModel::RadialConstant::parse_parameters (prm);
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~^~~~~~~~~~~~~~~~
/Users/hyunseong/aspect/source/gravity_model/radial_with_tidal_potential.cc:186:35: note: in instantiation of member function 'aspect::GravityModel::RadialWithTidalPotential<3>::parse_parameters' requested here
ASPECT_REGISTER_GRAVITY_MODEL(RadialWithTidalPotential,

I can push the ones with public variables with errors if you want it!

gassmoeller

Yes, much better structure. I left a number of comments, but they are mostly suggestions for improvements, I didnt see anything wrong.

Based on the error messages that you posted and your questions, I think it is worth for you to review the difference between static functions and non static functions of a class in c++. Here is a short summary: https://www.learncpp.com/cpp-tutorial/static-member-functions/.

In short, 'normal' member functions of a class (say class RadialConstant) require an object of this class to be created (in your case radialconstant). Then you can call the member functions like so: radialconstant.gravity_vector(...). Static member functions however do not require an object of that class, and are usually called through their class name (instead of a particular object name), like so: RadialConstant<dim>::declare_parameters(...).

gassmoeller · 2025-10-02T08:13:04Z

source/gravity_model/radial_with_tidal_potential.cc

+        return Tensor<1,dim>();
+
+      const double r = p.norm();
+      return -radialconstant.magnitude * p/r + tidal_gravity;


Ah, here is where you access the private member variable magnitude of the object radialconstant which is of the class type RadialConstant<dim>. The way to work around that is to call the public member function gravity_vector like so:

return radialconstant.gravity_vector(p) + tidal_gravity;

Then you do no longer access the private member variable magnitude, and instead call the public member function gravity_vector.

gassmoeller · 2025-10-02T08:13:36Z

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 forces.


Suggested change

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 with tidal potential' for tidal forces.

gassmoeller · 2025-10-02T08:13:53Z

include/aspect/gravity_model/radial.h

+
+        template <int dim2> friend class RadialWithTidalPotential;


Remove this

gassmoeller · 2025-10-02T08:15:45Z

include/aspect/gravity_model/radial_with_tidal_potential.h

+  {
+    /**
+     * 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.


Suggested change

* 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.

gassmoeller · 2025-10-02T08:16:18Z

include/aspect/gravity_model/radial_with_tidal_potential.h

+     * 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.


Suggested change

* 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.

gassmoeller · 2025-10-02T09:00:25Z

source/gravity_model/radial_with_tidal_potential.cc

+    void
+    RadialWithTidalPotential<dim>::parse_parameters (ParameterHandler &prm)
+    {
+      AssertThrow (dim==3, ExcMessage ("The 'radial with tidal potential' gravity model "


Ah, so you only allow this for 3D models. Please add this to the documentation at the bottom of this file. I will leave a comment above for how to handle this in the gravity_vector function.

gassmoeller · 2025-10-02T09:07:27Z

source/gravity_model/radial_with_tidal_potential.cc

+    template <int dim>
+    Tensor<1,dim>
+    RadialWithTidalPotential<dim>::gravity_vector (const Point<dim> &p) const
+    {


Since you want this to only work in 3D models, but ASPECT is compiled for 2D and 3D models, you will (or maybe already have) run into trouble with using types like Tensor<1,dim>. When ASPECT is compiled for dim==2 the compiler will notice you try to access elements that dont exist in 2D and likely complain. Do the following to work around this (this is taken from here):

Suggested change

template <int dim>

Tensor<1,dim>

RadialWithTidalPotential<dim>::gravity_vector (const Point<dim> &p) const

{

template <int dim>

Tensor<1,dim>

RadialWithTidalPotential<dim>::gravity_vector (const Point<dim> &p) const

{

// This plugin is not implemented for 2D models

AssertThrow(false, ExcNotImplemented());

return Tensor<1,dim>();

}

template <>

Tensor<1,3>

RadialWithTidalPotential<dim>::gravity_vector (const Point<3> &p) const

{

const unsigned int dim = 3;

This construct is called a template specialization (more information here: https://www.learncpp.com/cpp-tutorial/function-template-specialization/), and it provides different functions for your class for dim==2 (which will crash) and dim==3 (which will work as you expect it to).

Thank you! but I had to change a bit to make this code work. I hope the code below is the correct modification.

template <>
Tensor<1,3>
RadialWithTidalPotential<3>::gravity_vector (const Point<3> &p) const
{

I will use this later for Tensor<1,2> when Antoniette and I decide how to implement 2D case.
So far, we had discussion how to implement this plugin in 2D but we will make it as different PR.

Hi @gassmoeller , I have a bit of trouble for upper part of suggestion, which is template

Tensor<1,dim> RadialWithTidalPotential<dim>::gravity_vector (const Point<dim> &p) const { // This plugin is not implemented for 2D models AssertThrow(false, ExcNotImplemented()); return Tensor<1,dim>(); }

While in make, I see the warning that warning:

unused parameter 'p' [-Wunused-parameter]

Because tests considers warnings as errors, I fail tests. Is there a simple but reasonable way I can use Point 'p'? or the way to avoid this warning?

Yes, our usual way to avoid this warning is to comment the name (but not the type) of the parameter argument. This tells the compiler that this parameter exists (for some reasons or other), but that we are intentionally not using it in the function:

RadialWithTidalPotential<dim>::gravity_vector (const Point<dim> &/*p*/) const

gassmoeller · 2025-10-02T09:11:27Z

source/gravity_model/radial_with_tidal_potential.cc

+        {
+          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;


Ah, so you do the conversion already. Then disregard my earlier comment, and just fix the documentation of this parameter to say that it is given in degrees/year or degrees/second depending on the input parameter 'Use years instead of seconds'.

gassmoeller · 2025-10-02T09:11:52Z

source/gravity_model/radial_with_tidal_potential.cc

+          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;


but according to your formulas, shouldnt this parameter now be b and not b_NSR?

gassmoeller · 2025-10-02T09:12:18Z

source/gravity_model/radial_with_tidal_potential.cc

+        }
+      else
+        {
+          Assert (false, ExcMessage ("This initial condition can only be used if the geometry "


Suggested change

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 "

hyunseong96 added 3 commits June 19, 2025 23:14

initiate gravity model for tidal effect

f6a4b3d

modifying files

1ecc8f3

type correction

c8e0290

test apply diff

8ba55ba

adding contributors to change log

98784e4

hyunseong96 added 2 commits September 16, 2025 23:34

changing time parameter name

d814c23

applying test diff

8a8e39f

gassmoeller mentioned this pull request Sep 22, 2025

Split radial gravity plugins into separate files #6681

Open

gassmoeller reviewed Sep 22, 2025

View reviewed changes

modifying by suggestions

66a6297

hyunseong96 added 2 commits September 25, 2025 18:14

minor corrections

545332f

minor corrections about plugin

58b8278

updating angular rate

19db0c7

gassmoeller reviewed Oct 2, 2025

View reviewed changes

hyunseong96 added 8 commits October 6, 2025 15:02

default value change and the status working with friend class

4478a28

working with public

2a0042b

modifying template and tensor expressions

a6daae7

adding further description on semimajor axis

59b5abb

removing unused lines on radial

ed5d7b5

correction on unused parameter error

a268a78

indent apply

bf76c59

test diff correction

053571d



		template <int dim>
		class RadialLinearWithTidalPotential : public RadialLinear<dim>

	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 ("Rate of nonsynchronous rotation") / time_scale;

	"Mass of satellte's host planet. "
	"Mass of satellite's host planet. "

	"Length of semimajor axis of satellite. "
	"Length of the semimajor axis of the orbit of the modeled planet. "

	a_s = prm.get_double ("Semimajor axis of satellite");
	a_s = prm.get_double ("Semimajor axis of orbit");

		@@ -0,0 +1,4 @@
		Added: ASPECT now has a new gravity model plugin called 'Radial linear with tidal potnetial' for tidal forces.

	* 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.

	* 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.

	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 "

Adding tidal potential to gravity model #6614

Are you sure you want to change the base?

Adding tidal potential to gravity model #6614

Uh oh!

Conversation

hyunseong96 commented Jul 16, 2025 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

hyunseong96 commented Jul 16, 2025

Uh oh!

bangerth commented Jul 28, 2025

Uh oh!

hyunseong96 commented Jul 30, 2025

Uh oh!

hyunseong96 commented Sep 16, 2025

Uh oh!

gassmoeller left a comment

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

hyunseong96 commented Sep 25, 2025

Uh oh!

gassmoeller commented Sep 26, 2025

Uh oh!

hyunseong96 commented Oct 1, 2025

Uh oh!

gassmoeller left a comment

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

hyunseong96 Oct 7, 2025 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

hyunseong96 commented Jul 16, 2025 •

edited

Loading

hyunseong96 Oct 7, 2025 •

edited

Loading