Merge pull request #188 from johnmaxwilson/master

Improved meshing, updated manual for software citation

AUTHORS

@@ -1,5 +1,5 @@
+John M. Wilson		<jhnwilson@ucdavis.edu>
+Kasey W. Schultz    <kwschultz@ucdavis.edu>
 Eric M. Heien		<emheien@ucdavis.edu>
 Michael K. Sachs	<mksachs@ucdavis.edu>
-Kasey W. Schultz    <kwschultz@ucdavis.edu>
 John B. Rundle		<jbrundle@ucdavis.edu>
-John M. Wilson		<jhnwilson@ucdavis.edu>

CMakeLists.txt

@@ -22,7 +22,7 @@ PROJECT (VQ)
 LIST(APPEND CMAKE_MODULE_PATH ${CMAKE_CURRENT_SOURCE_DIR}/quakelib/cmake/)
 
 # Set the version number
-SET(VQ_VERSION_MAJOR 2)
+SET(VQ_VERSION_MAJOR 3)
 SET(VQ_VERSION_MINOR 1)
 SET(VQ_VERSION_SUBMINOR 0)
 SET(VQ_VERSION_STR "${VQ_VERSION_MAJOR}.${VQ_VERSION_MINOR}.${VQ_VERSION_SUBMINOR}")

PyVQ/pyvq/pyvq.py

@@ -54,7 +54,13 @@
 try:
     import cPickle as pickle
 except ImportError:
-    cPickle_available = False
+    cPickle_available = False    
+
+netCDF4_available = True
+try:
+    from netCDF4 import Dataset
+except ImportError:
+    netCDF4_available = False
 
 
 # ----------------- Global constants -------------------------------------------
@@ -2221,21 +2227,57 @@ def __init__(self, geometry, event_id, event, element_slips, LLD_file):
         for i in range(len(self.lons_1d)):
             self.grid_1d.append(self.convert.convert2xyz(quakelib.LatLonDepth(self.lats_1d[i],self.lons_1d[i])))
     #            
-    def compute_field(self):        
+    def compute_field(self, netCDF_arg):        
         self.lame_lambda = 3.2e10
         self.lame_mu     = 3.0e10
         self.field_1d    = self.slip_map.displacements(self.grid_1d, self.lame_lambda, self.lame_mu, 1e9)
-        outname = self.LLD_file.split(".tx")[0]+"_dispField_event"+str(self.event_id)+".txt"
-        outfile = open(outname,'w')
-        # Write the header with the number of points
-        outfile.write("#### number of points ####\n")
-        outfile.write("{}\n".format(len(self.field_1d)))
-        outfile.write("##########################\n")
-        for i in range(len(self.field_1d)):
-            outfile.write("{}\t{}\t{}\n".format(self.lats_1d[i], self.lons_1d[i], self.field_1d[i][2]))
-        outfile.close()
-        sys.stdout.write("\n---> Event displacements written to "+outname)
-        sys.stdout.write("\n")
+        outname = self.LLD_file.split(".tx")[0]+"_dispField_event"+str(self.event_id)
+
+        #move file extension inside these
+        if netCDF_arg:
+            if not netCDF4_available:
+                raise BaseException("\nCannot save as netCDF, netCDF4 module is not available.")
+            else:
+                #Heres the netCDF file builder
+                outname = outname+".nc"
+                lats = []
+                lons = []
+                uplift = []
+                for i, lon in enumerate(self.lons_1d):
+                    if self.lats_1d[i] not in lats: lats.append(self.lats_1d[i])
+                    if lon not in lons: lons.append(lon)
+                    uplift.append(self.field_1d[i][2])
+
+                uplift = np.array(uplift)
+                uplift = np.reshape(uplift, (np.size(lats), np.size(lons)))
+
+                uplift_dataset = Dataset(outname, 'w', format='NETCDF4')
+
+                uplift_dataset.createDimension('lat', len(lats))
+                uplift_dataset.createDimension('lon', len(lons))
+
+                lats_data   = uplift_dataset.createVariable('latitude', 'f4', ('lat',))
+                lons_data   = uplift_dataset.createVariable('longitude', 'f4', ('lon',))
+                uplift_data = uplift_dataset.createVariable('uplift', 'f4', ('lat','lon'))
+
+                lats_data[:]     = lats
+                lons_data[:]     = lons
+                uplift_data[:,:] = uplift
+
+                uplift_dataset.close()
+
+        else:
+            outname = outname+".txt"
+            outfile = open(outname,'w')
+            # Write the header with the number of points
+            outfile.write("#### number of points ####\n")
+            outfile.write("{}\n".format(len(self.field_1d)))
+            outfile.write("##########################\n")
+            for i in range(len(self.field_1d)):
+                outfile.write("{}\t{}\t{}\n".format(self.lats_1d[i], self.lons_1d[i], self.field_1d[i][2]))
+            outfile.close()
+            sys.stdout.write("\n---> Event displacements written to "+outname)
+            sys.stdout.write("\n")
 
 
 
