[cig-commits] r19583 - seismo/2D/SPECFEM2D/trunk/src/specfem2D

[dkomati1 at geodynamics.org]

Author: dkomati1
Date: 2012-02-02 14:05:13 -0800 (Thu, 02 Feb 2012)
New Revision: 19583

Modified:
   seismo/2D/SPECFEM2D/trunk/src/specfem2D/specfem2D.F90
Log:
added a comment from Christina Morency about elastic / poroelastic coupling


Modified: seismo/2D/SPECFEM2D/trunk/src/specfem2D/specfem2D.F90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/src/specfem2D/specfem2D.F90	2012-02-02 22:03:54 UTC (rev 19582)
+++ seismo/2D/SPECFEM2D/trunk/src/specfem2D/specfem2D.F90	2012-02-02 22:05:13 UTC (rev 19583)
@@ -3761,7 +3761,7 @@
 if(ATTENUATION_PORO_FLUID_PART .and. any_poroelastic .and. &
 (time_stepping_scheme == 2.or. time_stepping_scheme == 3)) &
     stop 'RK and LDDRK time scheme not supported poroelastic simulations with attenuation'
- 
+
 if(coupled_elastic_poro) then
 
     if(ATTENUATION_VISCOELASTIC_SOLID .or. ATTENUATION_PORO_FLUID_PART) &
@@ -6565,6 +6565,42 @@
 !*******************************************************************************
 !         assembling the displacements on the elastic-poro boundaries
 !*******************************************************************************
+
+!
+! Explanation of the code below, from Christina Morency and Yang Luo, January 2012:
+!
+! Coupled elastic-poroelastic simulations imply continuity of traction and
+! displacement at the interface.
+! For the traction we pass on both sides n*(T + Te)/2 , that is, the average
+! between the total stress (from the poroelastic part) and the elastic stress.
+! For the displacement, we enforce its continuity in the assembling stage,
+! realizing that continuity of displacement correspond to the continuity of
+! the acceleration we have:
+!
+! accel_elastic = rmass_inverse_elastic * force_elastic
+! accels_poroelastic = rmass_s_inverse_poroelastic * force_poroelastic
+!
+! Therefore, continuity of acceleration gives
+!
+! accel = (force_elastic + force_poroelastic)/
+!     (1/rmass_inverse_elastic + 1/rmass_inverse_poroelastic)
+!
+! Then
+!
+! accel_elastic = accel
+! accels_poroelastic = accel
+! accelw_poroelastic = 0
+!
+! From there, the velocity and displacement are updated.
+! Note that force_elastic and force_poroelastic are the right hand sides of
+! the equations we solve, that is, the acceleration terms before the
+! division by the inverse of the mass matrices. This is why in the code below
+! we first need to recover the accelerations (which are then
+! the right hand sides forces) and the velocities before the update.
+!
+! This implementation highly helped stability especially with unstructured meshes.
+!
+
     if(coupled_elastic_poro) then
       icount(:)=ZERO