[cig-commits] r22978 - in seismo/3D/SPECFEM3D/trunk: . utils

[ampuero at geodynamics.org]

Author: ampuero
Date: 2013-11-08 17:38:10 -0800 (Fri, 08 Nov 2013)
New Revision: 22978

Added:
   seismo/3D/SPECFEM3D/trunk/utils/FSEM3D_moment_area.m
Modified:
   seismo/3D/SPECFEM3D/trunk/todo_list_please_dont_remove.txt
   seismo/3D/SPECFEM3D/trunk/utils/FSEM3D_snapshot.m
Log:
matlab script to compute seismic potency and rupture area from dynamic rupture outputs

Modified: seismo/3D/SPECFEM3D/trunk/todo_list_please_dont_remove.txt
===================================================================
--- seismo/3D/SPECFEM3D/trunk/todo_list_please_dont_remove.txt	2013-11-04 08:16:42 UTC (rev 22977)
+++ seismo/3D/SPECFEM3D/trunk/todo_list_please_dont_remove.txt	2013-11-09 01:38:10 UTC (rev 22978)
@@ -462,9 +462,11 @@
 
 + fault_solver_dynamic:
       - many hard-coded ad hoc features need to be set instead as user options
+        (e.g. rate-and-state friction)
 
 + rate-and-state friction:
       - make it a user-friendly option (currently a flag in fault_solver)
+      - add documentation
 
 + cubit interface:
       - consolidate python scripts (eliminate redundancy)

Added: seismo/3D/SPECFEM3D/trunk/utils/FSEM3D_moment_area.m
===================================================================
--- seismo/3D/SPECFEM3D/trunk/utils/FSEM3D_moment_area.m	                        (rev 0)
+++ seismo/3D/SPECFEM3D/trunk/utils/FSEM3D_moment_area.m	2013-11-09 01:38:10 UTC (rev 22978)
@@ -0,0 +1,68 @@
+%FSEM3D_moment_area computes seismic moment and rupture area
+%
+% [M,A] = FSEM3D_moment_area(ngll,isnap, [data_dir,db_dir,Dmin])
+%
+% INPUTS	ngll 	number of GLL points per element edge (parameter NGLL of SPECFEM3D)
+%		isnap	last snapshot index in Snapshot*.bin file names (contains final slip)
+%		dat_dir	["OUTPUT_FILES"] directory containing the files Snapshot*.bin
+%		db_dir	["OUTPUT_FILES/DATABASES_MPI"] directory containing the files proc*fault_db.bin
+%		Dmin	minimum slip to define the rupture area
+%
+% OUTPUTS	Px	seismic potency along-strike (integral of along-strike slip)
+% 		Pz	seismic potency along-dip (integral of along-dip slip)
+%		A	area of regions with slip > Dmin
+%
+% WARNING: implemented only to process the first fault
+% 
+% Jean-Paul Ampuero	ampuero at gps.caltech.edu
+
+%		fault	[1] fault id
+
+function [Px,Pz,A] = FSEM3D_moment_area(ngll,isnap,data_dir,db_dir,Dmin) %,fault)
+
+if nargin<3, data_dir = 'OUTPUT_FILES'; end
+if nargin<4, db_dir = 'OUTPUT_FILES/DATABASES_MPI'; end
+if nargin<5, Dmin=1e-3; end 
+
+%if nargin<6, fault = 1; end
+fault=1;
+
+% read slip
+dat = FSEM3D_snapshot(isnap,data_dir,fault);
+Dx=dat.Dx;
+Dz=dat.Dz;
+clear dat
+
+list = dir([db_dir '/*fault_db.bin']); 
+nproc = length(list);
+
+Px=0;
+Pz=0;
+A=0;
+i0=0;
+
+for p=1:nproc
+
+  fid=fopen([db_dir '/' list(p).name]);
+  dat = fread(fid,3,'int');
+ % WARNING: reading the first fault only
+  dat = fread(fid,4,'int');
+  nspec=dat(2);
+  nglob=dat(3);
+  if (nspec==0), fclose(fid); continue; end % if this proc does not contain the fault, skip to next proc
+  dat = fread(fid,nspec*ngll^2+2,'int');
+  ibool = dat(2:end-1);
+  dat = fread(fid,nspec*ngll^2+2,'float');
+  jacw = dat(2:end-1)';
+  fclose(fid);
+
+  Dpx = Dx(i0+ibool);
+  Dpz = Dz(i0+ibool);
+  Px = Px + jacw*Dpx;
+  Pz = Pz + jacw*Dpz;
+  Dp2 = Dpx.*Dpx+Dpz.*Dpz;
+  A = A + jacw*(Dp2>Dmin^2);
+
+  i0 = i0+nglob;
+
+end

Modified: seismo/3D/SPECFEM3D/trunk/utils/FSEM3D_snapshot.m
===================================================================
--- seismo/3D/SPECFEM3D/trunk/utils/FSEM3D_snapshot.m	2013-11-04 08:16:42 UTC (rev 22977)
+++ seismo/3D/SPECFEM3D/trunk/utils/FSEM3D_snapshot.m	2013-11-09 01:38:10 UTC (rev 22978)
@@ -16,8 +16,8 @@
 %			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.
+% Jean-Paul Ampuero	ampuero at erdw.ethz.ch 
+% 19/01/2011: modified by Percy Galvez percy.galvez at sed.ethz.ch 
 %
 % WARNING : Works only for single precision snapshot files.