diff --git a/PyVQ/pyvq/pyvq.py b/PyVQ/pyvq/pyvq.py index 261fe20a..db905980 100755 --- a/PyVQ/pyvq/pyvq.py +++ b/PyVQ/pyvq/pyvq.py @@ -99,20 +99,26 @@ def calculate_averages(x,y,log_bin=False,num_bins=None): return x_ave, y_ave class SaveFile: - def event_plot(self, event_file, plot_type, min_mag): + def event_plot(self, event_file, plot_type, min_mag, min_year, max_year): min_mag = str(min_mag) # Remove any folders in front of model_file name if len(event_file.split("/")) > 1: event_file = event_file.split("/")[-1] + # Add tags to convey the subsets/cuts being made + add="" + if min_year is not None: add+="_yearMin"+str(int(min_year)) + if max_year is not None: add+="_yearMax"+str(int(max_year)) + if args.use_sections is not None: + for sec in args.use_sections: + add+="_"+geometry.model.section(sec).name() if min_mag is not None: # e.g. min_mag = 7.5, filename has '7-5' if len(min_mag.split(".")) > 1: - mag = "minMag_"+min_mag.split(".")[0]+"-"+min_mag.split(".")[1]+"_" + add += "_minMag_"+min_mag.split(".")[0]+"-"+min_mag.split(".")[1] else: - mag = "minMag_"+min_mag+"_" - else: - mag = "" - return plot_type+"_"+mag+event_file.split(".")[0]+".png" + add += "_minMag_"+min_mag + + return plot_type+add+"_"+event_file.split(".")[0]+".png" def field_plot(self, model_file, field_type, uniform_slip, event_id): # Remove any folders in front of model_file name @@ -134,11 +140,31 @@ def trace_plot(self, model_file): model_file = model_file.split("/")[-1] return "traces_"+model_file.split(".")[0]+".png" - def diagnostic_plot(self, event_file, plot_type): + def diagnostic_plot(self, event_file, plot_type, min_year=None, max_year=None, min_mag=None): # Remove any folders in front of model_file name if len(event_file.split("/")) > 1: event_file = event_file.split("/")[-1] - return plot_type+"_diagnostic_"+event_file.split(".")[0]+".png" + add="" + if min_year is not None: add+="_yearMin"+str(int(min_year)) + if max_year is not None: add+="_yearMax"+str(int(max_year)) + if args.use_sections is not None: + for sec in args.use_sections: + add+="_"+geometry.model.section(sec).name() + if min_mag is not None: + min_mag = str(min_mag) + # e.g. min_mag = 7.5, filename has '7-5' + if len(min_mag.split(".")) > 1: + add += "_minMag_"+min_mag.split(".")[0]+"-"+min_mag.split(".")[1] + else: + add += "_minMag_"+min_mag + + return plot_type+"_diagnostic"+add+"_"+event_file.split(".")[0]+".png" + + def event_kml_plot(self, event_file, event_id): + if len(event_file.split("/")) > 1: + event_file = event_file.split("/")[-1] + event_file = event_file.split("events")[-1] + return "event_"+str(event_id)+event_file.split(".")[0]+".kml" class MagFilter: @@ -265,6 +291,49 @@ def get_fault_traces(self): traces_lat_lon[sid] = [(lat,lon)] #break return traces_lat_lon + + def get_slip_rates(self, elements): + # Convert slip rates from meters/second to meters/(decimal year) + CONVERSION = 3.15576*pow(10,7) + return {id:self.model.element(id).slip_rate()*CONVERSION for id in elements} + + def get_slip_time_series(self, events, elements=None, min_year=None, max_year=None, DT=None): + # slip_time_series = dictionary indexed by block_id with entries being arrays of absolute slip at each time step + # Get slip rates for the elements + slip_rates = self.get_slip_rates(elements) + #Initialize blocks with 0.0 slip at time t=0.0 + slip_time_series = {id:[0.0] for id in elements} + # Grab the events data + event_years = events.event_years() + event_numbers = events.event_numbers() + #Initialize time steps to evaluate slip + time_values = np.arange(min_year+DT, max_year+DT, DT) + for k in range(len(time_values)): + if k>0: + # current time in simulation + right_now = time_values[k] + # back slip all elements by subtracting the slip_rate*dt + for block_id in slip_time_series.keys(): + last_slip = slip_time_series[block_id][k-1] + this_slip = slip_rates[block_id]*DT + slip_time_series[block_id].append(last_slip-this_slip) + # check if any elements slip as part of simulated event in the window of simulation time + # between (current time - DT, current time), add event slips to the slip at current time + # for elements involved + for j in range(len(event_numbers)): + evid = event_numbers[j] + ev_year = event_years[j] + if right_now-DT < ev_year <= right_now: + event_element_slips = events.get_event_element_slips(evid) + for block_id in event_element_slips.keys(): + try: + slip_time_series[block_id][k] += event_element_slips[block_id] + #sys.stdout.write("element {} slips {} in event {}\n".format(block_id,event_element_slips[block_id],evid)) + #sys.stdout.flush() + except KeyError: + pass # Ignore event elements that we are not asked for (in elements) + return slip_time_series + # ======= h5py I/O ============================================ def read_events_h5(sim_file, event_numbers=None): @@ -355,8 +424,11 @@ def event_rupture_areas(self): def event_magnitudes(self): return [self._events[evnum].getMagnitude() for evnum in self._filtered_events if not np.isnan(self._events[evnum].getMagnitude())] -# TODO: Handle NaN magnitudes on the C++ side + # TODO: Handle NaN magnitudes on the C++ side + def event_numbers(self): + return [evnum for evnum in self._filtered_events if not np.isnan(self._events[evnum].getMagnitude())] + def event_mean_slip(self): return [self._events[evnum].calcMeanSlip() for evnum in self._filtered_events if not np.isnan(self._events[evnum].getMagnitude())] @@ -1604,16 +1676,23 @@ def create_plot(self, plot_type, log_y, x_data, y_data, plot_title, x_label, y_l plt.savefig(filename,dpi=100) sys.stdout.write("Plot saved: {}\n".format(filename)) - def multi_line_plot(self, x_data, y_data, colors, labels, linewidths, plot_title, x_label, y_label, legend_str, filename): + def multi_line_plot(self, x_data, y_data, labels, linewidths, plot_title, x_label, y_label, legend_str, filename, colors=None, linestyles=None): fig = plt.figure() ax = fig.add_subplot(111) ax.set_xlabel(x_label) ax.set_ylabel(y_label) - ax.set_title(plot_title) - if not (len(x_data) == len(y_data) and len(x_data) == len(colors) and len(colors) == len(labels) and len(linewidths) == len(colors)): - raise "These lists must be the same length: x_data, y_data, colors, labels, linewidths." - for i in range(len(x_data)): - ax.plot(x_data[i], y_data[i], color=colors[i], label=labels[i], linewidth=linewidths[i]) + if linestyles is None: linestyles = ["-" for each in x_data] + fig.suptitle(plot_title, fontsize=10) + if colors is not None: + if not (len(x_data) == len(y_data) and len(x_data) == len(colors) and len(colors) == len(labels) and len(linewidths) == len(colors)): + raise "These lists must be the same length: x_data, y_data, colors, labels, linewidths." + for i in range(len(x_data)): + ax.plot(x_data[i], y_data[i], color=colors[i], label=labels[i], linewidth=linewidths[i], ls=linestyles[i]) + else: + if not (len(x_data) == len(y_data) and len(x_data) == len(labels) and len(linewidths) == len(y_data)): + raise "These lists must be the same length: x_data, y_data, labels, linewidths." + for i in range(len(x_data)): + ax.plot(x_data[i], y_data[i], label=labels[i], linewidth=linewidths[i], ls=linestyles[i]) ax.legend(title=legend_str, loc='best') plt.gca().get_xaxis().get_major_formatter().set_useOffset(False) plt.savefig(filename,dpi=100) @@ -1688,6 +1767,8 @@ def scatter_and_line(self, log_y, x_data, y_data, line_x, line_y, line_label, pl ax.legend(loc = "best") ax.get_xaxis().get_major_formatter().set_useOffset(False) + if args.zoom: plt.ylim(-5,5) + plt.savefig(filename,dpi=100) sys.stdout.write("Plot saved: {}\n".format(filename)) @@ -1905,7 +1986,7 @@ def plot_p_of_t_multi(self, events, filename, beta=None, tau=None, num_t0=4, num y_lab = r'P(t, t$_0$)' x_lab = r't = t$_0$ + $\Delta$t [years]' plot_title = "" - self.multi_line_plot(x_data, y_data, colors, labels, linewidths, plot_title, x_lab, y_lab, legend_string, filename) + self.multi_line_plot(x_data, y_data, labels, linewidths, plot_title, x_lab, y_lab, legend_string, filename, colors=colors) def plot_dt_vs_t0(self, events, filename, years_since=None): # Plot the waiting times corresponding to 25/50/75% conditional probabilities @@ -2113,6 +2194,19 @@ def wells_coppersmith(self, type, min_mag=None, max_mag=None, num=5): help="Plot normal stress changes for events") parser.add_argument('--event_mean_slip', required=False, action='store_true', help="Plot the mean slip for events") + parser.add_argument('--zoom', required=False, action='store_true', + help="Force zoomed bounds on scatter and line plots") + + # Geometry + parser.add_argument('--slip_rates', required=False, action='store_true', + help="Print element id and slip rate for all elements.") + parser.add_argument('--elements', type=int, nargs='+', required=False, + help="List of elements for filtering.") + parser.add_argument('--slip_time_series', required=False, action='store_true', + help="Return the slip time series for all specified --elements.") + parser.add_argument('--dt', required=False, type=float, help="Time step for slip rate plots, unit is decimal years.") + parser.add_argument('--event_kml', required=False, action='store_true', + help="Save a KML (Google Earth) file of the event elements, colored by event slip.") # Validation/testing arguments parser.add_argument('--validate_slip_sum', required=False, @@ -2175,17 +2269,22 @@ def wells_coppersmith(self, type, min_mag=None, max_mag=None, num=5): args.plot_mag_mean_slip = True args.wc94 = True + # Set a default lower magnitude limit of 5 + if args.min_magnitude is None: + args.min_magnitude = 5 # Set up filters event_filters = [] if args.min_magnitude or args.max_magnitude: event_filters.append(MagFilter(min_mag=args.min_magnitude, max_mag=args.max_magnitude)) + + if args.min_year or args.max_year: event_filters.append(YearFilter(min_year=args.min_year, max_year=args.max_year)) if args.min_slip or args.max_slip: # Setting default lower limit on mean slip of 1cm - if args.min_slip is None: args.min_slip = 0.01 + #if args.min_slip is None: args.min_slip = 0.01 event_filters.append(SlipFilter(min_slip=args.min_slip, max_slip=args.max_slip)) if args.min_event_num or args.max_event_num: @@ -2195,6 +2294,10 @@ def wells_coppersmith(self, type, min_mag=None, max_mag=None, num=5): if not args.model_file: raise "Must specify --model_file for --use_sections to work." event_filters.append(SectionFilter(geometry, args.use_sections)) + # Also grab all the elements from this section in case this is being used to grab element ids + if args.elements is None: + args.elements = [elem_num for elem_num in range(geometry.model.num_elements()) if geometry.model.element(elem_num).section_id() in args.use_sections] + if args.use_trigger_sections: if not args.model_file: @@ -2215,37 +2318,37 @@ def wells_coppersmith(self, type, min_mag=None, max_mag=None, num=5): args.event_shear_stress = True args.event_normal_stress = True args.event_mean_slip = True - args.plot_freq_mag = 3 + args.plot_freq_mag = 1 args.plot_mag_mean_slip = True args.plot_mag_rupt_area = True args.wc94 = True if args.plot_freq_mag: - filename = SaveFile().event_plot(args.event_file, "freq_mag", args.min_magnitude) + filename = SaveFile().event_plot(args.event_file, "freq_mag", args.min_magnitude, args.min_year, args.max_year) if args.plot_freq_mag == 1: UCERF2,b1 = False, False if args.plot_freq_mag == 2: UCERF2,b1 = False, True if args.plot_freq_mag == 3: UCERF2,b1 = True, False if args.plot_freq_mag == 4: UCERF2,b1 = True, True FrequencyMagnitudePlot().plot(events, filename, UCERF2=UCERF2, b1=b1) if args.plot_mag_rupt_area: - filename = SaveFile().event_plot(args.event_file, "mag_rupt_area", args.min_magnitude) + filename = SaveFile().event_plot(args.event_file, "mag_rupt_area", args.min_magnitude, args.min_year, args.max_year) MagnitudeRuptureAreaPlot().plot(events, filename, WC94=args.wc94) if args.plot_mag_mean_slip: - filename = SaveFile().event_plot(args.event_file, "mag_mean_slip", args.min_magnitude) + filename = SaveFile().event_plot(args.event_file, "mag_mean_slip", args.min_magnitude, args.min_year, args.max_year) MagnitudeMeanSlipPlot().plot(events, filename, WC94=args.wc94) if args.plot_prob_vs_t: - filename = SaveFile().event_plot(args.event_file, "prob_vs_time", args.min_magnitude) + filename = SaveFile().event_plot(args.event_file, "prob_vs_time", args.min_magnitude, args.min_year, args.max_year) ProbabilityPlot().plot_p_of_t(events, filename) if args.plot_prob_vs_t_fixed_dt: - filename = SaveFile().event_plot(args.event_file, "p_vs_t_fixed_dt", args.min_magnitude) + filename = SaveFile().event_plot(args.event_file, "p_vs_t_fixed_dt", args.min_magnitude, args.min_year, args.max_year) ProbabilityPlot().plot_conditional_fixed_dt(events, filename) if args.plot_cond_prob_vs_t: - filename = SaveFile().event_plot(args.event_file, "cond_prob_vs_t", args.min_magnitude) + filename = SaveFile().event_plot(args.event_file, "cond_prob_vs_t", args.min_magnitude, args.min_year, args.max_year) if args.beta: ProbabilityPlot().plot_p_of_t_multi(events, filename, beta=args.beta, tau=args.tau) else: ProbabilityPlot().plot_p_of_t_multi(events, filename) if args.plot_waiting_times: - filename = SaveFile().event_plot(args.event_file, "waiting_times", args.min_magnitude) + filename = SaveFile().event_plot(args.event_file, "waiting_times", args.min_magnitude, args.min_year, args.max_year) ProbabilityPlot().plot_dt_vs_t0(events, filename) if args.field_plot: type = args.field_type.lower() @@ -2320,22 +2423,64 @@ def wells_coppersmith(self, type, min_mag=None, max_mag=None, num=5): if args.small_model is None: args.small_model = False TP = TracePlotter(geometry, filename, use_sections=args.use_sections, small_model=args.small_model) + if args.slip_rates: + if args.elements is None: args.elements = geometry.model.getElementIDs() + slip_rates = geometry.get_slip_rates(args.elements) + for id in slip_rates.keys(): + sys.stdout.write("{} {}\n".format(id,slip_rates[id])) + + if args.slip_time_series: + if args.elements is None: raise "Must specify element ids, e.g. --elements 0 1 2" + if args.min_year is None: args.min_year = 0.0 + if args.max_year is None: args.max_year = 20.0 + if args.dt is None: args.dt = 0.5 # Unit is decimal years + if args.use_sections is not None: + if len(args.use_sections) > 1: + section_name = "" + for sec in args.use_sections: + section_name += geometry.model.section(sec).name()+", " + else: + section_name = geometry.model.section(args.use_sections[0]).name()+", " + time_series = geometry.get_slip_time_series(events, elements=args.elements, min_year=args.min_year, max_year=args.max_year, DT=args.dt) + if len(time_series.keys()) < 10: + labels = time_series.keys()+[""] + else: + labels = [None for each in range(len(time_series.keys())+1)] + x_data = [list(np.arange(args.min_year+args.dt, args.max_year+args.dt, args.dt)) for key in time_series.keys()]+[[args.min_year,args.max_year]] + linewidths = [0.8 for key in time_series.keys()]+[1] + styles = ["-" for key in time_series.keys()]+["--"] + y_data = time_series.values()+[[0,0]] + if args.use_sections is not None: + plot_title = "Slip time series for {}from years {} to {} with step {}\n{}".format(section_name, args.min_year,args.max_year,args.dt,args.event_file.split("/")[-1]) + else: + plot_title = "Slip time series for {} elements, from years {} to {} with step {}\n{}".format(len(args.elements), args.min_year,args.max_year,args.dt,args.event_file.split("/")[-1]) + filename = SaveFile().diagnostic_plot(args.event_file, "slip_time_series", min_year=args.min_year, max_year=args.max_year, min_mag=args.min_magnitude) + BasePlotter().multi_line_plot(x_data, y_data, labels, linewidths, plot_title, "sim time [years]", "cumulative slip [m]", "", filename, linestyles=styles) + + if args.event_kml: + if args.event_id is None or args.event_file is None or args.model_file is None: + raise "Must specify an event to plot with --event_id and provide an --event_file and a --model_file." + else: + event = events._events[args.event_id] + filename = SaveFile().event_kml_plot(args.event_file, args.event_id) + geometry.model.write_event_kml(filename, event) + # Generate stress plots if args.stress_elements: # TODO: check that stress_set is valid StressHistoryPlot().plot(stress_set, args.stress_elements) if args.num_sweeps: - filename = SaveFile().diagnostic_plot(args.event_file, "num_sweeps") + filename = SaveFile().diagnostic_plot(args.event_file, "num_sweeps", min_year=args.min_year, max_year=args.max_year, min_mag=args.min_magnitude) DiagnosticPlot().plot_number_of_sweeps(events, filename) if args.event_shear_stress: - filename = SaveFile().diagnostic_plot(args.event_file, "shear_stress") + filename = SaveFile().diagnostic_plot(args.event_file, "shear_stress", min_year=args.min_year, max_year=args.max_year, min_mag=args.min_magnitude) DiagnosticPlot().plot_shear_stress_changes(events, filename) if args.event_normal_stress: - filename = SaveFile().diagnostic_plot(args.event_file, "normal_stress") + filename = SaveFile().diagnostic_plot(args.event_file, "normal_stress", min_year=args.min_year, max_year=args.max_year, min_mag=args.min_magnitude) DiagnosticPlot().plot_normal_stress_changes(events, filename) if args.event_mean_slip: - filename = SaveFile().diagnostic_plot(args.event_file, "mean_slip") + filename = SaveFile().diagnostic_plot(args.event_file, "mean_slip", min_year=args.min_year, max_year=args.max_year, min_mag=args.min_magnitude) DiagnosticPlot().plot_mean_slip(events, filename) # Validate data if requested diff --git a/doc/graphics/Event_0_plot.jpg b/doc/graphics/Event_0_plot.jpg new file mode 100644 index 00000000..152ec488 Binary files /dev/null and b/doc/graphics/Event_0_plot.jpg differ diff --git a/doc/graphics/Event_4406_crop_logo.png b/doc/graphics/Event_4406_crop_logo.png new file mode 100644 index 00000000..764cc59c Binary files /dev/null and b/doc/graphics/Event_4406_crop_logo.png differ diff --git a/doc/graphics/VQ_Logo_fancy.png b/doc/graphics/VQ_Logo_fancy.png new file mode 100644 index 00000000..2762ea2f Binary files /dev/null and b/doc/graphics/VQ_Logo_fancy.png differ diff --git a/doc/graphics/greens_fitting.png b/doc/graphics/greens_fitting.png new file mode 100644 index 00000000..2957f1d9 Binary files /dev/null and b/doc/graphics/greens_fitting.png differ diff --git a/doc/graphics/mag_mean_slip_minMag_5-0_events_ALLCAL2_VQmeshed_3km_25kyr_dyn1-3_GreensLimits_VQstressDrops0-3_AllowNegSlips.png b/doc/graphics/mag_mean_slip_minMag_5-0_events_ALLCAL2_VQmeshed_3km_25kyr_dyn1-3_GreensLimits_VQstressDrops0-3_AllowNegSlips.png new file mode 100644 index 00000000..e893db49 Binary files /dev/null and b/doc/graphics/mag_mean_slip_minMag_5-0_events_ALLCAL2_VQmeshed_3km_25kyr_dyn1-3_GreensLimits_VQstressDrops0-3_AllowNegSlips.png differ diff --git a/doc/graphics/mag_rupt_area_minMag_5-0_events_ALLCAL2_VQmeshed_3km_25kyr_dyn1-3_GreensLimits_VQstressDrops0-3_AllowNegSlips.png b/doc/graphics/mag_rupt_area_minMag_5-0_events_ALLCAL2_VQmeshed_3km_25kyr_dyn1-3_GreensLimits_VQstressDrops0-3_AllowNegSlips.png new file mode 100644 index 00000000..e59f49a9 Binary files /dev/null and b/doc/graphics/mag_rupt_area_minMag_5-0_events_ALLCAL2_VQmeshed_3km_25kyr_dyn1-3_GreensLimits_VQstressDrops0-3_AllowNegSlips.png differ diff --git a/doc/graphics/slip_time_series_diagnostic_yearMin0_yearMax100_SAF-Mendo_Offs_events_ALLCAL2_VQmeshed_3km_dyn1-3_newStressDrops_GreenLimits_stressDropFactor0-3.png b/doc/graphics/slip_time_series_diagnostic_yearMin0_yearMax100_SAF-Mendo_Offs_events_ALLCAL2_VQmeshed_3km_dyn1-3_newStressDrops_GreenLimits_stressDropFactor0-3.png new file mode 100644 index 00000000..a3b49b39 Binary files /dev/null and b/doc/graphics/slip_time_series_diagnostic_yearMin0_yearMax100_SAF-Mendo_Offs_events_ALLCAL2_VQmeshed_3km_dyn1-3_newStressDrops_GreenLimits_stressDropFactor0-3.png differ diff --git a/doc/graphics/slip_time_series_diagnostic_yearMin0_yearMax200_events_ALLCAL2_VQmeshed_3km_dyn1-3_newStressDrops_GreenLimits_stressDropFactor0-3.png b/doc/graphics/slip_time_series_diagnostic_yearMin0_yearMax200_events_ALLCAL2_VQmeshed_3km_dyn1-3_newStressDrops_GreenLimits_stressDropFactor0-3.png new file mode 100644 index 00000000..43fafd27 Binary files /dev/null and b/doc/graphics/slip_time_series_diagnostic_yearMin0_yearMax200_events_ALLCAL2_VQmeshed_3km_dyn1-3_newStressDrops_GreenLimits_stressDropFactor0-3.png differ diff --git a/doc/vq.bib b/doc/vq.bib index 77a64605..24a168e8 100755 --- a/doc/vq.bib +++ b/doc/vq.bib @@ -158,7 +158,7 @@ @ARTICLE{2005AGUFMNG21A..03Y adsnote = {Provided by the SAO/NASA Astrophysics Data System} } -@article{Turcotte:2007up, +@article{Turcotte2007up, author = {Turcotte, D and Holliday, J and Rundle, J}, title = {{BASS, an alternative to ETAS}}, journal = {Geophys. Res. Lett}, diff --git a/doc/vq.tex b/doc/vq.tex index de2e0a9f..55d54b65 100755 --- a/doc/vq.tex +++ b/doc/vq.tex @@ -41,8 +41,8 @@ %%%%%%%%%%%%%%%%%%%%%%%%%%% Lyx Title Page \title{Virtual Quake User Manual} -\author{Eric M. Heien\\Michael K. Sachs\\Kasey W. Schultz\\Donald L. Turcotte\\John B. Rundle\\\\ \copyright University of California, Davis\\ -Version 1.1.0} +\author{Kasey W. Schultz\\Eric M. Heien\\Michael K. Sachs\\Mark R. Yoder\\John B. Rundle\\Donald L. Turcotte\\\\ \copyright University of California, Davis\\ +Version 2.0.0} %\date{\noindent \today} @@ -76,27 +76,29 @@ %COLOR AND CODENAME BLOCK% \begin{center} -\resizebox{\textwidth}{!}{\colorbox -% FILL: color of code name text box -% e.g. blue -{blue}{\fontfamily{\rmdefault}\selectfont \textcolor{white} { -% FILL: name of the code -% You may want to add \hspace to both sides of the codename to better center it, such as: -% \newcommand{\codename}{\hspace{0.1in}CodeName\hspace{0.1in}} -\hspace{0.02in}Virtual Quake\hspace{0.0in} -}}} +%\resizebox{\textwidth}{!}{\colorbox +%% FILL: color of code name text box +%% e.g. blue +%{blue}{\fontfamily{\rmdefault}\selectfont \textcolor{white} { +%% FILL: name of the code +%% You may want to add \hspace to both sides of the codename to better center it, such as: +%% \newcommand{\codename}{\hspace{0.1in}CodeName\hspace{0.1in}} +%\hspace{0.02in}Virtual Quake\hspace{0.0in} +%}}} +\includegraphics[width=\textwidth]{graphics/VQ_Logo_fancy.png} \end{center} %MAIN PICTURE% -\begin{textblock*}{0in}(0.8in,-0.05in) +\begin{textblock*}{0in}(0.0in,-0.0in) % FILL: image height % e.g. height=6.5in \begin{center} -\vspace{.5in} -\includegraphics[height=5.0in] +\vspace{.8in} +\includegraphics[height=4.52in] % FILL: image file name % e.g. cover_image.png -{graphics/overlay_fields_109382.png} +%{graphics/overlay_fields_109382.png} +{graphics/4406_cover.jpg} \end{center} \end{textblock*} @@ -105,7 +107,7 @@ \hfill{\Huge \fontfamily{\sfdefault}\selectfont User Manual \\ % FILL: manual version % e.g. 1.0 -\raggedleft \huge \fontfamily{\sfdefault}\selectfont Version {1.1}\\} +\raggedleft \huge \fontfamily{\sfdefault}\selectfont Version {2.0}\\} %AUTHOR(S) & WEBSITE% \null @@ -115,7 +117,7 @@ % FILL: author list % e.g. Author One\\Author Two\\Author Three\\ % be sure to have a newline (\\) after the final author -Eric M. Heien~~~Michael K. Sachs~~~Kasey W. Schultz\\Donald L. Turcotte~~~John B. Rundle\\ +Kasey W. Schultz~~~Eric M. Heien~~~Michael K. Sachs\\Mark R. Yoder~~~John B. Rundle~~~Donald L. Turcotte\\ } {\fontfamily{\sfdefault}\selectfont www.geodynamics.org} @@ -599,7 +601,7 @@ \subsection{Event Transition Time} $A$ based on model element strengths. During long term stress accumulation, an element is defined to fail at time $t_{f}$ when $CFF^{A}(t_{f})=0$, which is referred to as static failure. At this point the simulation -changes to the rupture event model described below. +changes to the rupture model described below. Given the change in stress over time it is relatively straightforward to calculate the time to failure for an element. Since effective long @@ -624,7 +626,7 @@ \subsection{Event Transition Time} when the next element will fail. -\section{Rupture Event Model\label{sec:Event_Model}} +\section{Rupture Model\label{sec:Event_Model}} A rupture event in VQ is comprised of multiple sweeps. During each sweep one or more elements fail and one or more elements slip. The sweeps continue as long as there are elements that fail - once no more elements fail the event is complete. @@ -664,7 +666,7 @@ \section{Rupture Event Model\label{sec:Event_Model}} This relationship between element slips during a sweep is determined for each ruptured element $A$ as: \begin{eqnarray} -0 = \sum (T_s^{AB} - \mu_s^B T_n^{AB}) s_B +\Delta\sigma_A - CFF_A = \sum (T_s^{AB} - \mu_s^B T_n^{AB}) s_B \end{eqnarray} Since non-ruptured blocks do not change their slip, they are excluded from the system. This system is then solved using Gaussian elimination. @@ -694,6 +696,43 @@ \section{Rupture Event Model\label{sec:Event_Model}} sweeps in an event, but do not accumulate additional stress after failing. +\subsection{Stress Drops \label{sec:stress_drops}} + +The stress drops for the fault elements are computed via Wells and Coppersmith 1994 scaling relations along with an analytical solution for the shear stress change per unit slip on a vertical rectangular fault. The simulation is pretty sensitive to these values; large stress drops tend to produce larger earthquakes with longer periods between events, while smaller stress drops produce many more small earthquakes. We provide a tuning parameter that globally adjusts the stress drops, \texttt{sim.friction.stress\_drop\_factor}. The default value is 0.3 and it's logarithmically scaled, and increasing this value by 0.1 multiplies the stress drops by 1.4, a 40\% increase. + +We begin by using scaling relations to determine a different characteristic magnitude $M^{char}$ and a characteristic slip $\Delta s^{char}$ for each element in a fault based on the fault's geometry from Wells and Coppersmith 1994 scaling relations + +\begin{eqnarray} + M^{char} = 4.07 + 0.98 \log_{10} (A) + \mbox{sim.friction.stress\_drop\_factor} +\end{eqnarray} + +where $A$ is the surface area of the fault in $km^2$. We then use the definition of moment magnitude to determine the mean slip $\Delta s^{char}$ + +\begin{eqnarray} + \Delta s^{char} = \frac{ 10^{ \frac{3}{2} (M^{char} + 10.7) } }{ 10^7 \mu A } +\end{eqnarray} + +where $\mu$ is the rigidity (shear modulus) of the elastic half space, and $A$ here is the area of the fault in $m^2$. + +We then use this slip to determine the characteristic stress drop for each fault element + +\begin{eqnarray} + \Delta\sigma = - \frac{2 \mu \Delta s^{char}}{(1-\nu) \pi R} \left( (1-\nu) \frac{L}{W} + \frac{W}{L} \right). +\end{eqnarray} + +where $\nu$ is the Poisson ratio of the elastic half space, $W$ is the down-dip width of the fault, $L$ is the length of the fault in the direction of slip, and + +\begin{eqnarray} + R = \sqrt{L^2 + W^2} . +\end{eqnarray} + + +\subsection{ETAS Aftershock Model} + +We include an Epidemic Type Aftershock Model in Virtual Quake, specifically the BASS variant of the model described in \cite{Turcotte2007up}. By setting the simulation parameter \texttt{sim.bass.max\_generations = 1} you can turn on the aftershock model. The aftershock model draws location from a prescribed Omori-type distribution, is mapped onto the closest simulation element to the randomly drawn location, and the aftershock is processed as a regular simulated earthquake and changes the local stress field accordingly. The aftershocks parameters are described in Section \ref{sec:etas_params}. + + + \section{Simulation Flow} \begin{figure}[th] @@ -1206,15 +1245,15 @@ \section{Advanced Usage} \subsection{Tuning Parameters\label{sec:tuning_parameters}} Fault simulations must initially be tuned to correctly simulate actual -earthquakes. The primary tuning parameter in Virtual Quake is +earthquakes. The primary tuning parameters in Virtual Quake are the dynamic triggering factor $\eta$, -defined in Section \ref{sec:Event_Model}. The dynamic triggering +defined in Section \ref{sec:Event_Model}, and the stress drop factor defined in Section \ref{sec:stress_drops}. The dynamic triggering factor is used to encourage rupture propagation during a simulated earthquake. This parameter acts to tune the rupture and slipping properties of faults without requiring field measurements of each fault's physical properties, a result of VQ's abstraction and generality. Currently, this parameter is set globally for the -fault model. +fault model. For parameter usage, see Section \ref{sec:friction_params}. %Figure \ref{fig:tuning_parameters} shows a comparison %of simulation results for a range of these parameter values. % @@ -1477,6 +1516,7 @@ \subsection{Simulation Parameter File} the simulator where the fault model files are located, sets various simulation variables, and specifies the output of the simulation. An example parameter file template is located in \texttt{examples/}. +For all possible simulation parameters, see Appendix \ref{cha:Appendix-A:-Input}. \subsection{Producing a Fault Model\label{sec:Using-Mesher}} @@ -1535,6 +1575,40 @@ \subsubsection{Parameter File} Edit the parameter file so the simulation parameters are set correctly for the current run. +\subsection{Bounding the Greens Functions} +In some cases, the mesher produces elements that have some small percentage of overlap. This overlap can cause the Greens functions (interaction coefficients) to take on extreme/non-physical values. If your simulation is producing anomalous earthquakes, or behaving oddly, try the following. + +After using a test simulation to save your Greens functions to a file (e.g. ``greens\_values.h5"), you can use the script in vq/PyVQ/pyvq/betas/ called greens.py. Open a \textbf{python environment} from this directory, or in a python script execute the following to fit a Gaussian profile to the Greens functions, plot the distribution of values, and return the 1 standard deviation bounds. The values shown in example output below are the 1 standard deviation bounds (max, min) for the Greens functions. Issue the plt.show() command to reveal the plot of the Greens functions and the best fit Gaussian, as shown in Figure \ref{fig:greens_fit}. + +\begin{verbatim} + $ python + Python 2.7.9 (default, Dec 13 2014, 15:13:49) + [GCC 4.2.1 Compatible Apple LLVM 6.0 (clang-600.0.56)] on darwin + Type "help", "copyright", "credits" or "license" for more information. + >>> from matplotlib import pyplot as plt + >>> import greens + >>> greens.plot_greens_hists(greens_fname="path/to/greens_values.h5", shear_normal="shear") + [....some output omitted...] + Greens range for 1 sigma ... x=[ 3661623.85 -831572.04] ... + >>> plt.show() + >>> greens.plot_greens_hists(greens_fname="path/to/greens_values.h5", shear_normal="normal") + [....some output omitted...] + Greens range for 1 sigma ... x=[ 61623.85 -31572.04] ... + >>> plt.show() +\end{verbatim} + +To implement these bounds for the shear stress Greens functions, in your simulation parameter file set the \\ \texttt{sim.greens.shear\_offdiag\_max} equal the the 1 standard deviation maximum value given by greens.py (3661623.85 in the example above), and set \texttt{sim.greens.shear\_offdiag\_min} equal the the 1 standard deviation minimum value (-831572.04 in the example above). You can do the same for the normal stress bounds, setting \texttt{sim.greens.normal\_offdiag\_max} and \texttt{sim.greens.normal\_offdiag\_min}. All parameters are given in Appendix \ref{sec:greens_parameters}. + + +\begin{figure} + \centering + \includegraphics[width=0.8\textwidth]{graphics/greens_fitting.png} + \caption{Example histogram of shear stress Greens function values with best fit Gaussian and bounds.} + \label{fig:greens_fit} +\end{figure} + + + \section{Single Element Tutorial\label{sec:Tutorial_single}} @@ -1876,15 +1950,13 @@ \subsection{Simulation Statistics Plots} To generate a frequency magnitude plot from your simulation file (hdf5), located at ``path/to/sim\_file.h5'', simply execute the shell command (in one line): \begin{verbatim} -$ python path/to/vq/pyvq/pyvq.py --event_file path/to/sim_file.h5 --event_file_type - 'hdf5' --plot_freq_mag 1 +$ python path/to/vq/pyvq/pyvq.py --event_file path/to/sim_file.h5 --plot_freq_mag 1 \end{verbatim} If you generated a text file from your simulation, simply use instead: \begin{verbatim} -$ python path/to/vq/pyvq/pyvq.py --event_file path/to/sim_file.txt --event_file_type - 'text' --plot_freq_mag 1 +$ python path/to/vq/pyvq/pyvq.py --event_file path/to/sim_file.txt --plot_freq_mag 1 \end{verbatim} This saves the plot to a file with the filename chosen from the name of the simulation file and the type of plot being drawn, Figure \ref{fig:freq_mag}. @@ -1898,32 +1970,74 @@ \subsection{Simulation Statistics Plots} To plot both the magnitude-rupture area and magnitude-mean slip distributions (Figure \ref{fig:mag_area_slip}) simply execute: \begin{verbatim} -$ python path/to/vq/pyvq/pyvq.py --event_file path/to/sim_file.h5 --event_file_type - 'hdf5' --plot_mag_rupt_area --plot_mag_mean_slip +$ python path/to/vq/pyvq/pyvq.py --event_file path/to/sim_file.h5 + --plot_mag_rupt_area --plot_mag_mean_slip \end{verbatim} \begin{figure} \centering{} -\includegraphics[width=0.49\textwidth]{graphics/mag_rupt_area_example.png} -\includegraphics[width=0.49\textwidth]{graphics/mag_mean_slip_example.png} +\includegraphics[width=0.49\textwidth]{graphics/mag_rupt_area_minMag_5-0_events_ALLCAL2_VQmeshed_3km_25kyr_dyn1-3_GreensLimits_VQstressDrops0-3_AllowNegSlips.png} +\includegraphics[width=0.49\textwidth]{graphics/mag_mean_slip_minMag_5-0_events_ALLCAL2_VQmeshed_3km_25kyr_dyn1-3_GreensLimits_VQstressDrops0-3_AllowNegSlips.png} \caption{\label{fig:mag_area_slip} \textbf{Left}: Example magnitude vs rupture area distribution. \textbf{Right}: Example magnitude vs mean slip distribution. Compare these to Wells and Coppersmith 1994 relations.} \end{figure} +\subsection{Overview Plots/Functions} +To list the largest $N$ (an integer) earthquakes (by magnitude) and print their properties, execute the following: + +\begin{verbatim} +$ python path/to/vq/pyvq/pyvq.py --event_file path/to/sim_file.h5 --summary N +\end{verbatim} + +To plot all the statistical plots for a simulation, e.g. the frequency-magnitude, the mean slip vs .magnitude, rupture area vs. magnitude, and other distributions, execute the following: + +\begin{verbatim} +$ python path/to/vq/pyvq/pyvq.py --event_file path/to/sim_file.h5 --all_stat_plots +\end{verbatim} + +or you can execute the following to plot all the statistical plots and some additional time series plots: + +\begin{verbatim} +$ python path/to/vq/pyvq/pyvq.py --event_file path/to/sim_file.h5 --diagnostics +\end{verbatim} + +\subsection{Time Series Plots} +To plot the slip time series of element 20 from year 0 to year 100 shown in figure \ref{fig:slip_time_series} (use the --dt parameter to set the time step in units of years): + +\begin{figure} +\centering{} +\includegraphics[width=0.49\textwidth]{graphics/slip_time_series_diagnostic_yearMin0_yearMax200_events_ALLCAL2_VQmeshed_3km_dyn1-3_newStressDrops_GreenLimits_stressDropFactor0-3.png} +\includegraphics[width=0.49\textwidth]{graphics/slip_time_series_diagnostic_yearMin0_yearMax100_SAF-Mendo_Offs_events_ALLCAL2_VQmeshed_3km_dyn1-3_newStressDrops_GreenLimits_stressDropFactor0-3.png} +\caption{\label{fig:slip_time_series} \textbf{Left}: Slip time series for element \#20 \textbf{Right}: Slip time series for all elements in section 0, SAF-Mendo\_Offs.} +\end{figure} + +\begin{verbatim} +$ python path/to/vq/pyvq/pyvq.py --event_file path/to/sim_file.h5 --slip_time_series + --elements 20 --min_year 0 --max_year 200 --dt 0.5 +\end{verbatim} + +And to plot the slip time series for all elements in section 0 (a.k.a. "fault 0") from year 0 to year 100 shown in figure \ref{fig:slip_time_series}: + +\begin{verbatim} +$ python path/to/vq/pyvq/pyvq.py --event_file path/to/sim_file.h5 --slip_time_series + --use_sections 0 --min_year 0 --max_year 100 --dt 0.5 +\end{verbatim} + + \subsection{Earthquake Probability Plots} To plot the conditional probabilities (Figure \ref{fig:probabilities}) of an earthquake with magnitude $\geq 7$ occurring on sections [4,5,6,7,8,9,10,11,12,13] (as defined in your fault model) as a function of time since the last earthquake on those sections: \begin{verbatim} -$ python path/to/vq/pyvq/pyvq.py --event_file path/to/sim_file.h5 --event_file_type - 'hdf5' --plot_cond_prob_vs_t --model_file path/to/model.txt --model_file_type +$ python path/to/vq/pyvq/pyvq.py --event_file path/to/sim_file.h5 + --plot_cond_prob_vs_t --model_file path/to/model.txt --model_file_type 'text' --min_magnitude 7.0 --use_sections 4 5 6 7 8 9 10 11 12 13 \end{verbatim} To plot the expected waiting times until the next earthquake with magnitude $\geq 7$ occurring on sections [4,5,6,7,8,9,10,11,12,13] as a function of time since the last earthquake on those sections: \begin{verbatim} -$ python path/to/vq/pyvq/pyvq.py --event_file path/to/sim_file.h5 --event_file_type - 'hdf5' --plot_waiting_times --model_file path/to/model.txt --model_file_type +$ python path/to/vq/pyvq/pyvq.py --event_file path/to/sim_file.h5 + --plot_waiting_times --model_file path/to/model.txt --model_file_type 'text' --min_magnitude 7.0 --use_sections 4 5 6 7 8 9 10 11 12 13 \end{verbatim} @@ -1936,7 +2050,29 @@ \subsection{Earthquake Probability Plots} \end{figure} -\subsection{Field Plots} +\subsection{Plotting Data on a Map and Event Field Plots} + +To plot the traces of your faults on a simple map: + +\begin{verbatim} +$ python path/to/vq/pyvq/pyvq.py --event_file path/to/sim_file.h5 --model_file + path/to/model.txt --traces +\end{verbatim} + +To create a Google Earth file (KML) of the elements involved in an event, colored by their event slips (shown in Figure \ref{fig:event_kml}): + +\begin{verbatim} +$ python path/to/vq/pyvq/pyvq.py --event_file path/to/sim_file.h5 --model_file + path/to/model.txt --event_kml --event_id 0 +\end{verbatim} + +\begin{figure} +\centering{} +\includegraphics[width=.9\textwidth]{graphics/Event_0_plot.jpg} +\caption{\label{fig:event_kml} Slips for event 0. The triggering element is shown with the white place marker. Color scale is white (almost zero slip) to red (max slip).} +\end{figure} + + QuakeLib uses Green's functions to compute co-seismic fields given a fault geometry and a slip distribution. For displacements we use Okada's equations and for gravity, potential, and geoid height changes we use Okubo's equations. @@ -1957,12 +2093,12 @@ \subsection{Field Plots} To compute the gravity changes and InSAR interferogram for event \#210 from your simulation (Figures \ref{fig:event_insar} and \ref{fig:event_gravity}): \begin{verbatim} -$ python path/to/vq/pyvq/pyvq.py --event_file path/to/sim_file.h5 --event_file_type - 'hdf5' --model_file path/to/model.txt --model_file_type 'text' --field_plot +$ python path/to/vq/pyvq/pyvq.py --event_file path/to/sim_file.h5 --model_file + path/to/model.txt --model_file_type 'text' --field_plot --field_type 'gravity' --event_id 210 -$ python path/to/vq/pyvq/pyvq.py --event_file path/to/sim_file.h5 --event_file_type - 'hdf5' --model_file path/to/model.txt --model_file_type 'text' --field_plot +$ python path/to/vq/pyvq/pyvq.py --event_file path/to/sim_file.h5 --model_file + path/to/model.txt --model_file_type 'text' --field_plot --field_type 'insar' --event_id 210 \end{verbatim} @@ -1979,6 +2115,17 @@ \subsection{Field Plots} \caption{\label{fig:event_gravity} Simulated gravity changes for magnitude 7.26 earthquake on the San Andreas Fault. } \end{figure} +To evaluate an event field at specified lat/lon points, you must provide the file containing lines of lat/lon pairs (--lld\_file), the result is written to a file: + +\begin{verbatim} +$ python path/to/vq/pyvq/pyvq.py --event_file path/to/sim_file.h5 --model_file + path/to/model.txt --model_file_type 'text' --field_eval + --field_type 'gravity' --event_id 210 --lld_file path/to/lld_file.txt +\end{verbatim} + + +\newpage + \section{PyVQ Plotting Arguments Grouped by Functionality\label{sec:pyvq_params}} This section summarizes the current functionality of the pyvq script (section \ref{sec:pyvq}). These parameters do not have default values, and are not used unless specified. @@ -1990,8 +2137,6 @@ \subsection{Model parameters} \hline \texttt{\small{--event\_file}} & The simulation output file.\tabularnewline \hline -\texttt{\small{--event\_file\_type}} & Either 'hdf5' or 'text'.\tabularnewline -\hline \texttt{\small{--sweep\_file}} & The simulation sweep file, required only if simulation file is text.\tabularnewline \hline \texttt{\small{--model\_file}} & The fault model file.\tabularnewline @@ -2000,7 +2145,7 @@ \subsection{Model parameters} \hline \end{tabular} -\subsection{Subsetting Parameters} +\subsection{Subsetting/Filtering Parameters} These parameters are used to select subsets of the faults or subsets of events. \\ \noindent % @@ -2008,10 +2153,20 @@ \subsection{Subsetting Parameters} \hline \texttt{\small{--use\_sections}} & Select a subset of the fault sections from the model file. To select sections 23, 43, and 101 it would be ``--use\_sections 23 43 101''.\tabularnewline \hline +\texttt{\small{--use\_trigger\_sections}} & Select a subset of events with triggering elements from the specified fault sections from the model file.\tabularnewline +\hline \texttt{\small{--min\_magnitude}} & Select a minimum magnitude for earthquakes to be analyzed. To select earthquakes only above magnitude 5.5 it would be ``--min\_magnitude 5.5''.\tabularnewline \hline \texttt{\small{--max\_magnitude}} & Select a maximum magnitude for earthquakes to be analyzed. Can be used in conjunction with --min\_magnitude to specify a range of magnitudes. \tabularnewline \hline +\texttt{\small{--min\_year}} & Select a minimum year for a time interval subset. Can be used in conjunction with --max\_year to specify a simulation time range. \tabularnewline +\hline +\texttt{\small{--max\_year}} & Select a maximum year for a time interval subset. Can be used in conjunction with --min\_year to specify a simulation time range. \tabularnewline +\hline +\texttt{\small{--min\_slip}} & Select a minimum mean slip threshold for plotting event data, units are meters. \tabularnewline +\hline +\texttt{\small{--max\_slip}} & Select a maximum mean slip threshold for plotting event data, units are meters. \tabularnewline +\hline \end{tabular} @@ -2020,6 +2175,8 @@ \subsection{\noindent Statistical Plotting Parameters} \noindent % \begin{tabular}{|>{\raggedright}p{2in}|>{\raggedright}p{4.1in}|} \hline +\texttt{\small{--summary n}} & List the n earthquakes with the largest magnitude, also prints their event mean slip, rupture area, etc.\tabularnewline +\hline \texttt{\small{--plot\_freq\_mag n}} & Frequency-Magnitude. n=1 Normal plot, n=2 adds a Gutenberg-Richter line with b=1, n=3 adds UCERF2 observed seismicity rates in California, n=4 adds both the b=1 line and observed California rates.\tabularnewline \hline \texttt{\small{--plot\_mag\_rupt\_area}} & Magnitude vs Rupture Area scaling. Compare to Wells and Coppersmith 1994.\tabularnewline @@ -2030,6 +2187,27 @@ \subsection{\noindent Statistical Plotting Parameters} \hline \texttt{\small{--wc94}} & Add Wells and Coppersmith 1994 scaling relation to Mean Slip or Rupture Area plots.\tabularnewline \hline +\texttt{\small{--diagnostics}} & Plot many different statistical plots and time series plots to characterize a simulation.\tabularnewline +\hline +\end{tabular} + +\subsection{\noindent Time Series Plotting Parameters} + +\noindent % +\begin{tabular}{|>{\raggedright}p{2in}|>{\raggedright}p{4.1in}|} +\hline +\texttt{\small{--slip\_time\_series}} & Plot the slip time series for a single element up to an entire fault, use with --elements or --use\_sections.\tabularnewline +\hline +\texttt{\small{--dt}} & Time step in decimal years for slip time series plots.\tabularnewline +\hline +\texttt{\small{--num\_sweeps}} & Plot, for each event, the number of event sweeps as a function of simulation time.\tabularnewline +\hline +\texttt{\small{--event\_mean\_slip}} & Plot, for each event, the mean slip as a function of simulation time.\tabularnewline +\hline +\texttt{\small{--event\_shear\_stress}} & Plot, for each event, the fractional change in shear stress as a function of simulation time.\tabularnewline +\hline +\texttt{\small{--event\_normal\_stress}} & Plot, for each event, the fractional change in normal stress as a function of simulation time.\tabularnewline +\hline \end{tabular} \subsection{\noindent Probability Plotting Parameters} @@ -2070,6 +2248,14 @@ \subsection{\noindent Field Plotting Parameters} \hline \texttt{\small{--angles}} & The observing angles for displacements and InSAR, the field is projected along this direction. Must specify azimuth and elevation angles in degrees. E.g. ``--angles 20.5 66.2''. \tabularnewline \hline +\texttt{\small{--traces}} & Plot the fault traces on a simple map.\tabularnewline +\hline +\texttt{\small{--field\_eval}} & Evaluate an event field at given lat/lon points. Must specify --field\_type, --lld\_file. Results are saved to file\tabularnewline +\hline +\texttt{\small{--lld\_file}} & Text file containing lat/lon pairs in columns. Used for field evaluations.\tabularnewline +\hline +\texttt{\small{--event\_kml}} & Create and save a Google Earth (KML) file that includes elements that slip during the specified event, colored by the amount of event slip. Must also specify --event\_id.\tabularnewline +\hline \end{tabular} @@ -2123,7 +2309,7 @@ \subsection{\noindent System Parameters} \end{tabular} -\subsection{Friction parameters} +\subsection{Friction parameters \label{sec:friction_params}} \begin{tabular}{|>{\raggedright}p{2.3in}|>{\raggedright}p{3.8in}|} \hline @@ -2132,10 +2318,14 @@ \subsection{Friction parameters} and result in larger earthquakes, while lower values indicate ruptures are less likely to propagate and result in smaller earthquakes.\tabularnewline \hline +\texttt{\small{sim.friction.compute\_stress\_drops}} & Boolean. Default is True. If true, compute the stress drops from the fault geometry and scaling relations. If false, must specify stress drops in the model file. See section \ref{sec:stress_drops}.\tabularnewline +\hline +\texttt{\small{sim.friction.stress\_drop\_factor}} & Default is 0.3, this factor is a stress drop multiplier on a log scale. Increasing by 0.1 implies a global stress drop increase of 40\%. Must use this with sim.friction.compute\_stress\_drops. See section \ref{sec:stress_drops}.\tabularnewline +\hline \end{tabular} -\subsection{Green's function parameters} +\subsection{Green's function parameters \label{sec:greens_parameters}} \begin{tabular}{|>{\raggedright}p{2.3in}|>{\raggedright}p{3.8in}|} \hline @@ -2167,6 +2357,22 @@ \subsection{Green's function parameters} \hline \texttt{\small{sim.greens.sample\_distance = 1000.0}} & When calculating the Green's function, take samples at this minimum distance between samples. This allows better convergence between models with few large elements and models with many small elements. If the element size is smaller than this value, it has no effect.\tabularnewline \hline +\texttt{\small{sim.greens.shear\_offdiag\_min}} & Double, no default value. If specified, this truncates the shear interaction Greens function values to some minimum. Sometimes required if the mesher puts fault elements too close together.\tabularnewline +\hline +\texttt{\small{sim.greens.shear\_offdiag\_max}} & Double, no default value. If specified, this truncates the shear interaction Greens function values to some maximum. Sometimes required if the mesher puts fault elements too close together.\tabularnewline +\hline +\texttt{\small{sim.greens.shear\_diag\_max}} & Double, no default value. If specified, this truncates the shear self-stress Greens function values to some maximum. Sometimes required if the mesher puts fault elements too close together.\tabularnewline +\hline +\texttt{\small{sim.greens.shear\_diag\_min}} & Double, no default value. If specified, this truncates the shear self-stress Greens function values to some minimum. Sometimes required if the mesher puts fault elements too close together.\tabularnewline +\hline +\texttt{\small{sim.greens.normal\_offdiag\_min}} & Double, no default value. If specified, this truncates the normal interaction Greens function values to some minimum. Sometimes required if the mesher puts fault elements too close together.\tabularnewline +\hline +\texttt{\small{sim.greens.normal\_offdiag\_max}} & Double, no default value. If specified, this truncates the normal interaction Greens function values to some maximum. Sometimes required if the mesher puts fault elements too close together.\tabularnewline +\hline +\texttt{\small{sim.greens.normal\_diag\_max}} & Double, no default value. If specified, this truncates the normal self-stress Greens function values to some maximum. Sometimes required if the mesher puts fault elements too close together.\tabularnewline +\hline +\texttt{\small{sim.greens.normal\_diag\_min}} & Double, no default value. If specified, this truncates the normal self-stress Greens function values to some minimum. Sometimes required if the mesher puts fault elements too close together.\tabularnewline +\hline \end{tabular} @@ -2211,7 +2417,7 @@ \subsection{File name parameters} %\end{tabular} -\subsection{BASS (Branching Aftershock Sequence) model parameters} +\subsection{BASS (Branching Aftershock Sequence) model parameters \label{sec:etas_params}} \begin{tabular}{|>{\raggedright}p{2.3in}|>{\raggedright}p{3.8in}|} \hline @@ -2232,7 +2438,7 @@ \subsection{BASS (Branching Aftershock Sequence) model parameters} \texttt{\small{sim.bass.d = 300}}{\small \par} -\texttt{\small{sim.bass.q = 1.35}} & Different parameters for BASS model. See paper \cite{Turcotte:2007up} +\texttt{\small{sim.bass.q = 1.35}} & Different parameters for BASS model. See paper \cite{Turcotte2007up} for details. Mm: minimum magnitude diff --git a/examples/CMakeLists.txt b/examples/CMakeLists.txt index e2271424..d49145ed 100755 --- a/examples/CMakeLists.txt +++ b/examples/CMakeLists.txt @@ -40,8 +40,8 @@ FOREACH(NPROC ${NUM_PROCS}) FOREACH(TAPER_IND RANGE ${NUM_TAPER}) LIST(GET TAPER_METHODS ${TAPER_IND} TAPER) - LIST(GET METHOD_MEANSLIP ${TAPER_IND} MEANSLIP) - LIST(GET METHOD_INTEREVENT ${TAPER_IND} INTEREVENT) + #LIST(GET METHOD_MEANSLIP ${TAPER_IND} MEANSLIP) + #LIST(GET METHOD_INTEREVENT ${TAPER_IND} INTEREVENT) FOREACH(RES ${RESOLUTIONS}) SET(TEST_DIR ${CMAKE_CURRENT_BINARY_DIR}/PROCS${NPROC}/${TAPER}/) FILE(MAKE_DIRECTORY ${TEST_DIR}) @@ -84,23 +84,24 @@ FOREACH(NPROC ${NUM_PROCS}) ENVIRONMENT "PYTHONPATH=${QUAKELIB_BINARY_DIR}/python/") # Check that the mean slip is near our expectations - IF (NOT ${MEANSLIP} EQUAL 0) - ADD_TEST(NAME test_slip_P${NPROC}_${TAPER}_${RES} WORKING_DIRECTORY ${TEST_DIR} - COMMAND ${PYTHON_EXECUTABLE} ${CHECK_SCRIPT} --event_file ${EVENT_FILE} --sweep_file ${SWEEP_FILE} --mean_slip ${MEANSLIP}) - SET_TESTS_PROPERTIES (test_slip_P${NPROC}_${TAPER}_${RES} PROPERTIES - DEPENDS run_P${NPROC}_${TAPER}_${RES} - ENVIRONMENT "PYTHONPATH=${QUAKELIB_BINARY_DIR}/python/") - ENDIF (NOT ${MEANSLIP} EQUAL 0) + # Schultz: For now removing until we come up with better tests + #IF (NOT ${MEANSLIP} EQUAL 0) + # ADD_TEST(NAME test_slip_P${NPROC}_${TAPER}_${RES} WORKING_DIRECTORY ${TEST_DIR} + # COMMAND ${PYTHON_EXECUTABLE} ${CHECK_SCRIPT} --event_file ${EVENT_FILE} --sweep_file ${SWEEP_FILE} --mean_slip ${MEANSLIP}) + # SET_TESTS_PROPERTIES (test_slip_P${NPROC}_${TAPER}_${RES} PROPERTIES + # DEPENDS run_P${NPROC}_${TAPER}_${RES} + # ENVIRONMENT "PYTHONPATH=${QUAKELIB_BINARY_DIR}/python/") + #ENDIF (NOT ${MEANSLIP} EQUAL 0) # Check that the mean interevent time is near our expectations - IF (NOT ${INTEREVENT} EQUAL 0) - - ADD_TEST(NAME test_interevent_P${NPROC}_${TAPER}_${RES} WORKING_DIRECTORY ${TEST_DIR} - COMMAND ${PYTHON_EXECUTABLE} ${CHECK_SCRIPT} --event_file ${EVENT_FILE} --sweep_file ${SWEEP_FILE} --mean_interevent ${INTEREVENT}) - SET_TESTS_PROPERTIES (test_interevent_P${NPROC}_${TAPER}_${RES} PROPERTIES - DEPENDS run_P${NPROC}_${TAPER}_${RES} - ENVIRONMENT "PYTHONPATH=${QUAKELIB_BINARY_DIR}/python/") - ENDIF(NOT ${INTEREVENT} EQUAL 0) + # Schultz: For now removing until we come up with better tests + #IF (NOT ${INTEREVENT} EQUAL 0) + # ADD_TEST(NAME test_interevent_P${NPROC}_${TAPER}_${RES} WORKING_DIRECTORY ${TEST_DIR} + # COMMAND ${PYTHON_EXECUTABLE} ${CHECK_SCRIPT} --event_file ${EVENT_FILE} --sweep_file ${SWEEP_FILE} --mean_interevent ${INTEREVENT}) + # SET_TESTS_PROPERTIES (test_interevent_P${NPROC}_${TAPER}_${RES} PROPERTIES + # DEPENDS run_P${NPROC}_${TAPER}_${RES} + # ENVIRONMENT "PYTHONPATH=${QUAKELIB_BINARY_DIR}/python/") + #ENDIF(NOT ${INTEREVENT} EQUAL 0) ENDFOREACH(RES ${RESOLUTIONS}) ENDFOREACH(TAPER_IND RANGE ${NUM_TAPER}) ENDFOREACH(NPROC ${NUM_PROCS}) diff --git a/quakelib/src/QuakeLibIO.cpp b/quakelib/src/QuakeLibIO.cpp index cb33dfa1..7424a5f2 100755 --- a/quakelib/src/QuakeLibIO.cpp +++ b/quakelib/src/QuakeLibIO.cpp @@ -1901,6 +1901,21 @@ int quakelib::ModelWorld::write_file_kml(const std::string &file_name) { out_file << "\n"; out_file << "\n"; + + // Loop thru elements to compute min/max slip rates for color bounds + double max_rate = 0; + double min_rate = 0; + for (eit=_elements.begin(); eit!=_elements.end(); ++eit) { + max_rate = fmax( max_rate, eit->second.slip_rate()); + min_rate = fmin( min_rate, eit->second.slip_rate()); + } + + /////// Bounds for color in blue to red range, white in the middle + double y_min = 0; + double y_max = 255; + double x_min = min_rate; + double x_max = max_rate; + // Go through the sections for (fit=_sections.begin(); fit!=_sections.end(); ++fit) { // And output the elements for each section @@ -1927,16 +1942,41 @@ int quakelib::ModelWorld::write_file_kml(const std::string &file_name) { lld[3] = lld[2]; lld[2] = c.convert2LatLon(xyz[2]+(xyz[1]-xyz[0])); } - - // Output the KML format polygon for this section + + // Compute blue to red color (RGB) + // Keep red scale (min is white, max is red) so red=255 + // Blue and green are equal and vary from (255 for min vals to 0 for max vals) + int blue, green; + int red = y_max; + if (eit->second.slip_rate() == 0) { + blue = 0; + green = 0; + red = 0; + } else { + int interp_color = (int) linear_interp(eit->second.slip_rate(), x_min, x_max, y_min, y_max); + blue = y_max - interp_color; + green = blue; + } + + // Output the KML format polygon for this element out_file << "\t\t\n"; out_file << "\t\t\n"; out_file << "Fault name: " << fit->second.name() << "\n"; + out_file << "Element #: " << eit->second.id() << "\n"; out_file << "Slip rate: " << c.m_per_sec2cm_per_yr(eit->second.slip_rate()) << " cm/year\n"; out_file << "Rake: " << c.rad2deg(eit->second.rake()) << " degrees\n"; out_file << "Aseismic: " << eit->second.aseismic() << "\n"; out_file << "\t\t\n"; - out_file << "\t\t\t#baseStyle\n"; + out_file << "\t\t\t\n"; + //out_file << "\t\t\t#baseStyle\n"; out_file << "\t\t\t\n"; out_file << "\t\t\t\t0\n"; out_file << "\t\t\t\trelativeToGround\n"; @@ -1967,6 +2007,198 @@ int quakelib::ModelWorld::write_file_kml(const std::string &file_name) { return 0; } +// Schultz: Using these functions for color interpolation and RGB to HEX color conversions +double quakelib::ModelWorld::linear_interp(const double &x, const double &x_min, const double &x_max, const double &y_min, const double &y_max) const { + return ((y_max - y_min)/(x_max - x_min) * (x - x_min)) + y_min; +} + +char* quakelib::ModelWorld::rgb2hex(const int r, const int g, const int b) const { + // Returning ABGR to work with KML format + char *buf; + size_t sz; + sz = snprintf(NULL, 0, "FF%02X%02X%02X", b,g,r); + buf = (char *)malloc(sz + 1); + snprintf(buf, sz+1, "FF%02X%02X%02X", b,g,r); + return buf; +} + +int quakelib::ModelWorld::write_event_kml(const std::string &file_name, const quakelib::ModelEvent &event) { + std::ofstream out_file; + std::map::const_iterator fit; + std::map::const_iterator eit; + LatLonDepth min_bound, max_bound, center; + Vec<3> min_xyz, max_xyz; + double dx, dy, range, mag, mean_slip, ev_year; + unsigned int i, trigger, ev_num; + ElementIDSet involved_elements; + ElementIDSet::iterator it; + + out_file.open(file_name.c_str()); + + get_bounds(min_bound, max_bound); + center = LatLonDepth(max_bound.lat() - (max_bound.lat()-min_bound.lat())/2, + max_bound.lon() - (max_bound.lon()-min_bound.lon())/2); + Conversion c(center); + min_xyz = c.convert2xyz(min_bound); + max_xyz = c.convert2xyz(max_bound); + dx = max_xyz[0]-min_xyz[0]; + dy = max_xyz[1]-min_xyz[1]; + range = fmax(dx, dy) * (1.0/tan(c.deg2rad(30))); + + double max_depth = fabs(min_bound.altitude()); + + out_file << "\n"; + out_file << "\n"; + out_file << "\n"; + out_file << "\n"; + out_file << "\t" << center.lat() << "\n"; + out_file << "\t" << center.lon() << "\n"; + out_file << "\t0\n"; + out_file << "\t" << range << "\n"; + out_file << "\t0\n"; + out_file << "\t0\n"; + out_file << "\tabsolute\n"; + out_file << "\n"; + out_file << "\n"; + + //// Get relevant event info + involved_elements = event.getInvolvedElements(); + trigger = event.getEventTrigger(); + mean_slip = event.calcMeanSlip(); + ev_year = event.getEventYear(); + ev_num = event.getEventNumber(); + + // Write event info above triggering element + out_file << "\t\n"; + out_file << "\t\t" << "Event "<< ev_num << "\n"; + out_file << "\t\t\n"; + out_file << "Trigger: Element #"<< trigger << " \n"; + out_file << "Mean slip [m] "<< mean_slip << " \n"; + out_file << "Event year "<< ev_year << " \n"; + out_file << "\t\t\n"; + out_file << "\t\t#sectionLabel\n"; + out_file << "\t\t\n"; + out_file << "\t\t\t1\n"; + out_file << "\t\t\trelativeToGround\n"; + // Find the deepest element for this section + UIndex best_vertex; + double min_altitude = DBL_MAX, cur_alt; + ModelVertex vert; + + eit = _elements.find(trigger); + for (i=0; i<3; ++i) { + cur_alt = _vertices.find(eit->second.vertex(i))->second.lld().altitude(); + if (cur_alt < min_altitude) { + min_altitude = cur_alt; + best_vertex = eit->second.vertex(i); + } + } + + vert = _vertices.find(best_vertex)->second; + out_file << "\t\t\t" << vert.lld().lon() << "," << vert.lld().lat() << "," << max_depth + vert.lld().altitude() << "\n"; + out_file << "\t\t\n"; + out_file << "\t\n"; + + + // Loop thru elements to compute min/max slip rates for color bounds + double max_slip = -DBL_MAX; + double min_slip = DBL_MAX; + for (it=involved_elements.begin(); it!=involved_elements.end(); ++it) { + max_slip = fmax( max_slip, event.getEventSlip(*it)); + min_slip = fmin( min_slip, event.getEventSlip(*it)); + } + + /////// Bounds for color in blue to red range, white in the middle + double y_min = -255; + double y_max = 255; + double x_max = fmax(fabs(max_slip), fabs(min_slip)); + double x_min = -x_max; + + out_file << "\n"; + // Go through the involved elements + for (it=involved_elements.begin(); it!=involved_elements.end(); ++it) { + LatLonDepth lld[4]; + unsigned int i, npoints; + + for (i=0; i<3; ++i) { + std::map::const_iterator vit; + vit = _vertices.find( _elements.find(*it)->second.vertex(i) ); + lld[i] = vit->second.lld(); + } + + // If this is a quad element, calculate the 4th implicit point + if (eit->second.is_quad()) { + Vec<3> xyz[3]; + + for (i=0; i<3; ++i) xyz[i] = c.convert2xyz(lld[i]); + + lld[3] = lld[2]; + lld[2] = c.convert2LatLon(xyz[2]+(xyz[1]-xyz[0])); + } + + // Compute blue to red color (RGB) + int blue, green, red; + int interp_color = (int) linear_interp(event.getEventSlip(*it), x_min, x_max, y_min, y_max); + if (interp_color > 0) { + // More red + red = y_max; + blue = y_max - interp_color; + green = blue; + } else { + blue = y_max; + red = y_max - fabs(interp_color); + green = red; + } + + // Output the KML format polygon for this element + out_file << "\t\t\n"; + out_file << "\t\t\n"; + out_file << "Element #: " << *it << "\n"; + out_file << "Event slip [m]: " << event.getEventSlip(*it) << "\n"; + out_file << "\t\t\n"; + out_file << "\t\t\t\n"; + out_file << "\t\t\t\n"; + out_file << "\t\t\t\t0\n"; + out_file << "\t\t\t\trelativeToGround\n"; + out_file << "\t\t\t\t\n"; + out_file << "\t\t\t\t\t\n"; + out_file << "\t\t\t\t\t\t\n"; + npoints = (eit->second.is_quad() ? 4 : 3); + + for (i=0; i\n"; + out_file << "\t\t\t\t\t\n"; + out_file << "\t\t\t\t\n"; + out_file << "\t\t\t\n"; + out_file << "\t\t\n"; + + } + + out_file << "\n"; + out_file << "\n"; + out_file << "\n"; + + out_file.close(); + + return 0; +} + + void quakelib::ModelSection::apply_remap(const ModelRemapping &remap) { _data._id = remap.get_section_id(_data._id); } diff --git a/quakelib/src/QuakeLibIO.h b/quakelib/src/QuakeLibIO.h index a1aa24d1..07068107 100755 --- a/quakelib/src/QuakeLibIO.h +++ b/quakelib/src/QuakeLibIO.h @@ -506,129 +506,6 @@ namespace quakelib { typedef std::set ElementIDSet; - class ModelWorld : public ModelIO { - private: - std::map _vertices; - std::map _elements; - std::map _sections; - LatLonDepth _base; - double _min_lat, _max_lat, _min_lon, _max_lon; - -#ifdef HDF5_FOUND - void read_section_hdf5(const hid_t &data_file); - void read_element_hdf5(const hid_t &data_file); - void read_vertex_hdf5(const hid_t &data_file); - - void write_section_hdf5(const hid_t &data_file) const; - void write_element_hdf5(const hid_t &data_file) const; - void write_vertex_hdf5(const hid_t &data_file) const; -#endif - - public: - ModelSection &new_section(void); - ModelElement &new_element(void); - ModelVertex &new_vertex(void); - - ModelSection §ion(const UIndex &ind) throw(std::domain_error); - ModelElement &element(const UIndex &ind) throw(std::domain_error); - ModelVertex &vertex(const UIndex &ind) throw(std::domain_error); - - siterator begin_section(void) { - return siterator(&_sections, _sections.begin()); - }; - siterator end_section(void) { - return siterator(&_sections, _sections.end()); - }; - - eiterator begin_element(const UIndex &fid=INVALID_INDEX) { - return eiterator(&_elements, _elements.begin(), fid); - }; - eiterator end_element(const UIndex &fid=INVALID_INDEX) { - return eiterator(&_elements, _elements.end(), fid); - }; - - UIndex next_section_index(void) const { - if (_sections.size()) return _sections.rbegin()->first+1; - else return 0; - }; - UIndex next_element_index(void) const { - if (_elements.size()) return _elements.rbegin()->first+1; - else return 0; - }; - UIndex next_vertex_index(void) const { - if (_vertices.size()) return _vertices.rbegin()->first+1; - else return 0; - }; - - size_t num_sections(void) const; - size_t num_faults(void) const; - size_t num_elements(const UIndex &fid=INVALID_INDEX) const; - size_t num_vertices(const UIndex &fid=INVALID_INDEX) const; - - LatLonDepth min_bound(const UIndex &fid=INVALID_INDEX) const; - LatLonDepth max_bound(const UIndex &fid=INVALID_INDEX) const; - - Vec<3> get_base(void) const { - return Vec<3>(_base.lat(), _base.lon(), _base.altitude()); - } - - std::vector get_latlon_bounds(void) const { - std::vector bounds(4); - bounds[0] = _min_lat; - bounds[1] = _max_lat; - bounds[2] = _min_lon; - bounds[3] = _max_lon; - return bounds; - } - - void insert(const ModelWorld &other_world); - void insert(const ModelSection &new_section); - void insert(const ModelElement &new_element); - void insert(const ModelVertex &new_vertex); - - void get_bounds(LatLonDepth &minimum, LatLonDepth &maximum) const; - - SimElement create_sim_element(const UIndex &element_id) const; - SlippedElement create_slipped_element(const UIndex &element_id) const; - - ElementIDSet getElementIDs(void) const; - ElementIDSet getVertexIDs(void) const; - ElementIDSet getSectionIDs(void) const; - - bool overwrite(const ModelRemapping &remap); - void apply_remap(const ModelRemapping &remap); - ModelRemapping remap_indices_contiguous(const UIndex &start_section_index, - const UIndex &start_element_index, - const UIndex &start_vertex_index) const; - ModelRemapping remove_duplicate_vertices_remap(void) const; - void clear(void); - - void reset_base_coord(const LatLonDepth &new_base); - - void create_section(std::vector &unused_trace_segments, - const std::vector &trace, - const LatLonDepth &base_coord, - const UIndex &fault_id, - const float &element_size, - const std::string §ion_name, - const std::string &taper_method, - const bool resize_trace_elements); - - int read_file_ascii(const std::string &file_name); - int write_file_ascii(const std::string &file_name) const; - - int read_file_trace_latlon(std::vector &unused_trace_segments, const std::string &file_name, const float &elem_size, const std::string &taper_method, const bool resize_trace_elements); - int write_file_trace_latlon(const std::string &file_name, const float &depth_along_dip); - - int read_file_hdf5(const std::string &file_name); - int write_file_hdf5(const std::string &file_name) const; - - int write_file_kml(const std::string &file_name); - - int read_files_eqsim(const std::string &geom_file_name, const std::string &cond_file_name, const std::string &fric_file_name); - int write_files_eqsim(const std::string &geom_file_name, const std::string &cond_file_name, const std::string &fric_file_name); - }; - // Class recording data associated with a block that slipped during an event // mu is static, but we retain it for use in calculating magnitude. struct SweepData { @@ -1248,6 +1125,134 @@ namespace quakelib { int read_file_ascii(const std::string &stress_index_file_name, const std::string &stress_file_name); }; + + class ModelWorld : public ModelIO { + private: + std::map _vertices; + std::map _elements; + std::map _sections; + LatLonDepth _base; + double _min_lat, _max_lat, _min_lon, _max_lon; + +#ifdef HDF5_FOUND + void read_section_hdf5(const hid_t &data_file); + void read_element_hdf5(const hid_t &data_file); + void read_vertex_hdf5(const hid_t &data_file); + + void write_section_hdf5(const hid_t &data_file) const; + void write_element_hdf5(const hid_t &data_file) const; + void write_vertex_hdf5(const hid_t &data_file) const; +#endif + + public: + ModelSection &new_section(void); + ModelElement &new_element(void); + ModelVertex &new_vertex(void); + + ModelSection §ion(const UIndex &ind) throw(std::domain_error); + ModelElement &element(const UIndex &ind) throw(std::domain_error); + ModelVertex &vertex(const UIndex &ind) throw(std::domain_error); + + siterator begin_section(void) { + return siterator(&_sections, _sections.begin()); + }; + siterator end_section(void) { + return siterator(&_sections, _sections.end()); + }; + + eiterator begin_element(const UIndex &fid=INVALID_INDEX) { + return eiterator(&_elements, _elements.begin(), fid); + }; + eiterator end_element(const UIndex &fid=INVALID_INDEX) { + return eiterator(&_elements, _elements.end(), fid); + }; + + UIndex next_section_index(void) const { + if (_sections.size()) return _sections.rbegin()->first+1; + else return 0; + }; + UIndex next_element_index(void) const { + if (_elements.size()) return _elements.rbegin()->first+1; + else return 0; + }; + UIndex next_vertex_index(void) const { + if (_vertices.size()) return _vertices.rbegin()->first+1; + else return 0; + }; + + size_t num_sections(void) const; + size_t num_faults(void) const; + size_t num_elements(const UIndex &fid=INVALID_INDEX) const; + size_t num_vertices(const UIndex &fid=INVALID_INDEX) const; + + LatLonDepth min_bound(const UIndex &fid=INVALID_INDEX) const; + LatLonDepth max_bound(const UIndex &fid=INVALID_INDEX) const; + + Vec<3> get_base(void) const { + return Vec<3>(_base.lat(), _base.lon(), _base.altitude()); + } + + std::vector get_latlon_bounds(void) const { + std::vector bounds(4); + bounds[0] = _min_lat; + bounds[1] = _max_lat; + bounds[2] = _min_lon; + bounds[3] = _max_lon; + return bounds; + } + + void insert(const ModelWorld &other_world); + void insert(const ModelSection &new_section); + void insert(const ModelElement &new_element); + void insert(const ModelVertex &new_vertex); + + void get_bounds(LatLonDepth &minimum, LatLonDepth &maximum) const; + + SimElement create_sim_element(const UIndex &element_id) const; + SlippedElement create_slipped_element(const UIndex &element_id) const; + + ElementIDSet getElementIDs(void) const; + ElementIDSet getVertexIDs(void) const; + ElementIDSet getSectionIDs(void) const; + + bool overwrite(const ModelRemapping &remap); + void apply_remap(const ModelRemapping &remap); + ModelRemapping remap_indices_contiguous(const UIndex &start_section_index, + const UIndex &start_element_index, + const UIndex &start_vertex_index) const; + ModelRemapping remove_duplicate_vertices_remap(void) const; + void clear(void); + + void reset_base_coord(const LatLonDepth &new_base); + + void create_section(std::vector &unused_trace_segments, + const std::vector &trace, + const LatLonDepth &base_coord, + const UIndex &fault_id, + const float &element_size, + const std::string §ion_name, + const std::string &taper_method, + const bool resize_trace_elements); + + int read_file_ascii(const std::string &file_name); + int write_file_ascii(const std::string &file_name) const; + + int read_file_trace_latlon(std::vector &unused_trace_segments, const std::string &file_name, const float &elem_size, const std::string &taper_method, const bool resize_trace_elements); + int write_file_trace_latlon(const std::string &file_name, const float &depth_along_dip); + + int read_file_hdf5(const std::string &file_name); + int write_file_hdf5(const std::string &file_name) const; + + int write_file_kml(const std::string &file_name); + int write_event_kml(const std::string &file_name, const ModelEvent &event); + + int read_files_eqsim(const std::string &geom_file_name, const std::string &cond_file_name, const std::string &fric_file_name); + int write_files_eqsim(const std::string &geom_file_name, const std::string &cond_file_name, const std::string &fric_file_name); + + double linear_interp(const double &x, const double &x_min, const double &x_max, const double &y_min, const double &y_max) const; + char* rgb2hex(const int r, const int g, const int b) const; + }; + } #endif diff --git a/src/misc/GreensFunctions.cpp b/src/misc/GreensFunctions.cpp index 9abf9d8f..4eac1ae2 100755 --- a/src/misc/GreensFunctions.cpp +++ b/src/misc/GreensFunctions.cpp @@ -84,7 +84,12 @@ void GreensFuncCalcStandard::CalculateGreens(Simulation *sim) { for (int r=0; rnumLocalBlocks(); ++r) { for (int c=0; cgetGlobalBID(r); - sim->setGreens(global_r, c, ssh[global_r][c], snorm[global_r][c]); + //// Schultz, excluding zero slip rate elements from sim by setting Greens to zero + if (sim->getBlock(global_r).slip_rate()==0 || sim->getBlock(c).slip_rate()==0) { + sim->setGreens(global_r, c, 0, 0); + } else { + sim->setGreens(global_r, c, ssh[global_r][c], snorm[global_r][c]); + } } } } @@ -157,7 +162,12 @@ void GreensFuncFileParse::CalculateGreens(Simulation *sim) { greens_file_reader->getGreensVals(gid, in_shear_green, in_normal_green); for (j=0; jsetGreens(gid, j, in_shear_green[j], in_normal_green[j]); + //// Schultz, excluding zero slip rate elements from sim by setting Greens to zero + if (sim->getBlock(gid).slip_rate()==0 || sim->getBlock(j).slip_rate()==0) { + sim->setGreens(gid, j, 0, 0); + } else { + sim->setGreens(gid, j, in_shear_green[j], in_normal_green[j]); + } } }