centrifugal force

[migueldvb]

I think that the centrifugal force is not included when the coriolis directive is enabled. I wonder if it makes sense to define a rotating_frame directive that includes the coriolis and centrifugal fictitious forces in gravity.F90? It would be nice to define the rotation frequency of the frame in the parameter file instead of using always the keplerian frequency for two point masses.

[gawrysz]

I think that if you plan to use rotating grids extensively, we should drop the CORIOLIS preprocesor symbol and include the forces related to non-inertial coordinates as a normal module. When rotation frequency is 0 the contribution should be just skipped to minimize overhead. There should also be an option to either give rotation frequency directly or calculate it by routine that implements gravity from binary system. This means that we should ensure proper precedence between initialization of gravity and coriolis, put proper values to the log files and make the code die whenever the combinatin of parameters and/or solvers is invalid or unimplemented.

[migueldvb]

I will make a module under src/base that includes the fictitious forces in rotating coordinates, should it be called src/base/non_inertial.F90?

[gawrysz]

Good idea.

[migueldvb]

I am not sure how to get the current grid container coordinates to calculate the centrifugal force. The rotacc array has and index for each fluid and the cg%x and cg%y have different shape so the compilation fails in line 131 and 133: https://github.com/migueldvb/piernik/blob/non-inertial/src/base/non_inertial.F90 Is there any example of an acceleration that uses the positions of the grid centers?

[gawrysz]

I guess it is best to add explicit loops over first index (fluids).

[migueldvb]

Thanks. This has been fixed now in the non-inertial branch.

When I use a namelist with an underscore as in NON_INERTIAL I get the error message below. This is not very important because removing the underscore in the namelist solves the problem.

mpif90  -fdefault-real-8 -ffree-form -std=gnu -fimplicit-none -ffree-line-length-none -g -O2 -fstack-arrays -c non_inertial.F90
non_inertial.F90:58.14:

      namelist // omega
              1
Error: Syntax error in NAMELIST statement at (1)
non_inertial.F90:70.25:

         write(nh%lun,nml=)
                         1
Error: Syntax error in WRITE statement at (1)
non_inertial.F90:74.26:

         read(unit=nh%lun, nml=, iostat=nh%ierrh, iomsg=nh%errstr)
                          1
Error: Syntax error in READ statement at (1)
non_inertial.F90:77.29:

         read(nh%cmdl_nml,nml=, iostat=nh%ierrh)
                             1
Error: Syntax error in READ statement at (1)
non_inertial.F90:80.25:

         write(nh%lun,nml=)
                         1

[Xarthisius]

The underscore is not the cause (there are several other namelists using it). My bet is that you defined precompiler macro NON_INERTIAL and it's being substituted to empty string.

[migueldvb]

Yes, that was the problem. Is it OK to have NON_INERTIAL as precompiler macro and a NONINERTIAL namelist containing the omega parameter or should we have different names?

[gawrysz]

Best would be to avoid using precompiler macros for something that is included in the code.

[migueldvb]

Now the fictitious acceleration is included in line 486 in rtvd.F90 by checking if the NON_INERTIAL macro is defined. https://github.com/migueldvb/piernik/blob/non-inertial/src/scheme/rtvd_split/rtvd.F90 If this macro is not used should we check if the value of omega is greater than zero and then apply the non-inertial forces?

[Xarthisius]

I was meant to rewrite rtvd to use function pointers, but I kept postponing it "till after my PhD". I guess I have no excuses now... @migueldvb keep implementing your module the way you're currently doing, i.e. using precompiler macro. I'll convert everything after we merge your work.

[Xarthisius]

As for the NONINERTIAL vs NON_INERTIAL for namelist and macro it's completely up to you. We don't have it standardized.

[migueldvb]

OK thank you. I have tested the advection problem using non-inertial module in my branch and it is not working well at the moment. The force is added after the sweep as in the coriolis force implementation but it does not have any effect. Here is a movie of the advection test in rotating frame. https://youtu.be/msZtP5Op8_E I wil try to debug this.

[migueldvb]

I have tried to inlcude the non-inertial force after the sweep in different places but it does not have any effect. This is the result of the advection test in Cartesian coordinates including non-inertial force: https://youtu.be/5xX1ESDoT70 The modified test and rtvd.F90 that I am working on are also included in the branch https://github.com/migueldvb/piernik/tree/non-inertial

[gawrysz]

For the cylindrical test of non-inertial forces we have 2 important cases:

• stationary pulse, which should be just advected around the grid (remember to set periodic domain in y-direction),
• a pulse flying in the opposite directon to the grid rotation, which should draw funny loops (cardioids?). Try e.g the following gnuplot code: se si sq ; se param ; vy=1; x0=1; om=-.5; p [-5:5][0:4][-2:2] x0*cos(om*t)-vy*t*sin(om*t),x0*sin(om*t)+vy*t*cos(om*t). On this trajectory the rectangular pulse should do one full rotation.

While the first test should run quite smoothly with minimal radial diffusion, the second one might be quite difficult for the current RTVD solver.

[migueldvb]

Thanks. I have tried to run a stationary pulse case in cylindrical coordinates and did not get advection. I think that I will work first on the cartesian advection_test problem. The pulse moves in a straight line accross the grid but it is not deflected either when the non-inertial forces are included in the same way as in the Coriolis module.

[gawrysz]

I think that it is worth to think a while about initproblem.F90 in advection_test. IIRC it has been extended for cylindrical tests in a quick and dirty way. Now when we add new module we should perhaps clarify in what coordinate system we define the initial conditons and how they have to be interpreted. In general both the coordinate system and things like angular frequency should go into initproblem and take part in defining initial velocity and pulse coordinates.

[migueldvb]

I have changed the initial position and velocities of the pulse to do some tests. I think that the angular frequency of the reference frame is not needed in initproblem.F90? Now the omega value is read in initpiernik.F90 and this is used to derive the non-inertial force that is applied after the source terms in rtvd.F90.

[gawrysz]

Well, everything depends on what interpretation we give to the initial position and velocitiy components of the pulse. Currently the pulse_vel parameter contains cartesian velocity. In rotating cylindrical coordinates, the initial velocity field must depend on angular frequency, if we want uniform velocity field without shearing. Also the calculate_error_norm will need angular frequency to correctly calculate the position and shape of the density pulse.