RSQSim Restart and Processing Instructions

Restarting an RSQSim simulation in order to extend a catalog from where it ended

Note: You'll need some of the R scripts from RSQSimPostProcess

1. Create a new directory for the simulation you are setting up, and copy over the input file, pbs script, fault file (e.g. *.flt, or zfault_Deepen.in), *.KZero, and *.neighbors files. For this example, we'll call the previous catalog "cat1" and the extension "cat2"

2. Rename cat1.in and cat1.pbs to match the new catalog name (cat2.in and cat2.pbs).

3. Update the catalog name and input file name in the new cat2.pbs, and the outNameInfix in the new cat2.in files.

4. Read in catalog 1:

eqs = readEqs("/path/to/eqs.cat1.out")

5. Create the new input files for shear stress, normal stress, and theta from the end of catalog 1 and fill those in for the variables initTauFname, initSigmaFname, and initThetaFname in your cat2.in file:

tmp = mkInitTauSigmaThetaSlipSpeed(eqs, Inf, writeTau = "final", writeSlipSpeed = 0,   
		                   initTauFile="cat2.initTau",
		                   initSigmaFile="cat2.initSigma",
	                           initThetaFile = "cat2.initTheta")

• If the final snapshots (.out.final files) weren’t written out for cat1, but snapshots were written out during the simulations, you can change writeTau = “final” to writeTau = 2 in order to use those snapshots (out.2 files). If no snapshots were written then you’re out of luck and can’t restart that simulation.

6. Get the start time for cat2 and copy the whole number (all 22 digits) into the variable tStart in the new cat2.in file:

tStart = format(tmp$t, digits=22) 	 

7. Read or create a list of pinned patches using the variable pin:

If there is a pin file for cat1:

pin = scan("../cat1.pin") 

Otherwise, the pinned file is just a list of 0's (not pinned) or 1's (pinned) for each patch in the fault model, so if there isn't one, you can just create one of 0's:

pin = rep(0, eqs$fault$np)

8. Update pin from the list of patches that got locked in cat1 and write out the new pin file for cat2:

system("grep Eliminating ../*e | awk '{print $4}' > ../cat1.locked")
pin2 = scan(file="../cat1.locked")			
pin[pin2] = 1						
write(pin, file="cat2.pin", ncol=1)				

9. Update the pinnedFname in your cat2.in file to cat2.pin

10. Update maxT in your cat2.in file to the length you want the extended catalog to be (in seconds).

Combining extended/restarted catalogs into one long catalog

Note: You'll need some of the R scripts from RSQSimPostProcess

1. Set up a list to tell R to only load the files you need (eList, pList, dList, and tList):

L = c("e", "p", "d", "t")

2. Set up the list of catalogs to combine (in the order that they ran):

eqFiles = c("home/user/name/catalog/cat1/eqs.cat1.out",
	    "home/user/name/catalog/cat2/eqs.cat2.out")

3. Read and combine them into one catalog:

eqs = readAndCombineEqfiles(eqFiles, returnLists=L)

4. Write out the new eqs.out, and List files:

writeEqsAndLists(eqs, outFnameInfix = paste("combinedCat", sep = ""))

Note: Reading in an eqs.RData file is faster so it’s helpful save that version too. The whole combined ‘eqs’ data structure will be stored in the combinedCat.RData file:

save(eqs, file = "combinedCat.RData")

Final Note: Use the R function load() to read the combinedCat.RData file back in (only use readEqs.R on the eqs.out files).

load("combinedCat.RData")

5. If you also saved and need to combine your transitions file:

transFiles = c(“home/user/name/catalog/cat1/trans.cat1.out",
	       "home/user/name/catalog/cat2/trans.cat2.out")

combineTransFiles(transfiles=transFiles, outtransfile = "trans.combinedCat.out")

Filtering a catalog by magnitude and saving a new filtered catalog

1. Read the full catalog in (or load the eqs.RData file):

eqs = readEqs("eqs.cat.out")			

2. Pick a minimum magnitude M for your new catalog (you'll get M=muse & greater):

muse = 7 

3. Get a list of the events >= muse

euse = which(eqs$M>muse)

4. Subset/filter the catalog by euse (set renumberEvents=FALSE if you want to keep the original event ID's):

eqsNew = subsetEqs(eqs,euse, renumberEvents=FALSE)

5. Save the new filtered catalog:

subName = "new_catalog_M7"				
writeEqsAndLists(eqsNew, outFnameInfix = paste(subName, sep = ""))

Filtering a catalog by event time and removing events in the simulation runup

Note: To save the filtered catalog, see above instructions.

1. Read the full catalog in (if you have already done this in your current R session, you don’t need to read it in again):

eqs = readEqs("eqs.cat.out")		

2. Check the range of event times in the catalog (in catalog years):

range(eqs$t0yr)

3. Choose the minimun and maximum catalog times you’re interested in:

Note: The simulation runup time depends on the slip rates of your faults (you want each one to rupture at least once). If you are running a simulation with the UCERF3 California fault model and you’re interested in the San Andreas Fault, a reasonable runup time is 5,000 years.

tmin = 5000
tmax = max(eqs$t0yr)

4. Get a list of all events in your specified range (in this example, after the runup time):

inRange = which(eqs$t0yr > tmin & eqs$t0yr <tmax)	

Note: You can combine the time and magnitude filtering into a single which function:

Mmin = 7
Mmax = 8

inRange = which(eqs$t0yr > tmin & eqs$t0yr <tmax & eqs$M >= Mmin & eqs$M <=Mmax)

Filtering a catalog by which fault(s) ruptured

1. Read the full catalog in (if you have already don’t this in your current R session, you don’t need to read it in again):

eqs = readEqs("eqs.cat.out")		

2. Make a list of the fault segments you are interested in:

Note: If you are using the UCERF3 fault model, you can pass eqs$fault$faultName to the function unique() to get a list of the options.

faults = which(eqs$fault$faultName %in% c("SanAndreas(MojaveN)", "SanAndreas(MojaveS)"))

See the next instruction set (Plotting the fault model in plotFault3d) to visually check that you’ve selected the correct fault segments.

3. Get a list of all of the events that rupture patches of those fault segments:

euse = unique(eqs$eList[which(eqs$pList %in% faults)])	

Note: You can easily filter the catalog further by the range of events in the previous instruction set:

use = euse[which(euse %in% inRange)]	

Plotting the fault model in plotFault3d

1. Go through steps 1 and 2 in the previous instruction set to get the list of fault patches you want to work with.

2. Initialize a variable for the plot colors:

Note: This will create list of NA (which will show up as clear when you plot it) that is the length of the number of fault patches in the model.

col = rep(NA, length(eqs$fault$np))

3. Set the color of the fault patches that match segments you selected:

col[faults]= "red"

4. Plot the fault model, colored by the color variable you created:

plotFault3d(eqs$fault, col = col, lwd = 0.5)

Note: The argument lwd sets the line width. To see the additional arguments you can pass to plotFault3d() use the args() function:

args(plotFault3d)