SOPALE : a finite-element computer code for the computation of visco-plastic creeping flows
with applications to tectonics
(short summary)





      Method, limitations, history



1.      What is SOPALE?


SOPALE is a computer code written in Fortran 77 which computes the flow of visco-plastic materials in the ‘creeping flow’ limit, that is, in the limit where viscosities are so high that inertial accelerations may be neglected. It is specialized to the field of geodynamics which studies the thermo-mechanical behavior of the crust, the lithosphere and upper mantle on geological time scales, for which it provides a forward evolution simulation tool.


2.      How does SOPALE work?


The material flow is represented as a sequence of quasi-static equilibria between body forces (due e.g. due the presence of an external uniform gravity field) and the viscous resistance opposed by internal stresses in the deforming material. The material deformation is caused by the combination of applied kinematical boundary conditions and body forces, and limited by the material viscosity. The equations of continuum mechanics are used to cast the equilibrium problem as a Stokes flow problem, which is a saddle-point system of partial differential equations to be solved for the velocity and pressure of the flow. The computed velocity field is then integrated in time to follow the trajectories of the material.


The code embeds a choice among a few parametric constitutive laws relating stresses and deformation (named rheological laws). These laws attempt to mimic experimental observations and depend on the category of material studied.

In SOPALE, they have been somewhat specialized, although in very generic form, to

give a coarse description of the mechanical properties of rocks (see section 5.). It should be emphasized that these rheological laws treat the geological material (rocks) as fluid and not solid. SOPALE and MICROFEM (see section 9) were probably among the first codes written in this perspective (for tectonic applications), which was chosen to allow very large deformation, while the traditional elastic or elasto-viscoplastic models were typically much more limited in the total deformation that could be computed.


Very large deformation is  one of the key characteristics of SOPALE.


In the context of geophysics, temperature is a primary factor affecting the value of viscosities, and the material flows in a non-uniform temperature field (named geotherm, for the earth). The Stokes model must be supplemented with a thermal model describing the conservation of heat. SOPALE uses the mechanical (Stokes) velocities as the transporting speed (one assumption of the physical model used)    and combines transport with heat diffusion, limited by the material thermal conductivity, and heat sources, representing for example the effect of radioactivity, or shear heating. Heat is therefore produced, diffused and transported, and the ratio of heat per thermal mass, or temperature field, by a feed-back effect on the rheologies, contributes to control the mechanical behavior.


This gives rise to a coupled ‘thermo-mechanical’ system.


Constitutive laws coded in SOPALE  for the mechanics and the thermics reflect the primary type of applications that SOPALE has been used for, that is problems in geophysics involving crustal and upper mantle dynamics over time scales typical of orogenesis (see section 3) . Similarly, all external processes (see section 6) applied to the evolving flow deal with specific aspects of tectonics or geomorphology.

The initial geometry of the domain to be deformed is typically a rectangular window, or a perturbation of such, with an aspect ratio selected as a function of the processes and scales considered.


The equations are transformed to a finite matrix form (in the sense of linear algebra), through finite-element discretization (FE for short).


The FE method approximates the material continuum by a grid (nodes and mesh) and uses a variational formulation of the equations of continuum mechanics. The FE method is applied to both the mechanical and thermal equilibrium models (see section 4)


SOPALE performs a finite sequence of incremental motion computations (or ‘time steps’) to move and deform the material. The velocities and temperature are computed with the FE method on a mildly deformed (eulerian) grid, while another (lagrangian) grid (or cloud of points) follows the material motion. The method is at least as old as Harlow and Welsh first particle-in-cell computations in Los Alamos (1965), and is widely used in computational fluid dynamics. One refers to such methods as Arbitrary-Lagrangian-Eulerian or ALE methods


-         an eulerian grid serves as the FE discretization grid, and carries the velocity and pressure computed from the Stokes equations. The same grid is used to solve the thermal equations. This grid conforms to the total domain boundary by adjusting its thickness (to match the current positions of the top and bottom boundaries)


