Skip to content

Commit 11bb687

Browse files
committed
Merge branch 'master' of ssh://github.com/magic-sph/magic
2 parents 090ba5b + 3e26d19 commit 11bb687

12 files changed

Lines changed: 269 additions & 73 deletions

File tree

samples/magic_wizard.py

Lines changed: 8 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -31,6 +31,7 @@
3131
import testGraphMovieOutputs.unitTest
3232
import testdtBMovieOutputs.unitTest
3333
import testTOGeosOutputs.unitTest
34+
import tides_radial_velocity.unitTest
3435

3536
__version__ = '1.0'
3637

@@ -363,6 +364,11 @@ def getSuite(startdir, cmd, precision, args):
363364
'outputFileDiff',
364365
'{}/testdtBMovieOutputs'.format(startdir),
365366
execCmd=cmd))
367+
# Tides: radial velocity
368+
suite.addTest(tides_radial_velocity.unitTest.TidesTest('outputFileDiff',
369+
'{}/tides_radial_velocity'.format(startdir),
370+
execCmd=cmd,
371+
precision=precision))
366372

367373
return suite
368374

@@ -409,7 +415,8 @@ def printLevelInfo():
409415
print(" Test r.m.s. force balance calculation ")
410416
print(" Test Graphic and Movie outputs ")
411417
print(" Test TO and Geos outputs ")
412-
418+
print(" Test Tides: radial velocity ")
419+
print(" ")
413420

414421
if __name__ == '__main__':
415422

samples/tides_radial_velocity/__init__.py

