Benchmarks for the Magma Dynamics Demonstration Suite (MADDs).

A new formulation for the equations of magma migration in viscous materials as originally derived by McKenzie is presented, as well as a set of well-understood special case problems that form a useful benchmark-suite for developing and testing new codes. upload file

The Magma Development team has finished the alpha release of the Magma Dynamics Demonstration Suite (MADDs). The initial code implements the zero porosity/no melting magma benchmark for mid-ocean ridge solid flows in 2D and 3D built on the Underworld framework. The purpose of this code is principally to validate accurate pressure solvers for Stokes flow in current CIG supported software. The stgMADDs source code is available in CIG's Mercurial Repository (geodynamics.org/hg).

For a first time download of the stgMADDs repository, do the following:

Create the topmost repository with:

   hg clone http://geodynamics.org/hg/magma/3D/stgMADDs

Then obtain all the other repositories using:

   ./obtainRepositories.py 

To push, you may have to use the ssh syntax, e.g.:

    hg push ssh://hg@geodynamics.org/hg/magma/3D/stgMADDs

Caveat Emptor: This is very much an alpha release code for experimentation with the accuracy of different mixed FEM pressure solvers. Questions, complaints and bug reports should be directed to cig-magma@geodynamics.org.

Details how to run the first milestone of the MADDs project in 2D and 3D and provides some results of these simulations. It also gives the rates of convergence of the pressure gradient solutions as the resolution is increased.

In order to run the simulations for milestone 1 of the MADDs project (in 2D), first:

   cd Magma/Models/Milestone1/Ridge2D_Quadratic

Then make a symbolic link to the executable binary (assuming the code has been successfully built):

   ln -s ../../../../build/bin/StGermain . 

The simulation may then be run (in parallel), passing in the respective XML file as input:

   mpiexec -np <# of procs> ./StGermain Ridge2D.xml

Alternatively, the 3D simulation may be run as:

   cd Magma/Models/Milestone1/Ridge3D_Quadratic
   ln -s ../../../../build/bin/StGermain .
   mpiexec -np <# of procs.> ./StGermain Ridge3D.xml

These simulations will produce graphical output of the velocity, pressure and pressure gradient solutions, as well as the analytic reference solutions and the element-wise normalised L2 error fields for the pressure and pressure gradients, as shown below. It will also generate text files to the output directory giving the node-wise results for the respective fields.

As the resolution is increased, the normalised global L2 errors are observed to decrease. This decrease is approximately linear for the 2D ridge model and slightly poorer for the 3D model. Graphs detailing the global errors as a function of resolution are given below.

Velocity, pressure and pressure gradients solutions and L2 errors for a 2D ridge model with 120 x 60 elements.

Global Pressure Gradient Errors for 2D Ridge Model. Normalised global L2 errors.

Velocity, pressure and pressure gradient solutions and L2 error fields for 3D ridge model.

Global Pressure Gradient Errors for 3D Ridge Model. Global normalised L2 pressure gradient errors at varying resolutions.

Details the results of the Milestone 2 simulations and analyzes the accuracy of the advection scheme.

The second milestone contains several different models. In each a porosity distribution has been supplied, which acts as a force term on the stokes equation. The porosity is advected according to a semi Lagrangian scheme, with no natural diffusion. Details of these models are given below.

The first model is a simple demonstration of how the porosity dependent force term drives the porosity distribution up through the domain. The domain is subjected to free slip boundary conditions, which also distort the porosity distribution as it is advected. This model is housed in the directory:

   Magma/Models/Milestone2/ZPNM_SemiLagrangianPorosity/

with the XML file to be passed at run time being the Porridge.xml file.

upload image - Gaussian Porosity Field Advection. Advection of Gaussian porosity field as a Stokes equation force term. The lower density porosity region rises due to gravity.

The second model is an extension of the first milestone, that is a ridge model with the initial boundary conditions loaded from a reference solution. However in this case a Gaussian porosity distribution has been added. The porosity distribution is distorted as it moves up through the domain in accordance with the velocity field generated from the boundary conditions. This model can be found in:

   Magma/Models/Milestone2/RidgeModelWithGaussianPorosity/

with the XML file to be passed being Ridge2D.xml.

upload image - Ridge Model with Gaussian Porosity Field. Stokes flow with 2D ridge model boundary conditions and Gaussian porosity initial distribution, driven by a porosity dependent force term.

