[cig-commits] r8423 - seismo/2D/SPECFEM2D/trunk/SPECFEM90

[walter at geodynamics.org]

Author: walter
Date: 2007-12-07 15:45:04 -0800 (Fri, 07 Dec 2007)
New Revision: 8423

Modified:
   seismo/2D/SPECFEM2D/trunk/SPECFEM90/constants.h
   seismo/2D/SPECFEM2D/trunk/SPECFEM90/specfem2D.f90
Log:
added threshold to check that run remains stable in FD and SEM codes.
also added flexible source location and orientation in FD code.


Modified: seismo/2D/SPECFEM2D/trunk/SPECFEM90/constants.h
===================================================================
--- seismo/2D/SPECFEM2D/trunk/SPECFEM90/constants.h	2004-05-29 12:16:02 UTC (rev 8422)
+++ seismo/2D/SPECFEM2D/trunk/SPECFEM90/constants.h	2007-12-07 23:45:04 UTC (rev 8423)
@@ -9,6 +9,9 @@
 ! mesh tolerance for fast global numbering
   double precision, parameter :: SMALLVALTOL = 0.000001d0
 
+! displacement threshold above which we consider the code became unstable
+  double precision, parameter :: STABILITY_THRESHOLD = 1.d+25
+
 ! input and output files
   integer, parameter :: IIN  = 40
 

Modified: seismo/2D/SPECFEM2D/trunk/SPECFEM90/specfem2D.f90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/SPECFEM90/specfem2D.f90	2004-05-29 12:16:02 UTC (rev 8422)
+++ seismo/2D/SPECFEM2D/trunk/SPECFEM90/specfem2D.f90	2007-12-07 23:45:04 UTC (rev 8423)
@@ -103,7 +103,7 @@
   integer ie,k,isource,iexplo
 
   integer ielems,iglobsource
-  double precision f0,t0,factor,a,angleforce,ricker
+  double precision f0,t0,factor,a,angleforce,ricker,displnorm_all
 
   double precision rsizemin,rsizemax,cpoverdxmin,cpoverdxmax, &
     lambdalSmin,lambdalSmax,lambdalPmin,lambdalPmax,vpmin,vpmax
@@ -970,10 +970,13 @@
 ! `corrector' update velocity
   veloc(:,:) = veloc(:,:) + deltatover2*accel(:,:)
 
-!
 !----  display max of norm of displacement
-!
-  if(mod(it,itaff) == 0) print *,'Max norm of displacement = ',maxval(sqrt(displ(1,:)**2 + displ(2,:)**2))
+  if(mod(it,itaff) == 0) then
+    displnorm_all = maxval(sqrt(displ(1,:)**2 + displ(2,:)**2))
+    print *,'Max norm of displacement = ',displnorm_all
+! check stability of the code, exit if unstable
+    if(displnorm_all > STABILITY_THRESHOLD) stop 'code became unstable and blew up'
+  endif
 
 ! store the seismograms
   do irec=1,nrec