-         A lagrangian grid or set of points (lagrangian  markers or particles) is moved along material trajectories. This grid records the physical history of material particles, in particular the pressure, temperature and accumulated deformation.

This grid is highly deformed and not suitable for a FE computation, which is why a separate eulerian grid is used.


To summarize, the computation initializes a domain geometry and spatial properties (fields), then computes a sequence of steps, each composed of the same subsequence: solution of a mechanical equilibrium (FE-Stokes), a solution of thermal ‘equilibrium’ (FE-Temperatures), and an update integrating in time the lagrangian or material fields, and computing a matching image of these fields on an update eulerian grid (ALE)



3.      What are the applications of SOPALE ?


The technology developed in SOPALE has been designed to solve problems for specific applications, and are therefore narrowly bound to them.  In this section we want to outline a few motivations which were central to the creation of the code and its development over the past 10 years.


SOPALE has been applied to the flow of geological materials (rocks, sediments) at crustal scale (around 30 kilometers) and upper mantle scale (depth of of order a few hundred kilometers). It has been used to study deformation subjected to compressive (e.g. fold and thrust belts) or extensive (rifting, lithospheric extension, necking) boundary conditions.


To briefly outline a few examples, SOPALE has been applied to


-understand the coupling between surface processes (e.g. erosion and deposition, in subaerial or submarine environments) and tectonics. This coupling is a direct result of the pressure sensitive nature of rock rheologies and leads to rock exhumation and metamorphic transformations which highly depend on the narrow interaction between the mechanical, thermal, and geomorphological components of the system (see section 9). The ability to represent such interactions in SOPALE is one feature which distinguishes it from many other public domain or commercial codes.


-understand halokinesis , in studies looking at the response of the gravitationally unstable (heavy plastic layer over a lighter, weak salt layer) to a given parametric progradation sequence.


-understand the mechanics of subduction and roll-back, at upper-mantle scale, for example the balance of  forces such as trench pull, slab push,  the influence of small scale convection, the role of crustal and mantle rheologies on slab break-off.


Applications have been developed primarily, but not exclusively (see section 9), at Dalhousie University in the Geodynamics Group (hereafter mentioned as DGG) led by Professor Chris Beaumont.


They have been published in a series of publicly available research articles. More details on these publications and the authors of various applications may be found on the Dalhousie Geodynamics Group website


Each application has typically enriched the SOPALE code with a module typical of the study at hand. Bonny LEE has implemented many of these features (see section 6 and 9), with the DGG geophysicists who would specify the type of law to be implemented.



4.      Details on the FE computations



·        Mechanical FE solution (Stokes flow)


The resolution of the Stokes equations uses the bilinear velocity/discontinuous pressure or Q1-P0 finite-element. This interpolation supplements the isoparametric 4-noded quadrilateral with a pressure degree-of-freedom representing the average pressure in the element. This ‘internal or bubble’ pressure unknown is eliminated by an approximation of the incompressibility condition named penalization. This approximation replaces the exact requirement of null velocity divergence by a linear relation between pressure and divergence, using a coefficient (the bulk viscosity) much larger than the shear viscosity.

The elimination of the pressure unknown replaces the original, saddle-point, Stokes problem, by an elliptical problem, which leads to a symmetric positive definite ‘stiffness’ or FE matrix. Such systems always admit a square root triangular matrix (the Cholesky factor, L) and can be solved, once L has been computed (Cholesky factorization), by 2 triangular matrix solves (upper and lower back-substitutions). This is the major benefit of the penalized approach, at the expense of approximate incompressibility, compared to the full indefinite solver in the velocity-pressure variables.


As viscosities are functions of the velocity gradient field (strain rate), through the rheological laws, the velocity is found by re-solving several times (non linear iteration loop) the equilibrium between external forces and internal (visco-plastic)

forces, with an update of viscosities from the last predicted strain rate to account for yield stresses and power-law rheologies (see section 5)