Whitespace-only changes.
Lines changed: 16 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,16 @@
1+
localhost
2+
localhost
3+
localhost
4+
localhost
5+
localhost
6+
localhost
7+
localhost
8+
localhost
9+
localhost
10+
localhost
11+
localhost
12+
localhost
13+
localhost
14+
localhost
15+
localhost
16+
localhost
Lines changed: 42 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,42 @@
1+
&grid
2+
n_r_max =33,
3+
n_cheb_max =31,
4+
n_phi_tot =192,
5+
n_r_ic_max =17,
6+
n_cheb_ic_max=15,
7+
minc =1,
8+
/
9+
&control
10+
mode =1,
11+
tag ="test",
12+
n_time_steps=5,
13+
dtmax =1.0D-8,
14+
alpha =0.6D0,
15+
runHours =12,
16+
runMinutes =00,
17+
time_scheme ='BPR353',
18+
l_non_rot =.true.,
19+
/
20+
&phys_param
21+
ra =1.0D5,
22+
pr =1.0D0,
23+
radratio =0.35D0,
24+
ktops =1,
25+
kbots =1,
26+
ktopv =2,
27+
kbotv =2,
28+
/
29+
&start_field
30+
l_start_file=.false.,
31+
/
32+
&output_control
33+
n_log_step =5,
34+
n_graphs =1,
35+
n_rsts =0,
36+
n_stores =0,
37+
runid ="Tides test",
38+
/
39+
&mantle
40+
amp_tide = 1.0d0,
41+
omega_tide = 1.0d0,
42+
/
Lines changed: 111 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,111 @@
1+
import unittest
2+
import numpy as np
3+
import glob
4+
import os
5+
import time
6+
import shutil
7+
import subprocess as sp
8+
9+
def cleanDir(dir):
10+
if os.path.exists('{}/pscond.dat'.format(dir)):
11+
os.remove('{}/pscond.dat'.format(dir))
12+
if os.path.exists('{}/scond.dat'.format(dir)):
13+
os.remove('{}/scond.dat'.format(dir))
14+
if os.path.exists('{}/run_magic.sh'.format(dir)):
15+
os.remove('{}/run_magic.sh'.format(dir))
16+
if os.path.exists('{}/run_magic_mpi.sh'.format(dir)):
17+
os.remove('{}/run_magic_mpi.sh'.format(dir))
18+
for f in glob.glob('{}/*_BIS'.format(dir)):
19+
os.remove(f)
20+
for f in glob.glob('{}/*.test'.format(dir)):
21+
os.remove(f)
22+
if os.path.exists('{}/stdout.out'.format(dir)):
23+
os.remove('{}/stdout.out'.format(dir))
24+
for f in glob.glob('{}/*.pyc'.format(dir)):
25+
os.remove(f)
26+
if os.path.exists('{}/__pycache__'.format(dir)):
27+
shutil.rmtree('{}/__pycache__'.format(dir))
28+
29+
30+
def get_ratios():
31+
"""
32+
Get the ratio of the temperature to radial velocity at CMB.
33+
"""
34+
from magic import MagicGraph
35+
G = MagicGraph(tag='test')
36+
37+
expr = ( 0.35 * (6/7) * np.sin(G.colatitude)**2 /
38+
( (2/7) * (3 * np.cos(G.colatitude)**2 - 1) +
39+
(8/7) * np.sin(G.colatitude)**2) )
40+
41+
maxratio = np.zeros(G.ntheta)
42+
for i in range(G.ntheta):
43+
maxratio[i] = G.entropy[:,i,0].max()/G.vr[:,i,0].max()
44+
45+
return expr, maxratio
46+
47+
class TidesTest(unittest.TestCase):
48+
49+
def __init__(self, testName, dir, execCmd='mpirun -n 8 ../tmp/magic.exe',
50+
precision=1e-8):
51+
super(TidesTest, self).__init__(testName)
52+
self.dir = dir
53+
self.precision = precision
54+
self.execCmd = execCmd
55+
self.startDir = os.getcwd()
56+
self.description = "Tides: radial velocity"
57+
58+
def list2reason(self, exc_list):
59+
if exc_list and exc_list[-1][0] is self:
60+
return exc_list[-1][1]
61+
62+
def setUp(self):
63+
# Cleaning when entering
64+
print('\nDirectory : {}'.format(self.dir))
65+
print('Description : {}'.format(self.description))
66+
self.startTime = time.time()
67+
cleanDir(self.dir)
68+
os.chdir(self.dir)
69+
cmd = '{} {}/input.nml'.format(self.execCmd, self.dir)
70+
sp.call(cmd, shell=True, stdout=open(os.devnull, 'wb'),
71+
stderr=open(os.devnull, 'wb'))
72+
73+
def tearDown(self):
74+
# Cleaning when leaving
75+
os.chdir(self.startDir)
76+
cleanDir(self.dir)
77+
78+
t = time.time()-self.startTime
79+
st = time.strftime("%M:%S", time.gmtime(t))
80+
print('Time used : {}'.format(st))
81+
82+
if hasattr(self, '_outcome'): # python 3.4+
83+
if hasattr(self._outcome, 'errors'): # python 3.4-3.10
84+
result = self.defaultTestResult()
85+
self._feedErrorsToResult(result, self._outcome.errors)
86+
else: # python 3.11+
87+
result = self._outcome.result
88+
else: # python 2.7-3.3
89+
result = getattr(self, '_outcomeForDoCleanups',
90+
self._resultForDoCleanups)
91+
92+
error = self.list2reason(result.errors)
93+
failure = self.list2reason(result.failures)
94+
ok = not error and not failure
95+
96+
if ok:
97+
print('Validating results.. OK')
98+
else:
99+
if error:
100+
print('Validating results.. ERROR!')
101+
print('\n')
102+
print(result.errors[-1][-1])
103+
if failure:
104+
print('Validating results.. FAIL!')
105+
print('\n')
106+
print(result.failures[-1][-1])
107+
108+
def outputFileDiff(self):
109+
datRef, datTmp = get_ratios()
110+
np.testing.assert_allclose(datRef, datTmp, rtol=self.precision,
111+
atol=1e-7)

src/Namelists.f90

Lines changed: 10 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -742,9 +742,16 @@ subroutine readNamelists(tscheme)
742742
l_AM=l_AM .or. l_correct_AMe .or. l_correct_AMz
743743
l_AM=l_AM .or. l_power
744744

745-
! Radial flow BC?
746-
if (ellipticity_cmb /= 0.0_cp .or. ellipticity_icb /= 0.0_cp .or. &
747-
& amp_tide /=0.0_cp ) l_radial_flow_bc=.true.
745+
! Radial flow BCs?
746+
l_vr_cmb=.false.
747+
l_vr_icb=.false.
748+
if (ellipticity_cmb /= 0.0_cp .or. amp_tide /=0.0_cp ) l_vr_cmb=.true.
749+
if (ellipticity_icb /= 0.0_cp ) l_vr_icb=.true.
750+
l_vr_bc = l_vr_cmb .or. l_vr_icb
751+
752+
if ( l_vr_bc .and. (ktops /= 1 .or. kbots /= 1) ) then
753+
call abortRun('! Radial flow boundary conditions not yet implemented for ktops /=1 and kbots /= 1')
754+
end if
748755

749756
!-- Heat boundary condition
750757
if ( impS /= 0 ) then

