Boussinesq dynamo: nondimensionalized using viscous timescale

'''

This code provides an example for using a custom non-dimensionalization
of Rayleigh in Boussinesq dynamo mode.

The non-dimensionalization used here is based on the Rayleigh default
viscous diffusion scaling for convection, thus providing a check
on the custom reference state.

The parameters correspond to a case listed in Table 2 of
Soderlund et al.: "The influence of magnetic fields in planetary
dynamo models", Earth Planet. Sci. Lett., v.333-334, p.9-20 (2012)

The numbers referenced below for the various functions and constants
refer to equation (5) in "Rayleigh_Output_Variables.pdf". Please refer to
that document for further details.

Requirements: (1) "rayleigh_diagnostics.py" ; and (2) "reference_tools.py"

'''

import numpy as np


import os, sys
sys.path.insert(0, os.path.abspath('../../'))

import post_processing.reference_tools as rt

# name of output file containing custom reference data
filename = 'custom_ref_viscous.dat'

# Non-dimensional input parameters

Ra = 1.12e5             # Rayleigh number
Pr = 1.0                # Prandtl number
Ek = 2.0e-3             # Ekman number
Pm = 5.0                # Magnetic Prandtl number
beta = 0.4              # Aspect ratio = r_inner/r_outer
gravity_power = 1.0     # power law variation of gravitational acceleration

# Create radial grid

# Numer of radial grid points for radial functions (f_1, f_2, etc.)
# Make large enough for accurate interpolation onto Chebyshev grid
nr = 2000

# non-dimensional r_inner
ri = beta/(1-beta)

# non-dimensional r_outer
ro=1.0/(1-beta)

# non-dimensional radial grid
radius=np.linspace(ri,ro,nr)

# Define the reference state functions and constants

ones = np.ones(nr,dtype='float64')
zeros = np.zeros(nr,dtype='float64')

# the function list below is default for Boussinesq
f_1 = ones
f_2 = (radius/radius[nr-1])**gravity_power
f_3 = ones
f_4 = ones
f_5 = ones
f_6 = zeros
f_7 = ones
f_8 = zeros
f_9 = zeros
f_10 = zeros
f_11 = zeros
f_12 = zeros
f_13 = zeros

c_1 = 2.0/Ek                # Coriolis force
c_2 = Ra/Pr                 # Buoyancy force
c_3 = 1.0/Ek                # Pressure gradient
c_4 = 1.0/(Ek*Pm)           # Lorentz force
c_5 = 1.0                   # Viscous force
c_6 = 1.0/Pr                # Thermal diffusion
c_7 = 1.0/Pm                # Ohmic diffusion
c_8 = 0.0                   # Viscous heating
c_9 = 0.0                   # Ohmic heating
c_10 = 0.0                  # Internal heating

# Set all of the functions and constants

my_ref = rt.equation_coefficients(radius)

# Set functions here
my_ref.set_function(f_1, 1)
my_ref.set_function(f_2, 2)
my_ref.set_function(f_3, 3)
my_ref.set_function(f_4, 4)
my_ref.set_function(f_5, 5)
my_ref.set_function(f_6, 6)
my_ref.set_function(f_7, 7)
my_ref.set_function(f_8, 8)
my_ref.set_function(f_9, 9)
my_ref.set_function(f_10, 10)
my_ref.set_function(f_11, 11)
my_ref.set_function(f_12, 12)
my_ref.set_function(f_13, 13)

# Set constants here
my_ref.set_constant(c_1, 1)
my_ref.set_constant(c_2, 2)
my_ref.set_constant(c_3, 3)
my_ref.set_constant(c_4, 4)
my_ref.set_constant(c_5, 5)
my_ref.set_constant(c_6, 6)
my_ref.set_constant(c_7, 7)
my_ref.set_constant(c_8, 8)
my_ref.set_constant(c_9, 9)
my_ref.set_constant(c_10, 10)

my_ref.write(filename)

print('Custom reference file', filename, 'was written successfully.')