In each non-linear iteration:


       -the body forces are computed. They combine the effect of temperature on rock densities (coefficient of thermal expansion at constant pressure, reference density,reference temperature), and ,for some applications, the weight of a water column loading the top ‘mechanical’ surface.


         -parametrized  ‘phase change’ may occur during iterations, depending on local thermodynamic (pressure/temperature) conditions.


         -discontinuous pressures and nodal pressures are computed from the last predicted velocities. An input flag allows the filtering of a spurious 0-energy mode asscciated with the Q1-P0 element (a mode named checkerboard mode or ‘element-by-element oscillating mode’). The mode is removed by L2 projection which amounts to computing the bilinear field closest in square energy norm to the discontinuous ‘elemental’ pressure field. Note that other quasi-singular modes , with a energy dissipation rate tending to zero in the limit of 0 grid size, exist for the Q1-P0 element, but are not filtered or stabilized in this code.


          -viscosities are updated, using discontinuous pressures and other element-averaged quantities (see section 5) controlling the rheologies


          -the FE stiffness matrix of each element is computed and assembled in the form required by the sparse linear solver BLKFCT (see section 10) (the right-hand of the linear system, a column vector, is also assembled). The code uses selective integration, and integrates the spherical part of the matrix with the one-Gauss-point integration rule, and the deviatoric part of the matrix with a four-Gauss-point rule.


Boundary conditions are enforced by penalization: at prescribed degrees of freedom, the diagonal of the matrix is augmented by a term much larger than the diagonal maximum, and a reaction using the same lagrangian multiplier and the prescribed value is added to the right-hand side vector. Boundary conditions are not discussed in detail here.


·        Thermal FE solution


The thermal FE calculation uses the standard 4 noded quadrilateral element with an streamline-upwind Galerkin (SUPG) stabilization. In the limit where heat transport is faster than heat diffusion, that is in the high Peclet number range, the Galerkin variational method becomes inadequate for the solution of the diffusion-advection equation. The presence of a dominant first-order advecting term destabilizes the Galerkin method. We use the SUPG method proposed by Hughes (and now widespread in the FE community) and compute upwind positions of the 4 standard Gauss integration points. The points are moved by a function of the local Peclet number (1/tanh(P)-1/P) which solves exactly the 1d advection-diffusion problem with uniform advecting velocity. This is equivalent to a Petrov-Galerkin method and stabilizes computation at high Peclet number.


A Cranck-Nicholson scheme is used for time discretaization and the advection matrix term responsible for the lack of symmetry of the FE matrix is reported in the right-hand side. The system matrix is therefore a combination of mass and diffusion, is symmetric positive definite, and uses the same Cholesky solver as the mechanical system. The right-hand side combines the complementary Cranck-Nicholson term (mass and diffusion), the updated effect of advection and the heat sources (e.g. radio-active production or frictional heating). Even for linear thermal materials (constant conductivity, thermal inertia) iterations are required to solve the (linear) non-symmetric problem.



·        Linear solver : BLKFCT (see section 10)




5.      Details on constitutive laws


Code users may prescribe several material zones, which will be advected (i.e moved with the material motion). New material may also we added as a mass flux, e.g. by filing depressions left by the tectonic deformation in the topography (top surface shape).


Each material zone is assigned a ‘physical material’, which is characterized by a particular thermal or mechanical constitutive law.


Mechanical constitutive laws are all visco-plastic. They are characterized by a bulk viscosity (a numerical parameters enforcing the penalization of incompressibility)  and an effective viscosity which depends on various controls (e.g. temperature, effective pressure. total accumulated strain). The deviatoric stress tensor and deviatoric strain rate tensor are proportional, have the same principal axes, and the effective viscosity measures their ratio.  The visco-plastic material are assumed to behave viscously until the norm of the deviatoric stress reaches a yield envelope, at which point deformation becomes visco-plastic and adjusts to the stress yield level.


