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

Author: tgastine

README.md

@@ -13,7 +13,7 @@
 * **MagIC** uses either Chebyshev polynomials or finite differences in the radial direction and spherical harmonic decomposition in the azimuthal and latitudinal directions.  MagIC supports several Implicit-Explicit time schemes where the nonlinear terms and the Coriolis force are treated explicitly, while the remaining linear terms are treated implicitly.
 
 
-* **MagIC** is written in Fortran and designed to be used on supercomputing clusters.  It thus relies on a hybrid parallelisation scheme using both [OpenMP](https://openmp.org/wp/) and [MPI](https://www.open-mpi.org/). Postprocessing functions written in python (requiring [matplotlib](https://matplotlib.org/) and [scipy](https://www.scipy.org/) are also provided to allow a useful data analysis.
+* **MagIC** is written in Fortran and designed to be used on supercomputing clusters.  It thus relies on a hybrid parallelisation scheme using both [OpenMP](https://openmp.org/wp/) and [MPI](https://www.open-mpi.org/). Postprocessing functions written in python (requiring [matplotlib](https://matplotlib.org/), [meson](https://mesonbuild.com/meson-python/) and [scipy](https://www.scipy.org/) are also provided to allow a useful data analysis.
 
 * **MagIC** is a free software. It can be used, modified and redistributed under the terms of the [GNU GPL v3 licence](https://www.gnu.org/licenses/gpl-3.0.en.html).
 
@@ -78,7 +78,7 @@ $ ./configure --enable-openmp --prefix=$HOME/local
 if FFTW is used or
 
 ```sh
-$ ./configure --enable-openmp --enable-ishioka --enable-magic-layout --prefix=$HOME/local --enable-mkl
+$ ./configure --enable-openmp --prefix=$HOME/local --enable-mkl
 ```
 
 if the MKL is used. Possible additional options may be required depending on the machine (check the website). Then compile and install the library
@@ -162,7 +162,7 @@ $ mpiexec -n 4 ./magic.exe input.nml
 
 ### 7) Data visualisation and postprocessing
 
-a) Set-up your PYTHON environment ([ipython](https://ipython.org/), [scipy](https://www.scipy.org/) and [matplotlib](https://matplotlib.org/) are needed)
+a) Set-up your PYTHON environment ([ipython](https://ipython.org/), [scipy](https://www.scipy.org/), [meson](https://mesonbuild.com/meson-python/) and [matplotlib](https://matplotlib.org/) are needed)
 
 b) Modify `magic.cfg` according to your machine in case the auto-configuration didn't work
 

bin/install-shtns.sh

@@ -6,7 +6,7 @@ if test ! -d $HOME/local; then
     mkdir $HOME/local
 fi
 
-wget https://bitbucket.org/nschaeff/shtns/downloads/shtns-$ver.tar.gz
+#wget https://bitbucket.org/nschaeff/shtns/downloads/shtns-$ver.tar.gz
 wget https://gricad-gitlab.univ-grenoble-alpes.fr/schaeffn/shtns/-/archive/v$ver/shtns-v$ver.tar.gz
 tar -xvf shtns-v$ver.tar.gz
 rm shtns-v$ver.tar.gz
@@ -44,7 +44,7 @@ else
    fi
 fi
 
-opts="--enable-magic-layout --prefix=$HOME/local"
+opts="--prefix=$HOME/local"
 if [[ -n $MKLROOT ]]
 then
    echo "MKL found, installing with MKL"

doc/magic_manual.pdf

@@ -91,7 +91,7 @@ the commonly used form
 
 .. math::
   \dfrac{\partial \vec{B}}{\partial t} = \vec{\nabla} \times \left( \vec{u}\times\vec{B} \right) 
-  + \lambda\,\vec{\Delta}\vec{B}.
+  + \lambda\,\vec{\nabla}^2\vec{B}.
   :label: eqInduction2
   
 The physical properties determining above equations are rotation rate
@@ -225,7 +225,7 @@ The first order Navier-Stokes equation (after to zero order hydrostatic referenc
    + \alpha \tilde{g}_o T' \dfrac{\vec{r}}{r_o}
    + \delta \tilde{g}_o \xi' \dfrac{\vec{r}}{r_o}
    + \dfrac{1}{\mu_0}(\vec{\nabla}\times\vec{B})\times\vec{B} 
-   + \tilde{\rho} \nu \Delta \vec{u}.
+   + \tilde{\rho} \nu \vec{nabla}^2 \vec{u}.
    :label: eqNSB
 
 Here :math:`u` and :math:`B` are understood as first order disturbances and
@@ -246,7 +246,7 @@ The first order energy equation becomes
 
 .. math::
   \tilde{\rho}\left(\dfrac{\partial T'}{\partial t}+\vec{u}\cdot \vec{\nabla} T' \right) = 
-  \kappa \Delta T'  + \epsilon,
+  \kappa \nabla^2 T'  + \epsilon,
   :label: eqEntropyB
 
 where we have assumed a homogeneous :math:`k` and neglected viscous and Ohmic
@@ -266,7 +266,7 @@ for chemical composition becomes
 
 .. math::
   \tilde{\rho}\left(\dfrac{\partial \xi'}{\partial t}+\vec{u}\cdot \vec{\nabla} \xi' \right) = 
-  \kappa_\xi \Delta \xi'  + \epsilon_\xi,
+  \kappa_\xi \nabla^2 \xi'  + \epsilon_\xi,
   :label: eqCompB
 
 where we have assumed a homogeneous :math:`k_\xi` and adjusted the definition of :math:`\epsilon_\xi`.
@@ -296,7 +296,7 @@ form:
    + \dfrac{Ra}{Pr} T' \dfrac{\vec{r}}{r_o}
    + \dfrac{Ra_\xi}{Sc} \xi' \dfrac{\vec{r}}{r_o}
    + \dfrac{1}{E Pm}(\vec{\nabla}\times\vec{B})\times\vec{B} 
-   + \Delta \vec{u},
+   + \vec{\nabla}^2 \vec{u},
    :label: eqNSBoussinesq
 
 where :math:`\vec{e}_z` is the unit vector in the direction of the rotation
@@ -471,6 +471,61 @@ Lorentz force in the Navier-Stokes equation. For small :math:`Di`  both heating
 terms can be neglected compared to entropy changes due to advection, an limit
 that is used in the Boussinesq approximation. 
 
+Equation for phase field model
+==============================
+
+To explore the uneven generation of topography associated with  the
+freezing and melting which occurs at the fluid-solid interface, MagIC can consider
+the evolution of motionless solid phase. To model the phase changes, MagIC relies
+on the phase field 
+formulation by `Beckermann et al. (1999) <https://doi.org/10.1006/jcph.1999.6323>`_ 
+combined with a volume penalization technique
+(`Hester et al. 2021 <https://doi.org/10.1016/j.jcp.2020.110043>`_).
+Practically, this method involves the
+time integration of a continuous scalar quantity :math:`\phi` which
+continuously varies from 0 in the liquid phase to 1 in the solid
+phase. A small dimensionless parameter :math:`\epsilon`, usually
+termed the Cahn number, then defines the ratio between
+the microscopic thickness of the transition between the two
+phases and the macroscopic domain size which is here
+the shell gap. Phase field methods represent a smoothed
+formulation of phase changes which are easier to implement
+numerically, especially when using pseudo-spectral methods,
+and converge to the exact moving boundary formulation
+in the limit of vanishing :math:`\epsilon`. The governing equation
+for the phase field reads
+
+.. math::
+   \frac{5}{6} St Pr \frac{\partial \phi}{\partial t} = a\nabla^2  \phi
+   -\frac{1}{\epsilon^2}\phi(1-\phi)[a(1-2\phi)+T-T_M],
+   :label: eqPhaseField
+
+In the above equation, :math:`T_M` denotes the melting temperature, while :math:`a`
+is an order one parameter of the phase field model which expresses the curvature
+dependence of the melting temperature `Beckermann et al. (1999) <https://doi.org/10.1006/jcph.1999.6323>`_.
+
+To account for the latent heat release during melting/freezing, an additional term
+needs to be added in the temperature equation. In the Boussinesq limit, this would
+yield (in absence of volumetric source terms)
+
+.. math::
+   \frac{\partial T}{\partial t}+\vec{u}\cdot\vec{\nabla} T = \frac{1}{Pr}\nabla^2 T
+   + St \dfrac{\partial \phi}{\partial t}.
+
+A penalty term of the form :math:`-\phi\vec{u}/(\epsilon^2\tau_p)` is also
+added to the r.h.s. of the Navier-Stokes equations :eq:`eqNSNd` to ensure that
+no flow is creeping into the solid. Practically, this term attenuates the
+velocity inside the solid phase, effectively treating it as a porous medium
+using a Darcy-Brinkman model.
+
+Although an optimal value for the parameter :math:`\tau_p` can be derived
+in simple configurations
+(`Hester et al. 2021 <https://doi.org/10.1016/j.jcp.2020.110043>`_),
+the value of :math:`\tau_p` needs to be
+adjusted for each simulation to ensure that the kinetic energy
+density of the solid phase remains smaller than that of the
+fluid phase.
+
 
 Dimensionless control parameters
 ================================
@@ -513,7 +568,15 @@ and the magnetic Prandtl number
    Pm = \frac{\nu_o}{\lambda_i}.
    :label: eqmaPrandtl
 
-In addition to these four numbers, the reference state is controlled by the geometry of
+In case, the phase field model is used, the Stefan number becomes also a control parameter defined by
+
+.. math::
+   St = \frac{{\cal L}}{c_p \Delta T},
+
+where :math:`{\cal L}` is the latent heat per unit mass associated with the solid-liquid
+transition.   
+
+In addition to these numbers, the reference state is controlled by the geometry of
 the spherical shell given by its radius ratio
 
 .. math::

doc/sphinx/equations.rst

@@ -17,8 +17,9 @@ several Implicit-Explicit time schemes where the nonlinear terms and the Corioli
 **MagIC** is written in Fortran and designed to be used on supercomputing
 clusters.  It thus relies on a hybrid parallelisation scheme using both `OpenMP
 <https://openmp.org/wp/>`_ and `MPI <https://www.open-mpi.org/>`_. Postprocessing
-functions written in python (requiring `matplotlib <https://matplotlib.org/>`_
-and `scipy <https://www.scipy.org/>`_) are also provided to allow a useful data
+functions written in python (requiring `matplotlib <https://matplotlib.org/>`_,
+`meson <https://mesonbuild.com/meson-python/>`_ and 
+`scipy <https://www.scipy.org/>`_) are also provided to allow a useful data
 analysis.
 
 .. figure:: figs/magic_occigen.png

doc/sphinx/foreword.rst

@@ -44,9 +44,10 @@ Welcome
 
      **MagIC** is written in Fortran and designed to be used on supercomputing
      clusters.  It thus relies on a hybrid parallelisation scheme using both `OpenMP
-     <http://openmp.org/wp/>`_ and `MPI <http://www.open-mpi.org/>`_. Postprocessing
-     functions written in python (requiring `matplotlib <http://matplotlib.org/>`_
-     and `scipy <http://www.scipy.org/>`_) are also provided to allow a useful data
+     <http://openmp.org/wp/>`_ and `MPI <https://www.open-mpi.org/>`_. Postprocessing
+     functions written in python (requiring `matplotlib <https://matplotlib.org/>`_,
+     `meson <https://mesonbuild.com/meson-python/>`_,
+     and `scipy <https://www.scipy.org/>`_) are also provided to allow a useful data
      analysis.
 
      .. figure:: figs/magic_occigen.png

doc/sphinx/index.rst

@@ -107,6 +107,15 @@ where :math:`u_h=(u_\theta^2+u_\phi^2)^{1/2}`, :math:`\delta h = \dfrac{r}{\sqrt
 
 where :math:`u_{F,r}` is the radial component of the fluid velocity and :math:`u_{A,r}=Br/\sqrt{E\,Pm}` is the radial Alven velocity. The denominator of the rightmost term accounts for the damping of the Alven waves.
 
+In case the phase field model is employed, the explicit treatment of the additional volume
+penalization term entering the Navier-Stokes equations yields an extra constraint
+on the time step size :math:`\delta t` compared to classical convection problems, such that
+
+.. math::
+   \delta t < 	C \tau_p \epsilon^2,
+
+where :math:`\epsilon` is the Cahn number.
+
 * **dtMax** (default :f:var:`dtMax=1e-4 <dtmax>`) is a  real. This is the maximum allowed time step :math:`\delta t`. If :math:`\delta t > \hbox{dtmax}`, the time step is decreased to at least dtMax (See routine `dt_courant`). Run is stopped if :math:`\delta t < \hbox{dtmin}` and :math:`\hbox{dtmin}=10^{-6}\,\hbox{dtmax}`.
 
 * **courfac** (default :f:var:`courfac=2.5 <courfac>`) is a real used to scale velocity in Courant criteria. This parameter corresponds to :math:`c_F` in the above equation.

doc/sphinx/inputNamelists/controlNamelist.rst

@@ -176,13 +176,13 @@ Phase field
   .. math::
      St = \frac{\mathcal{L}}{c_p\Delta T}
 
-* **tmelt** (default :f:var:`tmelt=0.0 <tmelt>`) is a real. This is the dimensionless melting temperature.
+* **tmelt** (default :f:var:`tmelt=0.0 <tmelt>`) is a real. This is the dimensionless melting temperature :math:`T_M` which enters :eq:`eqPhaseField`.
 
-* **epsPhase** (default :f:var:`epsPhase=0.01 <epsphase>`) is a real. This is the dimensionless interface thickness between the solid and the liquid phase (sometimes known as the Cahn number).
+* **epsPhase** (default :f:var:`epsPhase=0.01 <epsphase>`) is a real. This is the dimensionless interface thickness between the solid and the liquid phase (sometimes known as the Cahn number). This corresponds to :math:`\epsilon` in :eq:`eqPhaseField`.
 
-* **phaseDiffFac** (default :f:var:`phaseDiffFac=1.0 <phasedifffac>`) is a real. This is a coefficient that goes in front of the diffusion term in the phase field equation.
+* **phaseDiffFac** (default :f:var:`phaseDiffFac=1.0 <phasedifffac>`) is a real. This is a coefficient that goes in front of the diffusion term in the phase field equation. This corresponds to :math:`a` in :eq:`eqPhaseField`.
 
-* **penaltyFac** (default :f:var:`penaltyFac=1.0 <penaltyfac>`) is a real. This is coefficient used for the penalisation of the velocity field in the solid phase. The smaller the coefficient, the stronger the penalisation. Since this is a nonlinear term, it is handled explicitly and the time step size should be decreased with the square of :f:var:`penaltyfac`.
+* **penaltyFac** (default :f:var:`penaltyFac=1.0 <penaltyfac>`) is a real. This is coefficient :math:`\tau_p` used for the penalisation of the velocity field in the solid phase. The smaller the coefficient, the stronger the penalisation. Since this is a nonlinear term, it is handled explicitly and the time step size should be decreased with the square of :f:var:`penaltyfac`.
 
 
 Transport properties
@@ -259,6 +259,12 @@ Thermal diffusivity
   +----------------+----------------------------------------------------------------------------+
   | ``nVarDiff=4`` | polynomial-fit to an interior model of the Earth liquid core               |
   +----------------+----------------------------------------------------------------------------+
+  | ``nVarDiff=5`` | ...                                                                        |
+  +----------------+----------------------------------------------------------------------------+
+  | ``nVarDiff=6`` | tanh-like jump increase of thermal diffusivity in a stably-stratified-layer|
+  +----------------+----------------------------------------------------------------------------+
+  | ``nVarDiff=7`` | tanh-like jump increase of thermal diffusivity in a bottom stable layer    |
+  +----------------+----------------------------------------------------------------------------+
 
 .. _varnVarVisc:
 
@@ -277,6 +283,10 @@ Viscosity
   |                |    .. math:: \nu=\left(\frac{\tilde{\rho}(r)}                           |
   |                |              {\tilde{\rho}_i}\right)^{\alpha}                           |
   +----------------+-------------------------------------------------------------------------+
+  | ``nVarVisc=3`` | tanh-like increase of kinematic viscosity in a stably-stratified layer  |
+  +----------------+-------------------------------------------------------------------------+
+  | ``nVarVisc=4`` | tanh-like increase of kinematic viscosity in a bottom stable layer      |
+  +----------------+-------------------------------------------------------------------------+
 
   where :math:`\alpha` is an exponent set by the namelist input variable ``difExp``.
 
@@ -291,8 +301,6 @@ Anelastic liquid equations
 
   .. math:: \epsilon \simeq \frac{\Delta T}{T} \simeq \frac{\Delta s}{c_p}
 
-* **cmbHflux** (default :f:var:`cmbHflux=0.0 <cmbhflux>`) is a real. This is the CMB heat flux that enters the calculation of the reference state of the liquid core of the Earth, when the anelastic liquid approximation is employed.
-
 * **slopeStrat** (default :f:var:`slopeStrat=20.0 <slopestrat>`) is a real. This parameter controls the transition between the convective layer and the stably-stratified layer below the CMB.
 
 
@@ -379,7 +387,7 @@ Boundary conditions for chemical composition
 
   4. Imaginary amplitude (:math:`\sin` contribution)
 
-  For example, if the boundary condition should be a combination of an :math:`(\ell=1,m=0)` sherical harmonic with the amplitude 1 and an :math:`(\ell=2,m=1)` spherical harmonic with the amplitude (0.5,0.5) the respective namelist entry could read:
+  For example, if the boundary condition should be a combination of an :math:`(\ell=1,m=0)` spherical harmonic with the amplitude 1 and an :math:`(\ell=2,m=1)` spherical harmonic with the amplitude (0.5,0.5) the respective namelist entry could read:
 
 
   .. code-block:: fortran

doc/sphinx/inputNamelists/physNamelist.rst

@@ -85,13 +85,13 @@ or install it manually after downloading and extracting the latest version `here
 
 .. code-block:: bash
 
-   ./configure --enable-openmp --enable-ishioka --enable-magic-layout --prefix=$HOME/local
+   ./configure --enable-openmp --prefix=$HOME/local
 
 if FFTW is used or
 
 .. code-block:: bash
 
-   ./configure --enable-openmp --enable-ishioka --enable-magic-layout --prefix=$HOME/local --enable-mkl
+   ./configure --enable-openmp --prefix=$HOME/local --enable-mkl
 
 if the MKL is used. Possible additional options may be required depending on the machine (check the website). Then compile and install the library