src/nonlinear_bcs.f90

Lines changed: 5 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -10,6 +10,7 @@ module nonlinear_bcs
1010
use horizontal_data, only: cosTheta, sinTheta_E2, phi, sinTheta
1111
use constants, only: two, y10_norm, y11_norm, zero
1212
use useful, only: abortRun
13+
use special, only: l_vr_cmb, l_vr_icb
1314
use sht, only: spat_to_sphertor
1415

1516
implicit none
@@ -141,11 +142,14 @@ subroutine v_rigid_boundary(nR,omega,lDeriv,vrr,vtr,vpr,cvrr,dvrdtr, &
141142
!-- Local variables:
142143
real(cp) :: r2
143144
integer :: nPhi
145+
logical :: l_vr
144146

145147
if ( nR == n_r_cmb ) then
146148
r2=r_cmb*r_cmb
149+
l_vr=l_vr_cmb
147150
else if ( nR == n_r_icb ) then
148151
r2=r_icb*r_icb
152+
l_vr=l_vr_icb
149153
else
150154
write(output_unit,*)
151155
write(output_unit,*) '! v_rigid boundary called for grid'
@@ -155,7 +159,7 @@ subroutine v_rigid_boundary(nR,omega,lDeriv,vrr,vtr,vpr,cvrr,dvrdtr, &
155159

156160
!$omp parallel do default(shared)
157161
do nPhi=1,n_phi_max
158-
vrr(:,nPhi)=0.0_cp
162+
if ( .not. l_vr ) vrr(:,nPhi)=0.0_cp
159163
vtr(:,nPhi)=0.0_cp
160164
vpr(:,nPhi)=r2*rho0(nR)*sinTheta_E2(:)*omega
161165
if ( lDeriv ) then

src/preCalculations.f90

Lines changed: 8 additions & 10 deletions
Original file line numberDiff line numberDiff line change
@@ -44,11 +44,10 @@ module preCalculations
4444
use integration, only: rInt_R
4545
use useful, only: logWrite, abortRun
4646
use special, only: l_curr, fac_loop, loopRadRatio, amp_curr, Le, n_imp, &
47-
& l_imp, l_radial_flow_bc, &
48-
& ellipticity_cmb, ellip_fac_cmb, &
47+
& l_imp, ellipticity_cmb, ellip_fac_cmb, &
4948
& ellipticity_icb, ellip_fac_icb, &
5049
& tide_fac20, tide_fac22p, tide_fac22n, &
51-
& amp_tide
50+
& amp_tide, l_vr_bc
5251
use time_schemes, only: type_tscheme
5352

5453
implicit none
@@ -769,7 +768,7 @@ subroutine preCalc(tscheme)
769768
n_s_max = n_r_max+int(r_icb*n_r_max)
770769
n_s_max = int(sDens*n_s_max)
771770

772-
if (l_radial_flow_bc) then
771+
if (l_vr_bc) then
773772
if ( ellipticity_cmb /= 0.0_cp ) then
774773
ellip_fac_cmb=-two*r_cmb**3*ellipticity_cmb*omega_ma1*omegaOsz_ma1/6.0_cp
775774
else
@@ -784,12 +783,11 @@ subroutine preCalc(tscheme)
784783

785784
if ( amp_tide /= 0.0_cp ) then
786785
! 6 factor = l*(l+1) for l=2
787-
y20_norm = 0.5_cp * sqrt(5.0_cp/pi) * 6.0_cp
788-
y22_norm = 0.25_cp * sqrt(7.5_cp/pi) * 6.0_cp
789-
tide_fac20 = amp_tide / y20_norm * r_cmb**2 ! (2,0,1) mode of Ogilvie 2014
790-
tide_fac22p = half * amp_tide / y22_norm / sqrt(6.0_cp) * r_cmb**2 ! Needs a half factor, (2,2,1) mode
791-
tide_fac22n = -7.0_cp * tide_fac22p ! Half factor carried over, (2,2,3) mode,
792-
! has opposite sign to that of the other two (Polfliet & Smeyers, 1990)
786+
y20_norm = 0.5_cp * sqrt(5.0_cp/pi) * 6.0_cp / r_cmb**2
787+
y22_norm = 0.25_cp * sqrt(7.5_cp/pi) * 6.0_cp / r_cmb**2
788+
tide_fac22n = 0.5_cp * amp_tide / y22_norm ! Needs a half factor, (2,2,3) mode
789+
tide_fac22p = tide_fac22n / 7.0_cp ! (2,2,1) mode of Ogilvie 2014
790+
tide_fac20 = 4.0_cp/7.0_cp * amp_tide / y20_norm ! (2,0,1) mode of Ogilvie 2014
793791
else
794792
tide_fac20 = 0.0_cp
795793
tide_fac22p = 0.0_cp

src/rIter.f90

Lines changed: 20 additions & 11 deletions
Original file line numberDiff line numberDiff line change
@@ -54,7 +54,7 @@ module rIter_mod
5454
use RMS, only: get_nl_RMS, transform_to_lm_RMS, compute_lm_forces, &
5555
& transform_to_grid_RMS
5656
use probe_mod
57-
use special, only: ellip_fac_icb, l_radial_flow_bc
57+
use special, only: l_vr_cmb, l_vr_icb
5858

5959
implicit none
6060

@@ -568,17 +568,26 @@ subroutine transform_to_grid_space(this, nR, nBc, lViscBcCalc, lRmsCalc, &
568568
& this%gsa%dvtdpc, this%gsa%dvpdpc, l_R(nR))
569569
end if
570570
else if ( nBc == 2 ) then
571-
if ( nR == n_r_cmb .and. (.not. l_radial_flow_bc) ) then
571+
if ( nR == n_r_cmb ) then
572+
if ( l_vr_cmb ) then
573+
call torpol_to_spat(w_Rloc(:,nR), dw_Rloc(:,nR), z_Rloc(:,nR), &
574+
& this%gsa%vrc, this%gsa%vtc, this%gsa%vpc, l_R(nR))
575+
end if
572576
call v_rigid_boundary(nR, omega_ma, lDeriv, this%gsa%vrc, &
573-
& this%gsa%vtc, this%gsa%vpc, this%gsa%cvrc, &
574-
& this%gsa%dvrdtc, this%gsa%dvrdpc, &
575-
& this%gsa%dvtdpc,this%gsa%dvpdpc)
576-
else if ( nR == n_r_icb .and. ellip_fac_icb == 0.0_cp ) then
577-
call v_rigid_boundary(nR, omega_ic, lDeriv, this%gsa%vrc, &
578-
& this%gsa%vtc, this%gsa%vpc, &
579-
& this%gsa%cvrc, this%gsa%dvrdtc, &
580-
& this%gsa%dvrdpc, this%gsa%dvtdpc, &
581-
& this%gsa%dvpdpc)
577+
& this%gsa%vtc, this%gsa%vpc, this%gsa%cvrc, &
578+
& this%gsa%dvrdtc, this%gsa%dvrdpc, &
579+
& this%gsa%dvtdpc, this%gsa%dvpdpc)
580+
end if
581+
582+
if ( nR == n_r_icb ) then
583+
if ( l_vr_icb ) then
584+
call torpol_to_spat(w_Rloc(:,nR), dw_Rloc(:,nR), z_Rloc(:,nR), &
585+
& this%gsa%vrc, this%gsa%vtc, this%gsa%vpc, l_R(nR))
586+
end if
587+
call v_rigid_boundary(nR, omega_ic, lDeriv, this%gsa%vrc, &
588+
& this%gsa%vtc, this%gsa%vpc, this%gsa%cvrc, &
589+
& this%gsa%dvrdtc, this%gsa%dvrdpc, &
590+
& this%gsa%dvtdpc, this%gsa%dvpdpc)
582591
end if
583592
if ( lDeriv ) then
584593
call torpol_to_spat(dw_Rloc(:,nR), ddw_Rloc(:,nR), dz_Rloc(:,nR), &

src/special.f90

Lines changed: 3 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -47,8 +47,9 @@ module special
4747
real(cp), public :: tide_fac20 ! Pre-factor for tidal forcing for (2,0)
4848
real(cp), public :: tide_fac22p ! Pre-factor for tidal forcing for (2,2)
4949
real(cp), public :: tide_fac22n ! Pre-factor for tidal forcing for (2,-2)
50-
logical, public :: l_radial_flow_bc ! A flag that tells whether radial flow BCs are being used
51-
50+
logical, public :: l_vr_cmb ! Radial flow BCs at CMB
51+
logical, public :: l_vr_icb ! Radial flow BCs at ICB
52+
logical, public :: l_vr_bc ! Radial flow BCs
5253
real(cp), public :: ampForce ! Amplitude of external body force
5354

5455
end module special

0 commit comments

Comments
 (0)