Plane strain assumptions are built-in the model (the model plane , typically vertical, is assumed to contain the minimum and maximum principal axis, no deformation occurs in the direction  perpendicular to the plane).


The viscous behavior sums (as in parallel rheological circuits) the forces from 3 creeping mechanisms :

non-linear or power-law creep, grain-size sensitive creep, and Peierls creep. The forces are driven by the ambient strain rate and the various controlling fields (see below). The plastic yield stress obeys a Drucker-Prager law.


These laws give crude approximations to the complex rheological behavior of rocks


Controls on the local rheology include:

-         various parametric (non dynamic) forms of the fluid pressure. Fluid pressure may be prescribed e.g. by a user-defined function of strain (interpreted as having an effect on porosity) multiplied by the fluid hydrostatic pressure.

Terzaghi’s effective stress principle is implemented in various forms (e.g. the Drucker-Prager yield stress is an affine function of effective pressure)

-         parametric user-defined dependence of cohesion and angle of friction on total strain (implying strain hardening or strain softening

-         temperature (through pre-exponential activation energy factors e.g. in thermally activated dislocation creep) or pressure (activation volume)

-         the total accumulated strain (which integrates in time the lagrangian strain rate invariant, and should be considered as a measure of dissipated work, or of accumulated damage, rather than a geometrical deformation measure.

-         grain size




6.      Details on external controls


Users of SOPALE can provide information relative to the non-modeled part of the geological/geophysical system, by prescribing various forces, mass, or geometric loads in the form of imposed functions of time and space. The code response to these inputs needs to be checked and interpreted, as no consistency is required for these parameters.


Users may ,for example :


·        Subject the model to surface processes, which will either redistribute mass in a diffusive (curvature-dependent), or slope-dependent manner, or add mass to the system by partial or total filling of accommodation space below a reference height. (e.g. to mimic basin fill). The added material (sediment) may then be reworked by tectonics and lead to structures similar to e.g. inverted basins.

Sediment addition may take the form of uniform agradtion or a gaussian-shape progradation (progressing seaward)

·        Subject the model to a waterload, e.g. at an oceanic margin. This will later the surface processes (subaerial versus submarine), the fluid hydrostatic pressure and the total weight applied to the model supporting base

·        Subject the model to various boundary kinematical boundary conditions (velocities, with of without mass influx) and thermal boundary conditions (given temperature or heat flux)

·        Provide an external support to the model in the form of a broken beam, made of 2 elastic members meeting at a common point. This allows users to rest the model on an underlying elastic support (e.g. a piece of mantle) floating in an inviscid fluid. Models combining Airy and flexural isostasy are therefore available, and various rotation/deflection/force/torque boundary conditions may be applied to the beam ends. No opening is allowed however between the beams, which are assumed to join continuously. Flexural rigidities may be assigned per column.



7.      Limitations of SOPALE


What SOPALE can do is best documented in the published research papers using it (see section 3). We want here to explain what the code cannot do or is not guaranteed to do, and why. This will help preventing a misuse of the code, warning code users, and circumscribing the limit of the technology, which might be useful in issues involving distinctions between code improvements and new technologies.


We see 3 major sources of limitations to the present technology: mathematical model limitations, physical model limitations and algorithmic model limitations:


·        Mathematical model limitations


The rheologies used in SOPALE belong to the family of non-associated plastic models, for which the plastic potential and plastic yield functions differ. The standard tools of mathematical analysis (variational inequations and convex analysis) fail to provide a method for proving the existence and unicity of solutions in this case. The numerical model used in SOPALE can only be

seen as a pragmatic extension of a technique fully proven for other (e.g. linear, or associative plastic) rheologies, without much (if any) theoretical support/justification for the computed solution for the general non-linear case.


 The mathematical problem at hand is very challenging for it combines incompressibility, wide range of viscosities, and a yield stress which becomes null (singular) at the surface (for Coulomb plastic materials which have no cohesion). Most (all?) existing codes adopting the same method, face the same problem. This problem is aggravated by a lack of coherence in the formulation itself, i.e. the continuum mechanics description is unable to provide a consistent definition for faults (stress and velocity characteristics differ).


SOPALE can be used to approximately model  rigid-plastic systems. When the system is built with either one physical zone operating in the viscous range, or with a mix of viscous and plastic layers, the viscous part of the system ‘gives a time scale ‘ for evolution and  constrains the system to have a well defined solution. Hence we discuss here the more difficult question of rigid-plastic systems.


The traditional approach of solution uses a system of hyperbolic equations (Hilda Geiringer’s famous equations) which is not a boundary value problem, but a propagation problem. The theoretical question is: what validity does the plasticity penalization, boundary value method, used in SOPALE (at the continuum level) have for such problems? We note that rigid-plasticity  may be seen as the limit of a Bingham fluid rheology when the fluid viscosity is null. Bingham fluids use the Von Mises potential as a minimum stress below which no flow occurs. Duvaut and Lions have proven the existence and unicity of the interior Bingham fluid evolution problem, using the method of variational inequalities. However in SOPALE rigidity is not modeled precisely : it is approximated (regularized) by a user-chosen large viscosity below yield.


Known solutions for nonlinear plastic problems are rare and hard to find.  This prevents the calibration of numerical errors for problems involving these rheologies, that is for quasi-all problems in the range of applications of SOPALE.


It is therefore impossible to provide an estimation of numerical errors, and the code provides no indication of accuracy (even if it has of course been tested in a series of simple benchmark problems: Rayleigh-Taylor, thermal convection, critical wedge average slope , punching problem, advection-diffusion, simple corner flows…)


SOPALE can only be used as a qualitative tool indicating the trends of the system ‘s evolution, or identifying the major controls of this evolution.



·        Physical model limitations


We have selected in this section a few examples showing some limitations of the physical model used in the code. These examples are far from exhaustive, and given only for illustration.


The code uses a plane strain assumption and models a two-dimensional vertical plane, without deformation in the direction orthogonal to the plane. It ignores the earth curvature even on the scale of 500-1000 kilometers. It isolates a ‘slice’ of material which is part of a larger system, therefore requires external assumptions

regarding the fluxes/forces applied to the modeled ‘physical window’.


Consequently, many external parameters are introduced (see section 6), without any restriction imposed by a system of equations (these ad-hoc parameters will give opportunities to users to bias/modify the results in various directions).


This implies that good users should be ‘field-specialists’


They should be aware of the meaningful ‘forcing’ time scales and spatial scales (e.g. in surface processes), and capable of interpreting the pertinence of numerical results.


Rheologies involve the use of the effective pressure supported by the rock matrix However, fluid pressures are parametrized without an explicit solution of the equations of Darcy, (involving explicit porosity, permeability model etc.)


Rheologies are fluid and do not capture the elastic behavior of rocks.


Material phase changes are allowed but there is no associated production or absorption of heat..


Grain size sensitive creep rheologies are available but no model of grain size evolution is available.


All models are isotropic, and remain so even for material which have been highly sheared along a constant plane direction (it is questionable that this material would not develop an orthotropic symmetry, even at large scale).


The role of hydration / dehydration is not taken into account in any rheology.

However this role is important:


-in the dynamics of subduction

 The water cycle is perceived to play a key role in the lubrification of plate tectonics and dehydration of the water-loaded subducting plate, or dehydration related to partial melting is believed to have a major role on viscosities and on the time of transport of diapiric intrusions


-in halokinesis: the rheology of evaporates and salt rocks is sensitive to water and the cycle of water may play a role in the timing of salt diapirism. We also note that no sediment lithification modeled is included in the code. Other research work (e.g. Tuncay and Ortoleva) have attempted to use ‘finer scale’ rheological models which couple the physics and chemistry of sediment diagenesis with the evolving (dynamic) rheology of the sediments. Quantitative timing of halokinesis may be a very tricky business, because the Rayleigh-Taylor instability will amplify any model deficiency!


-erosion, deposition, progradation are all geometric postprocessors of the surface shape. They typically do not conserve mass and there is no attempt to understand where material is eroded from or where it is transported. Chronostratigraphy is imposed as a user-defined sequence of sedimentation rates and material labelling corresponding to a predetermined cycle, irrespective of the evolution of the model. These postprocessors are therefore different from numerical models (such as Cascade) used in the Geophysics community to compute the evolution of the surface geomorphology.



Hence, again, computations are rather generic in the code, and should only be seen as qualitative.



·        Algorithmic model limitations


-limitations of the discretization used :


The code uses a very primitive mesher, which only allows for a tensor product (given number of rows and columns) rectangular mesh. The top and bottom mesh boundary may be initially perturbed. They follow at each time the domain boundary, as prescribed by their net motion (e.g. tectonic minus erosional). This is done by ‘thickening’ or thinning elements. This may lead to a loss of accuracy (element distorsion).


The eulerian grid only moves its nodes vertically (by thickening) to follow the outer shape.

There is no addition of nodes, no adaptation of the grid where gradients would be higher, or near physical material boundaries.

Thickening involves a loss of resolution due to discretisation of a greater height column by a constant number of elements. Fortunately, in crustal-scale tectonic applications, angles are quiet small the thickening bounded by a factor around 2. However in other models the loss of accuracy does become very important.


Another deficient aspect of the discretization is the lack of  any form of discretization for material zones boundaries, i.e. the limit between a physical material and another physical material (e.g. a layer of salt and a layer of sediments) is simply not represented. Each material zone is represented by an array of particles and  patching is used to avoid cases where an eulerian grid would not contain any lagrangian particle. Much better algorithms (e.g. using a lagrangian Delaunay grid) clearly exist and have been tried, but not in SOPALE.


Some eulerian finite-elements may contain particles belonging to different material zones. In this case a crude majority rule is used. This rule does not take into account any of the actual flow stress that the participant materials would actually compute !



-treatment of the coupling between temperature and rheologies


The code never solves a coupled equilibrium involving velocities, pressure and temperature. It simply alternates a mechanical and a thermal solution. A fully coupled treatment, or an iterative solution of this coupling has been implemented by other experts in the field (e.g. Peter Van Keken) and would compute a more accurate equilibrium


-limitation of the linear solver used

Better, faster sparse solver exist. If a parallel hardware platform can be used, the model spatial resolution may be increased (FE grids of 1000x1000, or equivalent) should be feasible on  32 processors machines (e.g. IBM  p690 or Linux clusters)


            This increased resolution may be used e.g. to include other physical scales in the model and/or to increase or test the model accuracy




One of the ‘co-inventors’ (PF) thinks that new technologies and continued research work will be needed to address most of the limitations mentioned in this section. 


8.      Where is the code documentation?


In addition to the present document, users of SOPALE may read the following manuals :

  1. The sofware manual explains the code flowchart, with an emphasis on the concepts used in SOPALE.
  2. The source manual explains how the source code is organized and discusses the static memory allocation of variables
  3. The theory manual explains the numerical techniques used.
  4. The user manual explains how to use the software (inputs,run,postprocessing)

Documents 1 to 3 have been written by Philippe Fullsack and document 4 by Bonny Lee and Mai Nguyen.



9.      What is the code history? Who contributed what to the present code?


SOPALE was written by Philippe Fullsack (PF)  as a simplification of a more ambitious parallel FE code (OPALE) based on the same principles. OPALE was itself a more sophisticated version of a precursor code (MICROFEM) but was fully parallel (while SOPALE is serial) . Hence SOPALE stands for simplified-OPALE, and OPALE stands for OPtimized ALE code (optimized in the sense of a more powerful, parallel code).  SOPALE owes its success to its simplicity, and to the use of a (then) state-of-the-art sparse Cholesky solver (see below). OPALE and SOPALE were developed in the years 1995-1196, 10 years ago.


In DGG, MICROFEM was used mostly for collaborative work between Rebecca Jamieson, a metamorphic petrologist, and Chris Beaumont. The code got specialized to computations of pressure-temperature-paths, by prescribing a subduction geometry and embedding the mechanical,crustal scale problem into a larger, subduction-scale, thermal problem. This specialization also rang the knell of MICROFEM. PF deliberately chose to develop OPALE and SOPALE in a more generic form, starting with models which would involve a single eulerian grid, (and exploring problems with weaker kinematical controls). This gave a second chance to the method.


This example, chosen to illustrate some of the driving motivations for SOPALE, also shows how building some applications, in SOPALE, may require specific modifications


Similar computer codes applying the concepts of continuum mechanics to tectonic problems were developed in DGG by Jean BRAUN (visco-elastic code NONSAP) and Sean WILLETT (who used a fluid  Taylor-Hood FE method, which was proven to be inaccurate for coulomb plastic calculations! )


Bonny Lee has clearly contributed the most to the maintenance, documentation and coordination of code developments over the past 10 years. She has been implementing on a regular basis modifications requested by Chris Beaumont and absorbed modifications into the reference version.


Erik Demaine has been working with DGG. He developed an X11 graphic package (MOZART) that was used by SOPALE users years ago.



Over the years, other post-processors or visualization tools have been developed  by DGG members to exploit and interpret the results (e.g. IDL scripts by Sergei Medvedev and Mai Nguyen, or exports to MATLAB)


No post-processing tool is included in the technology discussed here. However a commercialization of SOPALE owe to address the question of outputs and post-processing. Indeed the binary output files produced by the code are of no direct use and need further processing for the interpretation of results.  **


Erik worked with PF on developing in-house codes for parallel Cholesky calculations, and suggested using the Cholesky solver BLKFCT written by Barry Peyton and Esmond G. Ng  

As far as PF remembers, Erik got the Fortran code directly from them, with no need for a licence, in 1994 or 1995 (?).  (see section 10)



Ritske Huismans, Susan Buiter, Steve Ings, Lykke Gemmer, all have worked with, debugged, modified, SOPALE or contributed to interfaces with other packages.


Steve Ings, Lykke Gemmer have tested the code for the problem of incipient normal faulting of a plastic sediment overburden above a weak viscous salt layer. They have designed models combining progradation with Rayleigh-Taylor instability.


Ritske Huismans has contributed to the deposition model, addition of new particles, record of stratigraphy, corrections of bugs in the thermal calculation.


Susan Buiter has contributed to sandbox extension benchmarks and compared SOPALE to other codes with similar functionalities, particularly in view of the brittle deformation pattern computed. This work is associated to the ongoing project of development of a long-term crustal dynamics software (GALE) which would implement principles and methods similar to SOPALE’s in the Computational Infrastructure for Geodynamics project .She created a copy of the beam code written by Philippe Fullsack for microfem  for use in SOPALE.


Russ Pysklywec has spent time testing the thermo-mechanical coupling in simple Boussinesq thermal convection experiments with constant and variable viscosity profiles.


Susan Ellis has been associated with the precursor MICROFEM code, and with all concepts involved in ALE visco-plastic fluid calculations. She may have used (?) SOPALE (Bern Institute of Geology) with Adrian Pfiffner




Christina Morency is developing now a version of SOPALE involving a coupling of the Stokes equations with the equations of Darcy (poro-visco-plasticity).This is placed  in a comment bracket because her modifications are still at the development stage, and not yet part of the reference version maintained by Bonny Lee


‘Real life users’, who attempt to use the code to tackle the understanding of actual geological/tectonic settings, have included, e.g.:


-for the ROCKIES :                     Glen Stockmal

-for the PYRENEES :                  Josep Munoz

-for the HIMALAYA :                 Rebecca Jamieson

-for the SWISS ALPS :                Adrian Pfiffner


No attempt is made here to analyze the code potential or deficiencies in these applications, which is best left to direct discussions with users.


Code users have received lots of help for running models or post-processing results from Sergei Medvedev and Mai Nguyen.





Probably all users in the following list have contributed to some development, concept, testing, or exploration of new applications for SOPALE. 



The order of the list is arbitrary and we leave to the assessors of the technology the role of determining further each individual’s contribution to the present technology.

Links and e-mail addresses are provided for contacts.






10.  External libraries used by SOPALE



                             The supernodal sparse Cholesky solver BLKFCT


Esmond G. Ng and Barry W. Peyton wrote in collaboration with J.W.H. Liu a supernodal sparse Cholesky solver, when they both were working at Oak Ridge National Laboratory (ORNL). Liu contributed the multiple mimimum degree reordering algorithm used in their package)


We were told (private correspondence with the authors) that we could use the name BLKFCT, but have also seen this solver referred to as the ORNL sparse matrix solver.  


This solver was adapted to the type of problem we solve (sparse symmetric positive definite systems) and addressed the question of memory hierarchy by allowing tuning of parameters such as cache-size.

This however proved less useful on our computers than the efficiency of supernodal factorization which ‘vectorizes’ the supernode factorization  and allows better computational intensity, as defined by the ratio of floating points computation per memory access.


The same solver is used 10 years later; and it has extermely useful to SOPALE and its users.



** SOPALE uses their sequential Cholesky solver but they also have developed more other efficient linear algebra tools to facilitate various applications **


However, we ignore if BLKFCT can be packaged in a distribution of SOPALE under licence (commercial or else). Any non-academic use of BLKFCT should be discussed with the authors mentioned above!


**We wish to acknowledge the gracious donation of BLFCT for public domain use. Such tools simply make research feasible. BLKFCT is used in some public domain packages (e.g. the interior point package LIPSOL  under Gnu general public licence, or SCILAB, another free, public domain software )**




** Esmond G. Ng ( short biography )now works for the US Department of Energy at the Lawrence Berkeley National Laboratory in San Francisco, CA, where he leads the scientific computing group (SCG). He is a partner in the Scientific Discovery Through Advanced Computing program ( )


It is worth noting that Esmond Ng has been working for on the solution of Stokes's pressure system within N3S using supernodal Cholesky factorization. Philippe Fullsack worked on iterative methods in the same Navier-Stokes code (N3S) , before coming to Canada in 1988. This gives credit to direct solvers as a viable alternative to iterative Krylov solvers**




                                                  SOPALE has also borrowed:


-         from the Fortran book of Numerical Recipes published by Cambridge University Press  (William H. Press and Saul A. Teukolsky)  ( Numerical Recipes ) the cubic spline interpolation routines  (num_receip_spline, num_recip_splint are copies of the original : spline routines published in this book)

-         from the ‘FE bible’ (The finite-element method by O.C. Zienkiewicz and Robert L. Taylor ) the routines for the computation of FE shape functions and gauss integration (pgauss1 and shape1 are simplified copies of the code given in appendix written by Robert Taylor. His code eventually evolved to the public domain FEAP program)

-        from the  MICROFEM code, concepts and the FE beam deflection routine.


11.  Similar projects, and a few experts


** Experts in the field of visco-plastic geophysical flows: Uli Christensen, Yuri Podaldchikov, Peter Van Keken, Louis Moresi, Jean Braun, Haro Schmeling


FE codes similar to SOPALE (fluid, ALE, visco-palstic, have been developed by Louis Moresi (CITCOM, ELLIPSIS), Boris Kaus SloMo, Peter Van Thienen, Peter Van Keken, Alexei Poliakov and Yuri Podlachikov, Harlow and Welsh …


Mathematicians experts include J.L. Lions, R.Duvaut , D. Reddy, D. Griffiths, D. Silvester ..**


document written by Philippe Fullsack (PF) ( March 26,27 2006)