As well as the models, a test is also supplied for validating the accuracy of the semi Lagrangian advection scheme. This involves an initial porosity distribution (either Gaussian or a diagonal line step function, as determined from the XML), which is subjected to a static shearing velocity field. The porosity distribution is subjected to the velocity field for a finite number of time steps (as determined from the XML input file), before the velocity field is reversed and the distribution advected back for the same number of time steps. The normalised global L2 error between the initial and final distributions is then calculated. This test is housed in the directory:

   Magma/Models/Milestone2/tests/

And may be run using the input file testSemiLagrangianAdvection.xml.

Running this simulation at varying resolutions, the convergence of the errors for the advection scheme were determined using both the Gaussian distribution and diagonal step function as initial conditions. The errors were recorded for schemes which used a cubic method as well as a quadratic method based on the element shape functions for interpolating the value at the end point of the characteristic. As can be observed from the convergence errors plots below, the quadratic and cubic interpolation method schemes converged at comparable (less than linear) rates using the step function initial condition, with an improvement in those results obtained using the cubic interpolation method. When the Gaussian initial distribution was applied (which was easier to solve accurately on account of the smoother gradients involved), both interpolation methods converged at a rate much better than linear, with the cubic interpolation method proving far superior to the quadratic method.

upload file - Semi Lagrangian Advection Scheme Test - Step Function. Diagonal step function initial distribution subjected to a shearing velocity field.

upload file - Semi Lagrangian Advection Scheme Test - Gaussian Distribution. Gaussian initial distribution subjected to a shearing velocity field.

upload file - Error Convergence for Advection Scheme - Step Function IC. Normalised global L2 errors for semi Lagrangian advection scheme with a diagonal step function initial condition as a function of resolution.

upload file - Error Convergence for Advection Scheme - Gaussian IC. Normalised global L2 errors for semi Lagrangian advection scheme with Gaussian initial distribution as a function of resolution.

Milestone 3 Results Details the results for the third milestone, in which the melt velocity was determined given the existing solid velocity and pressure fields. Melt Model - 2D Ridge with Constant Porosity Solid and melt velocity, pressure and pressure gradient fields for 2D ridge model with constant porosity. Melt velocity magnitudes are significantly larger near the point of discontinuity due to their proportionality to the pressure gradients, which are largest at these points. Melt Model - Gaussian Porosity Driven Flow Solid and melt velocity, pressure and pressure gradient fields for Stokes flow driven by a Gaussian initial porosity distribution. Milestone 4 Results and Analysis Discussion of the system being modeled, and details of how to run the model with different initial conditions in 2 and 3D. 2D Solitary Wave A 2D solitary wave with a wave speed of 7 rising through a solid with a constant speed of -2. The wave shows no visible diffusive behavior. Noisy 1D Solitary Wave Initial Condition Initial condition of a vertically changing 1D solitary wave with a certain amount of introduced noise, which allows 2D solitary waves to emerge over time. Emerging 2D Solitary Waves Solitary waves emerging from a noisy 1D solitary wave initial condition. Emergent 2D Solitary Waves Solitary waves having emerged from a noisy 1D solitary wave initial distribution. Emergent 3D Solitary Waves. 3D Solitary Waves emerging from a noisy 1D Solitary Wave initial distribution Milestone 5 Results and Analysis Results and analysis for the isoviscous McKenzie equations (with melting) driven by a corner flow velocity BC. Isoviscous McKenzie System with Corner Flow BC - 1 After 1 time step Isoviscous McKenzie System with Corner Flow BC - 50 After 50 time steps Isoviscous McKenzie System with Corner Flow BC - 3200 After 3200 time steps. Velocity - x component x component of the velocity field for the 3D isoviscous McKenzie model with ridge BCs at time step 150. Velocity - y component y component of the velocity field for the 3D isoviscous McKenzie model with ridge BCs at time step 150. Dynamic (Stokes) pressure Dynamic pressure due to viscous shear for the 3D isoviscous McKenzie model with ridge BCs at time step 150. Porosity Porosity field for the 3D isoviscous McKenzie model with ridge BCs at time step 150. Compaction pressure Compaction pressure due to compressibility of the solid phase for the 3D isoviscous McKenzie model with ridge BCs at time step 150. Melt fraction Melt fraction field representing the melt to solid phase of the 3D isoviscous McKenzie model with ridge BCs at time step 150. Melt velocity - x component x component of the melt velocity field for the 3D isoviscous McKenzie model with ridge BCs at time step 150. Melt velocity - y component y component of the Melt velocity field for the 3D isoviscous McKenzie model with ridge BCs at time step 150.