[cig-commits] r22132 - seismo/3D/SPECFEM3D/trunk/utils
surendra at geodynamics.org
surendra at geodynamics.org
Thu May 23 12:19:41 PDT 2013
Author: surendra
Date: 2013-05-23 12:19:41 -0700 (Thu, 23 May 2013)
New Revision: 22132
Added:
seismo/3D/SPECFEM3D/trunk/utils/FSEM3D_snapshot.m
seismo/3D/SPECFEM3D/trunk/utils/critical_timestep.m
Log:
Added post- and pre- processing utilities for faults with split nodes
Added: seismo/3D/SPECFEM3D/trunk/utils/FSEM3D_snapshot.m
===================================================================
--- seismo/3D/SPECFEM3D/trunk/utils/FSEM3D_snapshot.m (rev 0)
+++ seismo/3D/SPECFEM3D/trunk/utils/FSEM3D_snapshot.m 2013-05-23 19:19:41 UTC (rev 22132)
@@ -0,0 +1,61 @@
+%FSEM3D_SNAPSHOT reads fault data at a given time
+%
+% d = FSEM3D_snapshot(isnap, [dir, fault])
+%
+% INPUTS isnap snapshot index, as in Snapshot*.bin file names
+% dir ["."] directory containing the SPECFEM3D output data Snapshot*.bin
+% fault [1] fault id
+%
+% OUTPUTS d structure containing fault fields:
+% X,Y,Z fault node coordinates (km)
+% Dx,Dz slip (m)
+% Vx,Vz slip rate (m/s)
+% Tx,Ty,Tz stress change (MPa)
+% S slip "state" variable in friction law (m)
+% Sg strength relative to initial stress (MPa)
+% Trup rupture time (s)
+% Tpz process zone time, when slip=Dc (s)
+%
+% Jean-Paul Ampuero ampuero at erdw.ethz.ch modified by
+% Percy Galvez percy.galvez at sed.ethz.ch 19/01/2011.
+%
+% WARNING : Works only for single precision snapshot files.
+
+function d = FSEM3D_snapshot(isnap,DATA_DIR,fault)
+
+NDAT = 14;
+
+if nargin<2, DATA_DIR = '.'; end
+if nargin<3, fault = 1; end
+
+BinFile = sprintf('%s/Snapshot%u_F%u.bin',DATA_DIR,isnap,fault);
+
+if ~exist(BinFile,'file'), error(sprintf('File %s does not exist',BinFile)), end
+fid=fopen(BinFile);
+BinRead = fread(fid,[1,inf],'single')' ;
+fclose(fid);
+
+BinRead = reshape( BinRead(:),length(BinRead)/NDAT, NDAT);
+BinRead = BinRead(2:end-1,:);
+
+% Reorder fault nodes (lexicographic z,x)
+%[LOC,IND] = sortrows( BinRead(:,[1 3]) );
+%BinRead = BinRead(IND,:);
+
+d.X = BinRead(:,1)/1e3; % in km
+d.Y = BinRead(:,2)/1e3; % in km
+d.Z = BinRead(:,3)/1e3; % in km
+d.Dx = BinRead(:,4);
+d.Dz = BinRead(:,5);
+d.Vx = BinRead(:,6);
+d.Vz = BinRead(:,7);
+d.Tx = BinRead(:,8); % in MPa
+d.Ty = BinRead(:,9);
+d.Tz = BinRead(:,10); % in MPa
+d.S = BinRead(:,11);
+d.Sg = BinRead(:,12); % in MPa
+d.Trup = BinRead(:,13);
+d.Tpz = BinRead(:,14);
+
+
+return
Added: seismo/3D/SPECFEM3D/trunk/utils/critical_timestep.m
===================================================================
--- seismo/3D/SPECFEM3D/trunk/utils/critical_timestep.m (rev 0)
+++ seismo/3D/SPECFEM3D/trunk/utils/critical_timestep.m 2013-05-23 19:19:41 UTC (rev 22132)
@@ -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);
More information about the CIG-COMMITS
mailing list