@@ -3115,6 +3157,7 @@ def leonard_2010(self, type, min_mag, max_mag, num=5):
             help="Levels for contour plot.")
     parser.add_argument('--small_model', required=False, action='store_true', help="Small fault model, used to specify map extent.")    
     parser.add_argument('--field_eval', required=False, action='store_true', help="Evaluate an event field at specified lat/lon. Must provide the file, --lld_file")
+    parser.add_argument('--netCDF', required=False, action='store_true', help="Store field evaluation in netCDF format")
     parser.add_argument('--lld_file', required=False, help="File containing lat/lon columns to evaluate an event field.")
 
     # ---------  Greens function plotting arguments -----------
@@ -3601,7 +3644,7 @@ def leonard_2010(self, type, min_mag, max_mag, num=5):
             sys.stdout.write(" Loaded slips for {} elements :".format(len(ele_slips.keys()))) 
         sys.stdout.flush()
         FE = FieldEvaluator(geometry, args.event_id, event, ele_slips, args.lld_file)
-        FE.compute_field()
+        FE.compute_field(args.netCDF)
     if args.greens:
         # Set default values
         if args.dip is None: sys.exit("Must specify --dip")

doc/vq.tex

@@ -41,8 +41,8 @@
 
 %%%%%%%%%%%%%%%%%%%%%%%%%%% Lyx Title Page
 \title{Virtual Quake User Manual}
-\author{Kasey W. Schultz\\Eric M. Heien\\Michael K. Sachs\\Mark R. Yoder\\John M. Wilson\\John B. Rundle\\Donald L. Turcotte\\\\ \copyright University of California, Davis\\
-Version 2.2.0}
+\author{John M. Wilson\\Kasey W. Schultz\\Eric M. Heien\\Michael K. Sachs\\Mark R. Yoder\\John B. Rundle\\Donald L. Turcotte\\\\ \copyright University of California, Davis\\
+Version 3.1.0}
 %\date{\noindent \today}
 
 
@@ -107,7 +107,7 @@
 \hfill{\Huge \fontfamily{\sfdefault}\selectfont User Manual \\
 % FILL: manual version
 % e.g. 1.0
-\raggedleft \huge \fontfamily{\sfdefault}\selectfont Version {3.0.0}\\}
+\raggedleft \huge \fontfamily{\sfdefault}\selectfont Version {3.1.0}\\}
 
 %AUTHOR(S) & WEBSITE%
 \null
@@ -117,7 +117,7 @@
 % FILL: author list
 % e.g. Author One\\Author Two\\Author Three\\
 % be sure to have a newline (\\) after the final author
-Kasey W. Schultz~~~Eric M. Heien~~~Michael K. Sachs~~~John M. Wilson\\Mark R. Yoder~~~John B. Rundle~~~Donald L. Turcotte\\
+John M. Wilson~~~Kasey W. Schultz~~~Eric M. Heien~~~Michael K. Sachs\\Mark R. Yoder~~~John B. Rundle~~~Donald L. Turcotte\\
 }
 {\fontfamily{\sfdefault}\selectfont www.geodynamics.org}
 
@@ -193,11 +193,16 @@ \section*{Citation}
 practice by citing the appropriate peer reviewed papers and making
 appropriate acknowledgements. 
 
