[cig-commits] r19359 - in seismo/3D/FAULT_SOURCE/branches/new_fault_db: . utils

[ampuero at geodynamics.org]

Author: ampuero
Date: 2012-01-13 17:49:16 -0800 (Fri, 13 Jan 2012)
New Revision: 19359

Added:
   seismo/3D/FAULT_SOURCE/branches/new_fault_db/utils/
   seismo/3D/FAULT_SOURCE/branches/new_fault_db/utils/critical_timestep.m
Modified:
   seismo/3D/FAULT_SOURCE/branches/new_fault_db/README_SPECFEM3D_FAULT
Log:
doc Kelvin Voigt viscosity

Modified: seismo/3D/FAULT_SOURCE/branches/new_fault_db/README_SPECFEM3D_FAULT
===================================================================
--- seismo/3D/FAULT_SOURCE/branches/new_fault_db/README_SPECFEM3D_FAULT	2012-01-13 22:57:20 UTC (rev 19358)
+++ seismo/3D/FAULT_SOURCE/branches/new_fault_db/README_SPECFEM3D_FAULT	2012-01-14 01:49:16 UTC (rev 19359)
@@ -187,7 +187,7 @@
 	has a strict format:
 
 	Line 1: Number of faults (NF)
-	Lines 2 to NF+1: Kelvin Voigt damping (in seconds) for each fault.
+	Lines 2 to NF+1: Kelvin Voigt damping (in seconds) for each fault. (See below how to set this parameter)
         Line NF+2: Type of simulation (1=dynamic , 2 = kinematic)       
 	Line NF+3: Number of time steps between updates of the time series outputs at selected
 		fault points (see DATA/FAULT/FAULT_STATIONS), usually a large number (100s or 1000s). 
@@ -262,7 +262,31 @@
                 sets a constant value (val) within a sphere with center (xc,yc,zc) and radius r. 
 
 
+Note: How to set the Kelvin-Voigt damping parameter: 
+      The purpose of the Kelvin-Voigt viscosity is to damp spurious oscillations
+      generated by the fault slip at frequencies that are too high to be resolved by the mesh.
+      The viscosity "eta" (in seconds) depends on the size of the elements on the fault.
+      Let "h_fault" be the average linear size of the elements on the fault plane (they usually have similar size).
+      Eta must be a small fraction of the critical time step in an elastic medium 
+      for a (hypothetical) element of cubic shape with size equal to h_fault.
+      This (hypothetical) timestep value "dtc_fault" must be computed before simulation
+      using the matlab function utils/critical_timestep.m
+      Note that in general the critical timestep "dtc_bulk" of the whole simulation is smaller, 
+      because elements off the fault might be smaller or more distorted than element faces on the fault.
+      Set eta in Par_file_faults equal to (0.1 to 0.3)*critical_timestep(cp,h,ngll) .
+      A larger eta damps high-frequencies more aggresively but also affects lower frequencies
+      and affects rupture speed.
 
+      Warning: viscosity degrades numerical stability:
+      The critical timestep in a viscous simulation needs to be smaller than what
+      you would usually expect in a purely elastic simulation.
+      The critical timestep for a Kelvin-Voigt material is
+        dtc_kv = eta*( sqrt(1+dtc_bulk^2/eta^2)-1 )
+      where dtc_bulk is the critical timestep for a purely elastic medium.
+      To get the value of dtc_bulk, first run a test simulation with eta=0,
+      and look for the "maximum suggested time step" in OUTPUT_FILES/output_mesher.txt. 
+
+
 DATA/FAULT/FAULT_STATIONS	Stations in the fault plane.
 	Line 1: number of stations.
 	Line 2 to end: 5 columns: X, Y, Z (-depth), station name, fault-id

Added: seismo/3D/FAULT_SOURCE/branches/new_fault_db/utils/critical_timestep.m
===================================================================
--- seismo/3D/FAULT_SOURCE/branches/new_fault_db/utils/critical_timestep.m	                        (rev 0)
+++ seismo/3D/FAULT_SOURCE/branches/new_fault_db/utils/critical_timestep.m	2012-01-14 01:49:16 UTC (rev 19359)
@@ -0,0 +1,31 @@
+% CRITICAL_TIMESTEP	Computes the critical timestep in 3D, assuming:
+%                       purely elastic medium
+%                       the leapfrog time scheme
+%                       a cube element (same size for all faces)
+%
+% SYNTAX	dtc = critical_timestep(csp,h,ngll)
+%
+% INPUT		cp	P wave speed (in m/s)
+%		h	element size (in m)
+%		ngll	number of GLL nodes per spectral element edge (2 to 20)
+%
+% OUTPUT	dtc	critical timestep
+%
+function dtc = critical_timestep(csp,h,ngll)
+
+DIM = 3; % dimension
+
+if ngll<2 | ngll>20, error('ngll must be from 2 to 20'); end
+
+% tabulated values of critical frequency (non-dimensional, 1D)
+% Omega_max(p) = Omega_max(ngll-1)
+Omega_max = [2.0000000e+00 4.8989795e+00 8.6203822e+00 1.3540623e+01 1.9797952e+01 2.7378050e+01 ...
+     3.6256848e+01 4.6421894e+01 5.7867306e+01 7.0590158e+01 8.4588883e+01 9.9862585e+01 ...
+     1.1641072e+02 1.3423295e+02 1.5332903e+02 1.7369883e+02 1.9534221e+02 2.1825912e+02 2.4244948e+02];
+
+% stability factor por leapfrog timescheme
+C = 2;
+
+% critical time step, 
+% assumes a cube element
+dtc = C*h/csp/sqrt(DIM)/Omega_max(ngll-1);