Merge branch 'master' of ssh://github.com/magic-sph/magic

python/magic/bLayers.py

@@ -5,7 +5,10 @@
 from magic import MagicRadial, matder, intcheb, MagicSetup, scanDir, AvgField
 from magic.setup import labTex
 from scipy.signal import argrelextrema
-from scipy.integrate import simps
+try:
+    from scipy.integrate import simps
+except:
+    from scipy.integrate import simpson as simps
 from scipy.interpolate import splrep, splev
 
 def getAccuratePeaks(rad, uh, uhTop, uhBot, ri, ro):

python/magic/butterfly.py

@@ -8,8 +8,11 @@
 import scipy.interpolate as S
 from magic.plotlib import cut
 from magic.setup import labTex
-from scipy.integrate import simps
 from .npfile import *
+try:
+    from scipy.integrate import simps
+except:
+    from scipy.integrate import simpson as simps
 
 
 class Butterfly:

python/magic/compSims.py

@@ -1,11 +1,13 @@
 # -*- coding: utf-8 -*-
 from magic import *
-from scipy.integrate import trapz
 from matplotlib.ticker import ScalarFormatter
 import matplotlib.pyplot as plt
 import numpy as np
 import os
-
+try:
+    from scipy.integrate import trapz
+except:
+    from scipy.integrate import trapezoid as trapz
 
 class CompSims:
     """

python/magic/libmagic.py

@@ -1,10 +1,13 @@
 # -*- coding: utf-8 -*-
 import scipy.interpolate as S
-from scipy.integrate import simps
 import numpy as np
 import glob, os, re, sys
 import subprocess as sp
 from .npfile import *
+try:
+    from scipy.integrate import simps
+except:
+    from scipy.integrate import simpson as simps
 
 
 def selectField(obj, field, labTex=True, ic=False):
@@ -561,7 +564,7 @@ def fd_grid(nr, a, b, fd_stretching=0.3, fd_ratio=0.1):
         dr_after = dr_before
         for i in range(1, nr):
             rr[i]=rr[i-1]-dr_before
-    
+
     else:
         n_boundary_points = int( float(nr-1)/(2.*(1+ratio1)) )
         ratio1 = float(nr-1)/float(2.*n_boundary_points) -1.
@@ -574,7 +577,7 @@ def fd_grid(nr, a, b, fd_stretching=0.3, fd_ratio=0.1):
             dr_before = dr_before*dr_after
         dr_before = 1./(float(n_bulk_points)+ \
                     2.*dr_after*((1-dr_before)/(1.-dr_after)))
-    
+
         for i in range(n_boundary_points):
             dr_before = dr_before*dr_after
 

python/magic/surf.py

@@ -7,7 +7,10 @@
 import matplotlib.pyplot as plt
 import os
 import numpy as np
-from scipy.integrate import trapz
+try:
+    from scipy.integrate import trapz
+except:
+    from scipy.integrate import trapezoid as trapz
 
 
 class Surf:
@@ -869,7 +872,7 @@ def avg(self, field='vphi', levels=defaultLevels, cm=None,
 
             vs = self.gr.vr * np.sin(th3D) + self.gr.vtheta * np.cos(th3D)
             vz = self.gr.vr * np.cos(th3D) - self.gr.vtheta * np.sin(th3D)
-            vortz = 1./s3D*(-phideravg(vs, self.gr.minc) + 
+            vortz = 1./s3D*(-phideravg(vs, self.gr.minc) +
                             sderavg(s3D*self.gr.vphi, self.gr.radius))
 
             data = vortz * vz

python/magic/thHeat.py

@@ -3,7 +3,11 @@
 import numpy as np
 from magic import scanDir, MagicSetup, Movie, chebgrid, rderavg, AvgField
 import os, pickle
-from scipy.integrate import simps, trapz
+try:
+    from scipy.integrate import simps, trapz
+except:
+    from scipy.integrate import simpson as simps
+    from scipy.integrate import trapezoid as trapz
 
 json_model = {'phys_params': [],
               'time_series': { 'heat': ['topnuss', 'botnuss']},

src/Makefile

@@ -16,7 +16,7 @@ USE_PRECOND=yes
 # USE_FFTLIB can be MKL, FFTW or JW
 USE_FFTLIB=MKL
 #USE_DCTLIB can be MKL, FFTW or JW
-USE_DCTLIB=JW
+USE_DCTLIB=MKL
 #USE_LAPACKLIB can be MKL, LAPACK or JW
 USE_LAPACKLIB=MKL
 #Use shtns for Legendre/Fourier transforms
@@ -295,16 +295,15 @@ endif
 
 OBJS   := $(addsuffix .o, $(basename $(RED_SOURCES)))
 OBJS += truncation.o
-OBJS += c_utils.o
 
 .SUFFIXES: 
 
-ifeq ($(COMPILER),gnu)
-ifeq ($(USE_MPI),yes)
-all: mpi.mod
-	make $(OUT)
-endif
-endif
+#ifeq ($(COMPILER),gnu)
+#ifeq ($(USE_MPI),yes)
+#all: mpi.mod
+#	make $(OUT)
+#endif
+#endif
 
 $(OUT): $(OBJS)
 	$(FC) $(FFLAGS) -o $@ $^ $(LIBS)
@@ -363,7 +362,7 @@ get_nl.o              : truncation.o horizontal.o logic.o time_schemes.o\
                         phys_param.o radial.o blocking.o mem_alloc.o
 out_movie_file.o      : truncation.o output_data.o blocking.o\
                         radial.o horizontal.o movie.o fields.o\
-                        num_param.o $(FFT_OBJS) logic.o\
+                        num_param.o $(FFT_OBJS) logic.o outMisc.o\
                         out_dtB_frame.o parallel.o communications.o
 storeCheckPoints.o    : truncation.o phys_param.o num_param.o logic.o\
                         output_data.o init_fields.o dt_fieldsLast.o\
@@ -381,7 +380,8 @@ power.o               : truncation.o blocking.o horizontal.o\
                         phys_param.o num_param.o radial.o logic.o\
                         output_data.o outRot.o integration.o useful.o\
                         mem_alloc.o mean_sd.o
-integration.o         : $(DCT_OBJS) precision_mod.o constants.o
+integration.o         : $(DCT_OBJS) precision_mod.o constants.o chebyshev.o\
+                        communications.o
 out_coeff.o           : logic.o precision_mod.o parallel.o blocking.o\
                         truncation.o communications.o\
                         mem_alloc.o radial.o output_data.o phys_param.o
@@ -395,9 +395,9 @@ magic.o               : truncation.o num_param.o parallel.o logic.o\
 step_time.o           : truncation.o phys_param.o updateB.o updateZ.o\
                         num_param.o radial_data.o updateWP.o time_schemes.o\
                         logic.o output_data.o output.o movie.o updateS.o\
-                        timing.o courant.o signals.o updateXI.o\
+                        timing.o courant.o signals.o updateXI.o updatePHI.o\
                         radialLoop.o nonlinear_bcs.o updateWPS.o\
-                        LMLoop.o dt_fieldsLast.o c_utils.o useful.o probes.o
+                        LMLoop.o dt_fieldsLast.o useful.o probes.o
 phys_param.o          : precision_mod.o
 char_manip.o          : precision_mod.o
 constants.o           : precision_mod.o
@@ -408,7 +408,7 @@ chebyshev_polynoms.o  : constants.o logic.o precision_mod.o num_param.o
 rIteration.o          : precision_mod.o time_schemes.o truncation.o radial_data.o
 rIter.o               : precision_mod.o parallel.o truncation.o logic.o radial_data.o\
                         radial.o constants.o get_nl.o get_td.o TO.o dtB.o\
-                        out_graph_file.o nonlinear_bcs.o\
+                        out_graph_file.o nonlinear_bcs.o out_movie_file.o\
                         courant.o outRot.o fields.o time_schemes.o phys_param.o\
                         ${SHT_OBJS} outGeos.o
 radialLoop.o          : truncation.o time_schemes.o precision_mod.o mem_alloc.o\
@@ -420,9 +420,9 @@ LMLoop.o              : truncation.o blocking.o parallel.o time_array.o\
                         updateZ.o updateWP.o  debugging.o
 mpi_transpose.o       : precision_mod.o truncation.o radial_data.o\
                         mem_alloc.o parallel.o blocking.o
-debugging.o           : precision_mod.o
+debugging.o           : precision_mod.o precision_mod.o
 timing.o              : parallel.o precision_mod.o
-radial_derivatives.o  : $(DCT_OBJS) constants.o precision_mod.o
+radial_derivatives.o  : $(DCT_OBJS) constants.o precision_mod.o chebyshev.o
 output.o              : truncation.o blocking.o phys_param.o\
                         num_param.o logic.o output_data.o radial.o\
                         horizontal.o constants.o fields.o radial_spectra.o\
@@ -456,17 +456,20 @@ updateWPS.o           : truncation.o blocking.o num_param.o time_array.o\
                         phys_param.o radial.o horizontal.o logic.o\
                         $(DCT_OBJS) radial_derivatives.o\
                         mem_alloc.o init_fields.o RMS.o time_schemes.o
+updatePHI.o           : truncation.o blocking.o num_param.o time_array.o\
+                        phys_param.o radial.o horizontal.o logic.o\
+                        $(DCT_OBJS) radial_derivatives.o init_fields.o\
+                        mem_alloc.o time_schemes.o
 get_td.o              : truncation.o blocking.o horizontal.o\
                         phys_param.o num_param.o radial.o logic.o\
-                        RMS.o RMS_helpers.o fields.o\
-                        cutils_iface.o mem_alloc.o
+                        RMS.o RMS_helpers.o fields.o mem_alloc.o
 store_movie_IC.o      : truncation.o blocking.o logic.o radial_data.o\
                         movie.o radial.o horizontal.o\
                         $(SHT_OBJS) $(FFT_OBJS) out_movie_file.o phys_param.o
 radial_spectra.o      : truncation.o blocking.o horizontal.o radial.o\
                         num_param.o output_data.o logic.o useful.o\
                         char_manip.o radial_data.o LMmapping.o constants.o
-fft_fac.o             : constants.o precision_mod.o
+fft_fac.o             : constants.o precision_mod.o mem_alloc.o
 cosine_transform_even.o: truncation.o fft_fac.o constants.o useful.o\
                         mem_alloc.o
 readCheckPoints.o     : truncation.o blocking.o phys_param.o time_array.o\
@@ -505,7 +508,8 @@ updateB.o             : truncation.o blocking.o horizontal.o fields.o\
                         RMS_helpers.o mem_alloc.o
 outGeos.o             : truncation.o horizontal.o radial.o constants.o\
                         output_data.o logic.o radial_data.o parallel.o\
-                        communications.o integration.o mem_alloc.o logic.o
+                        communications.o integration.o mem_alloc.o logic.o\
+                        movie.o
 preCalculations.o     : truncation.o phys_param.o num_param.o constants.o\
                         radial.o horizontal.o init_fields.o time_schemes.o\
                         blocking.o logic.o output_data.o\
@@ -530,7 +534,7 @@ radial.o              : truncation.o radial_data.o $(LAPACK_OBJS)\
                         chebyshev_polynoms.o radial_derivatives.o\
 			cosine_transform_even.o mem_alloc.o useful.o radial_scheme.o\
                         chebyshev.o finite_differences.o
-output_data.o         : precision_mod.o
+output_data.o         : precision_mod.o mem_alloc.o
 horizontal.o          : truncation.o phys_param.o num_param.o precision_mod.o\
                         radial.o logic.o blocking.o plms.o $(FFT_OBJS)\
                         mem_alloc.o
@@ -539,7 +543,7 @@ init_fields.o         : truncation.o blocking.o radial.o horizontal.o\
                         constants.o logic.o $(FFT_OBJS) $(SHT_OBJS)\
                         useful.o phys_param.o mpi_transpose.o\
                         $(DCT_OBJS) mem_alloc.o parallel.o 
-mean_sd.o             : precision_mod.o mem_alloc.o matrices.o
+mean_sd.o             : precision_mod.o mem_alloc.o matrices.o communications.o
 movie.o               : truncation.o parallel.o radial_data.o\
                         output_data.o logic.o radial.o mem_alloc.o\
                         horizontal.o char_manip.o useful.o
@@ -562,7 +566,7 @@ multistep_schemes.o   : parallel.o precision_mod.o num_param.o constants.o\
                         time_array.o
 dt_fieldsLast.o       : truncation.o blocking.o precision_mod.o time_array.o\
                         logic.o mem_alloc.o constants.o
-fields.o              : truncation.o blocking.o radial_data.o\
+fields.o              : truncation.o blocking.o radial_data.o communications.o\
                         precision_mod.o	mem_alloc.o logic.o
 special.o             : precision_mod.o mem_alloc.o truncation.o
 probes.o              : parallel.o precision_mod.o truncation.o radial.o num_param.o\

src/Namelists.f90

@@ -134,7 +134,8 @@ subroutine readNamelists(tscheme)
       &    omega_ma1,omegaOsz_ma1,tShift_ma1,             &
       &    omega_ma2,omegaOsz_ma2,tShift_ma2,             &
       &    amp_mode_ma,omega_mode_ma,m_mode_ma,           &
-      &    mode_symm_ma,ellipticity_cmb
+      &    mode_symm_ma,ellipticity_cmb, amp_tide,        &
+      &    omega_tide
 
       namelist/inner_core/sigma_ratio,nRotIc,rho_ratio_ic, &
       &    omega_ic1,omegaOsz_ic1,tShift_ic1,              &
@@ -741,6 +742,10 @@ subroutine readNamelists(tscheme)
       l_AM=l_AM .or. l_correct_AMe .or. l_correct_AMz
       l_AM=l_AM .or. l_power
 
+      ! Radial flow BC?
+      if (ellipticity_cmb /= 0.0_cp .or. ellipticity_icb /= 0.0_cp .or. &
+      &   amp_tide /=0.0_cp ) l_radial_flow_bc=.true.
+
       !-- Heat boundary condition
       if ( impS /= 0 ) then
          rad=pi/180
@@ -1223,6 +1228,8 @@ subroutine writeNamelists(n_out)
       write(n_out,'(''  m_mode_ma       ='',i4,'','')')  m_mode_ma
       write(n_out,'(''  mode_symm_ma    ='',i4,'','')')  mode_symm_ma
       write(n_out,'(''  ellipticity_cmb ='',ES14.6,'','')') ellipticity_cmb
+      write(n_out,'(''  amp_tide        ='',ES14.6,'','')') amp_tide
+      write(n_out,'(''  omega_tide      ='',ES14.6,'','')') omega_tide
       write(n_out,*) "/"
 
       write(n_out,*) "&inner_core"
@@ -1354,6 +1361,8 @@ subroutine defaultNamelists
       radratio    =0.35_cp
       dilution_fac=0.0_cp    ! centrifugal acceleration
       ampForce    =0.0_cp    ! External body force amplitude
+      amp_tide    =0.0_cp    ! Amplitude of tide
+      omega_tide  =0.0_cp    ! Frequency of tide
 
       !-- Phase field
       tmelt       =0.0_cp    ! Melting temperature

src/fields.f90

@@ -7,7 +7,7 @@ module fields
    use precision_mod
    use mem_alloc, only: bytes_allocated
    use constants, only: zero
-   use physical_parameters, only: ampForce
+   use special, only: ampForce
    use truncation, only: lm_max, n_r_max, lm_maxMag, n_r_maxMag, &
        &                 n_r_ic_maxMag, fd_order, fd_order_bound
    use logic, only: l_chemical_conv, l_finite_diff, l_mag, l_parallel_solve, &

src/init_fields.f90

@@ -40,7 +40,8 @@ module init_fields
        &                          impXi, n_impXi_max, n_impXi, phiXi,     &
        &                          thetaXi, peakXi, widthXi, osc, epscxi,  &
        &                          kbotxi, ktopxi, BuoFac, ktopp, oek,     &
-       &                          epsPhase, ampForce
+       &                          epsPhase
+   use special, only: ampForce
    use algebra, only: prepare_mat, solve_mat
    use cosine_transform_odd
    use dense_matrices

src/phys_param.f90

@@ -108,12 +108,6 @@ module physical_parameters
    real(cp) :: penaltyFac     ! Factor that enters the penalty method in the NS equations
    real(cp) :: tmelt          ! Melting temperature
 
-   real(cp) :: ellipticity_cmb    ! Ellipticity of CMB, used for libration
-   real(cp) :: ellipticity_icb    ! Ellipticity of ICB, used for libration
-   real(cp) :: ellip_fac_cmb      ! d/d\phi (Y22) * d/dt exp(i\omega t) at CMB
-   real(cp) :: ellip_fac_icb      ! d/d\phi (Y22) * d/dt exp(i\omega t) at ICB
-
-   real(cp) :: ampForce           ! Amplitude of external body force
    real(cp) :: gammatau_gravi     ! Constant in front of the core/mantle gravitationnal torque (see Aubert et al. 2013)
 
 end module physical_parameters

src/preCalculations.f90

@@ -39,13 +39,16 @@ module preCalculations
        &                          ktops, kbots, interior_model, r_LCR,     &
        &                          n_r_LCR, mode, tmagcon, oek, Bn, imagcon,&
        &                          ktopxi, kbotxi, epscxi, epscxi0, sc, osc,&
-       &                          ChemFac, raxi, Po, prec_angle,           &
-       &                          ellipticity_cmb, ellip_fac_cmb,          &
-       &                          ellipticity_icb, ellip_fac_icb
+       &                          ChemFac, raxi, Po, prec_angle
    use horizontal_data, only: horizontal
    use integration, only: rInt_R
    use useful, only: logWrite, abortRun
-   use special, only: l_curr, fac_loop, loopRadRatio, amp_curr, Le, n_imp, l_imp
+   use special, only: l_curr, fac_loop, loopRadRatio, amp_curr, Le, n_imp, &
+       &              l_imp, l_radial_flow_bc,                             &
+       &              ellipticity_cmb, ellip_fac_cmb,                      &
+       &              ellipticity_icb, ellip_fac_icb,                      &
+       &              tide_fac20, tide_fac22p, tide_fac22n,                &
+       &              amp_tide
    use time_schemes, only: type_tscheme
 
    implicit none
@@ -75,6 +78,7 @@ subroutine preCalc(tscheme)
       character(len=76) :: fileName
       character(len=80) :: message
       real(cp) :: mom(n_r_max)
+      real(cp) :: y20_norm, y22_norm
 
       !-- Determine scales depending on n_tScale,n_lScale :
       if ( n_tScale == 0 ) then
@@ -763,18 +767,32 @@ subroutine preCalc(tscheme)
          call abortRun('LCR not compatible with imposed field!')
       end if
 
-      if ( ellipticity_cmb /= 0.0_cp ) then
-         ellip_fac_cmb=-two*r_cmb**3*ellipticity_cmb*omega_ma1*omegaOsz_ma1
-      else
-         ellip_fac_cmb=0.0_cp
-      end if
+      if (l_radial_flow_bc) then
+         if ( ellipticity_cmb /= 0.0_cp ) then
+            ellip_fac_cmb=-two*r_cmb**3*ellipticity_cmb*omega_ma1*omegaOsz_ma1
+         else
+            ellip_fac_cmb=0.0_cp
+         end if
 
-      if ( ellipticity_icb /= 0.0_cp ) then
-         ellip_fac_icb=-two*r_icb**3*ellipticity_icb*omega_ic1*omegaOsz_ic1
-      else
-         ellip_fac_icb=0.0_cp
-      end if
+         if ( ellipticity_icb /= 0.0_cp ) then
+            ellip_fac_icb=-two*r_icb**3*ellipticity_icb*omega_ic1*omegaOsz_ic1
+         else
+            ellip_fac_icb=0.0_cp
+         end if
 
+         if ( amp_tide /= 0.0_cp ) then
+            y20_norm = 0.5_cp  * sqrt(5.0_cp/pi)
+            y22_norm = 0.25_cp * sqrt(7.5_cp/pi)
+            tide_fac20  = amp_tide / y20_norm * r_cmb**2 ! (2,0,1) mode of Ogilvie 2014
+            tide_fac22p = half * amp_tide / y22_norm / sqrt(6.0_cp) * r_cmb**2 ! Needs a half factor, (2,2,1) mode
+            tide_fac22n = -7.0_cp * tide_fac22p                                ! Half factor carried over, (2,2,3) mode,
+                                                                               ! has opposite sign to that of the other two (Polfliet & Smeyers, 1990)
+         else
+            tide_fac20  = 0.0_cp
+            tide_fac22p = 0.0_cp
+            tide_fac22n = 0.0_cp
+         end if
+      end if
    end subroutine preCalc
 !-------------------------------------------------------------------------------
    subroutine preCalcTimes(time,tEND,dt,n_time_step,n_time_steps)