-The Virtual Quake development team asks that you cite 
+To cite this software:
+\begin{itemize}
+\item Wilson J.M., Schultz K.W., Heien E.M., Sachs M.K., Yoder M.R., Rundle J.B., Turcotte D.L. (2017), Virtual Quake v3.1.0 [software], Computational Infrastructure for Geodynamics, Available from: geodynamics.org, doi: tempVQ7, url: \url{https://geodynamics.org/cig/software/vq/}
+\end{itemize}
+
+Additionally, the Virtual Quake development team asks that you cite 
 the following:
 \begin{itemize}
+\item Schultz, K. W. and Sachs, M. K. and Heien, E. M. and Rundle, J. B. and Turcotte, D. L. and Donnellan, A. (2016) "Simulating Gravity Changes in Topologically Realistic Driven Earthquake Fault Systems: First Results", Pure Appl. Geophys. \href{http://dx.doi.org/10.1007/s00024-014-0926-4}{doi:10.1007/s00024-014-0926-4}
 \item Schultz, K. W. and Sachs, M. K. and Heien, E. M. and Yoder, M. R. and Rundle, J. B. and Turcotte, D. L. and Donnellan, A. (2015) "Virtual Quake: Statistics, Co-Seismic Deformations and Gravity Changes for Driven Earthquake Fault Systems,", International Association of Geodesy Symposia \href{http://dx.doi.org/10.1007/1345_2015_134}{doi:10.1007/1345\_2015\_134}
-\item Schultz, K. W. and Sachs, M. K. and Heien, E. M. and Rundle, J. B. and Turcotte, D. L. and Donnellan, A. (2014) "Simulating Gravity Changes in Topologically Realistic Driven Earthquake Fault Systems: First Results", Pure Appl. Geophys. \href{http://dx.doi.org/10.1007/s00024-014-0926-4}{doi:10.1007/s00024-014-0926-4}
 \item Eric M. Heien, Michael Sachs, "Understanding Long-Term Earthquake Behavior through Simulation," Computing in Science and Engineering, vol. 14, no. 5, pp. 10-20, Sept.-Oct. 2012, \href{http://dx.doi.org/10.1109/MCSE.2012.39}{doi:10.1109/MCSE.2012.39}
 \item Sachs, M.K., Heien, E.M., Turcotte, D.L., Yikilmaz, M.B., Rundle,
 J.B., Kellogg, L.H. ''Virtual California Earthquake Simulator''
@@ -207,7 +212,7 @@ \section*{Citation}
 
 And to cite the manual:
 \begin{itemize}
-\item K. W. Schultz, E. M. Heien, M. K. Sachs, M. R. Yoder, J. M. Wilson,  J. B. Rundle, D. L. Turcotte (2016), Virtual Quake User Manual, Version 2.2.0. Davis, California, USA: Computational Infrastructure for Geodynamics. URL:  \href{http://geodynamics.org/cig/software/github/vq/v2.2.0/vq_manual_2.2.0.pdf}{geodynamics.org/cig/software/github/vq/v2.2.0/vq\_manual\_2.2.0.pdf}
+\item Wilson J.M., Schultz K.W., Heien E.M., Sachs M.K., Yoder M.R., Rundle J.B., Turcotte D.L. (2017), Virtual Quake User Manual, Version 3.1.0. Davis, California, USA: Computational Infrastructure for Geodynamics. URL:  \href{https://geodynamics.org/cig/software/vq/vq_manual_3.1.0.pdf}{https://geodynamics.org/cig/software/vq/vq\_manual\_3.1.0.pdf}
 \end{itemize}
 
 The developers also request that in your oral presentations and in

helper_scripts/param_generation.py

@@ -4,12 +4,21 @@
 import readline
 readline.parse_and_bind("tab: complete")
 
+print "%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%"
 print "\nThis script will help you generate a Virtual Quake parameter file.\n"
 print "You will be asked about commonly modified parameters.\n"
 print "If you would like to make further customizations, please consult the user manual and the example parameter file to make manual adjustements to your paramter file.\n"
 print "\n"
+print "After running this script, move the \"generated_parameters.d\" file to the location of your VQ run.\n"
+print "A simulation can then be initiated with $[Path to the vq program] ./generated_parameters.d"
+print "%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%"
+print "\n"
+
+SIM_LOCATION = raw_input("What is the full path of the full path of the directory where you'll be running this simulation?\n")
+while not os.path.isdir(SIM_LOCATION):
+        SIM_LOCATION = raw_input("Can't find directory, please input the full path of the directory where you'll be running this simulation\n")
 
-year_input = raw_input("How many years would you like to simulate?\n")
+year_input = raw_input("\nHow many years would you like to simulate?\n")
 while True:
     try:
         END_YEAR = float(year_input)
@@ -19,42 +28,52 @@
     else:
         break
 
-greens_exist = raw_input("Have you already computed the Green's functions? (y/n)\n")
+greens_exist = raw_input("\nHave you already computed the Green's functions? (y/n)\n")
 while greens_exist not in ['y','Y','n','N']:
     greens_exist = raw_input("Please enter \'y\' or \'n\'\n")
 if greens_exist in ['y', 'Y']:
     GREENS_METHOD = "file"
-    GREENS_PATH = raw_input("What is the full path of the file where the Green's functions are stored?\n")
-    while not os.path.isfile(GREENS_INPUT):
-        GREENS_PATH = raw_input("Can't find file, please input the full path of the Green's function file\n")
+    GREENS_IOMETH = "sim.greens.input"
+    GREENS_PATH = raw_input("\nWhat is the path (relative to your simulation directory) of the file where the Green's functions are stored?\n")
+    while not os.path.isfile(SIM_LOCATION+GREENS_PATH):
+        GREENS_PATH = raw_input("Can't find file, please input the path (relative to your simulation directory) of the Green's function file\n")
 else:
     GREENS_METHOD = "standard"
-    GREENS_PATH = raw_input("Please enter the full path of the file where you'd like to save the Green's functions\n")
-    while not os.path.isdir(os.path.dirname(GREENS_PATH)):
-        GREENS_PATH = raw_input("Can't find directory, please input the full path where you'd like to save the Green's functions\n")
-
-MODEL_FILE = raw_input("What is the full path of the fault model file for this simulation?\n")
-while not os.path.isfile(MODEL_FILE):
-        MODEL_FILE = raw_input("Can't find file, please input the full path of the model file\n")
-
+    GREENS_IOMETH = "sim.greens.output"
+    GREENS_PATH = raw_input("\nPlease enter the path (relative to your simulation directory) of the file where you'd like to save the Green's functions\n")
+    while not os.path.isdir(os.path.dirname(SIM_LOCATION+GREENS_PATH)):
+        GREENS_PATH = raw_input("Can't find directory, please input the path (relative to your simulation directory) where you'd like to save the Green's functions\n")
+
+MODEL_FILE = raw_input("\nWhat is the path (relative to your simulation directory) of the fault model file for this simulation?\n")
+while not os.path.isfile(SIM_LOCATION+MODEL_FILE):
+        MODEL_FILE = raw_input("Can't find file, please input the path (relative to your simulation directory) of the model file\n")
+
+if '.txt' in MODEL_FILE:
+    MODEL_TYPE = 'text'
+elif '.h5' in MODEL_FILE:
+    MODEL_TYPE = 'hdf5'
+else:
+    MODEL_TYPE = raw_input("\nWhat is the file type for your fault model file? (hdf5 or text)\n")
+    while MODEL_TYPE not in ['hdf5', 'text']:
+        MODEL_TYPE = raw_input("Please enter \'hdf5\' or \'text\'\n")
 
-EVENTS_FILE = raw_input("What is the full path of the file where you would like to save the events for this simulation?\n")
-while not os.path.isdir(os.path.dirname(EVENTS_FILE)):
-    EVENTS_FILE = raw_input("Can't find directory, please input the full path where you'd like to save the events\n")
+EVENTS_FILE = raw_input("\nWhat is the path (relative to your simulation directory) of the file where you would like to save the events for this simulation?\n")
+while not os.path.isdir(os.path.dirname(SIM_LOCATION+EVENTS_FILE)):
+    EVENTS_FILE = raw_input("Can't find directory, please input the path (relative to your simulation directory) where you'd like to save the events\n")
 
 
 
-parameter_list = ['sim.time.end_year', 'sim.system.progress_period', 'sim.greens.use_normal', 'sim.greens.method', 'sim.greens.input', 'sim.file.input',
+parameter_list = ['sim.time.end_year', 'sim.system.progress_period', 'sim.greens.use_normal', 'sim.greens.method', GREENS_IOMETH, 'sim.file.input',
                   'sim.file.input_type', 'sim.file.output_event', 'sim.file.output_event_type', 'sim.friction.dynamic', 'sim.friction.dynamic_stress_drops',
                   'sim.greens.shear_offdiag_min', 'sim.greens.shear_offdiag_max', 'sim.greens.normal_offdiag_min', 'sim.greens.normal_offdiag_max']
 
-value_list = [END_YEAR, 5, 'true', GREENS_METHOD, GREENS_PATH, MODEL_FILE, 'hdf5', EVENTS_FILE, 'hdf5', 0.2, 'true',
+value_list = [END_YEAR, 5, 'true', GREENS_METHOD, GREENS_PATH, MODEL_FILE, MODEL_TYPE, EVENTS_FILE, 'hdf5', 0.2, 'true',
               -745348, 3590676, -197265.50322689, 202035.55337244]
 
 
 param_file = open("generated_parameters.d", "w")
 
-print "\"Parameters saved as generated_parameters.d.\""
+print "Parameters saved as \"generated_parameters.d\".  Remember to move this file to your simulation location before running"
 
 
 for i, param in enumerate(parameter_list):