Ctrl+K

Special Features of the Engine#

There are a few less frequently used features of the engine that are not documented in the general user’s manual, since their usage is quite specific. They are documented here.

Sensitivity analysis#

Running a sensitivity analysis study means to run multiple calculations by changing a parameter and to study how the results change. For instance, it is interesting to study the random seed dependency when running a calculation using sampling of the logic tree, or it is interesting to study the impact of the truncation level on the PoEs. The engine offers a special command to run a sensitivity analysis with respect to one (or even more than one) parameter. Consider for instance the MultiPointClassicalPSHA demo; suppose you are interesting in studying what happens when changing the truncation_level parameter; then you can run a sensitivity analysis as follows:

$ oq sensitivity_analysis job.ini truncation_level=[2,3] | bash

The engine will spawn two calculations with descriptions as follows:

Multipoint demo [truncation_level=2]
Multipoint demo [truncation_level=3]

The engine is also able to run sensitivity analysis on file parameters. For instance if you want to run a classical_risk calculation starting from three different hazard inputs you can write:

$ oq sensitivity_analysis job.ini hazard_curves_file=["hazard1.csv","hazard2.csv","hazard3.csv"] | bash

Notice that the current approach (since engine 3.21) works by generating bash code instead of directly running calculations on the local machine, as in previous versions of the engine. It means that if you are using a supercomputer, by configuring correctly the submit_cmd and other parameters in openquake.cfg, it is possible to use multiple nodes of your infrastructure. To see the generated code, just remove the “| bash” pipe at the end of the command.

Here is an example changing two parameters at the same time:

$ oq sensitivity_analysis job.ini truncation_level=[2,3] area_source_discretization=[10,20]
#!/bin/bash
job_id=(`oq create_jobs 4`)
oq run job.ini -p truncation_level=2 area_source_discretization=10 job_id=${job_id[0]}
oq run job.ini -p truncation_level=2 area_source_discretization=20 job_id=${job_id[1]}
oq run job.ini -p truncation_level=3 area_source_discretization=10 job_id=${job_id[2]}
oq run job.ini -p truncation_level=3 area_source_discretization=20 job_id=${job_id[3]}
oq collect_jobs ${job_id[@]}

The exact details of the script may change across versions of the engine.

Ruptures in CSV format#

Since engine v3.10 there is a way to serialize ruptures in CSV format. The command to give is:

$ oq extract "ruptures?min_mag=<mag>" <calc_id>`

For instance, assuming there is an event based calculation with ID 42, we can extract the ruptures in the datastore with magnitude larger than 6 with oq extract "ruptures?min_mag=6" 42: this will generate a CSV file. Then it is possible to run multi-rupture scenario calculations starting from that file by simply setting rupture_model_file = ruptures-min_mag=6_42.csv in the job.ini file. The format is provisional and may change in the future, but it will stay a CSV with JSON fields. Here is an example for a planar rupture, i.e. a rupture generated by a point source:

#,,,,,,,,,,"trts=['Active Shallow Crust']"
seed,mag,rake,lon,lat,dep,multiplicity,trt,kind,mesh,extra
24,5.050000E+00,0.000000E+00,0.08456,0.15503,5.000000E+00,1,Active Shallow Crust,ParametricProbabilisticRupture PlanarSurface,"[[[[0.08456, 0.08456, 0.08456, 0.08456]], [[0.13861, 0.17145, 0.13861, 0.17145]], [[3.17413, 3.17413, 6.82587, 6.82587]]]]","{""occurrence_rate"": 4e-05}"

The format is meant to support all kind of ruptures, including ruptures generated by simple and complex fault sources, characteristic sources, nonparametric sources and new kind of sources that could be introduced in the engine in the future. The header will be the same for all kind of ruptures that will be stored in the same CSV. Here is description of the fields as they are named now (engine v3.11):

seed

the random seed used to compute the GMFs generated by the rupture

mag

the magnitude of the rupture

rake

the rake angle of the rupture surface in degrees

lon

the longitude of the hypocenter in degrees

lat

the latitude of the hypocenter in degrees

dep

the depth of the hypocenter in km

multiplicity

the number of occurrences of the rupture (i.e. number of events)

trt

the tectonic region type of the rupture; must be consistent with the trts listed in the pre-header of the file

kind

a space-separated string listing the rupture class and the surface class used in the engine

mesh

3 times nested list with lon, lat, dep of the points of the discretized rupture geometry for each underlying surface

extra

extra parameters of the rupture as a JSON dictionary, for instance the rupture occurrence rate

Notice that using a CSV file generated with an old version of the engine is inherently risky: for instance if we changed the ParametricProbabilisticRupture class or the PlanarSurface classes in an incompatible way with the past, then a scenario calculation starting with the CSV would give different results in the new version of the engine. We never changed the rupture classes or the surface classes, but we changed the seed algorithm often, and that too would cause different numbers to be generated (hopefully, statistically consistent). A bug fix or change of logic in the calculator can also change the numbers across engine versions.

The minimum_distance parameter#

GMPEs often have a prescribed range of validity. In particular they may give unexpected results for points too close to ruptures. To avoid this problem the engine recognizes a minimum_distance parameter: if it is set, then for distances below the specified minimum distance, the GMPEs return the ground-motion value at the minimum distance. This avoids producing extremely large (and physically unrealistic) ground-motion values at small distances. The minimum distance is somewhat heuristic. It may be useful to experiment with different values of the minimum_distance, to see how the hazard and risk change.

The max_sites_disagg#

There is a parameter in the job.ini called max_sites_disagg, with a default value of 10. This parameter controls the maximum number of sites on which it is possible to run a disaggregation. If you need to run a disaggregation on a large number of sites you will have to increase that parameter. Notice that there are technical limits: trying to disaggregate 100 sites will likely succeed, trying to disaggregate 100,000 sites will most likely cause your system to go out of memory or out of disk space, and the calculation will be terribly slow. If you have a really large number of sites to disaggregate, you will have to split the calculation and it will be challenging to complete all the subcalculations.

The parameter max_sites_disagg is extremely important not only for disaggregation, but also for classical calculations. Depending on its value and then number of sites (N) your calculation can be in the few sites regime or the many sites regime.

In the few sites regime (N <= max_sites_disagg) the engine stores information for each rupture in the model (in particular the distances for each site) and therefore uses more disk space. The problem is mitigated since the engine uses a relatively aggressive strategy to collapse ruptures, but that requires more RAM available.

In the many sites regime (N > max_sites_disagg) the engine does not store rupture information (otherwise it would immediately run out of disk space, since typical hazard models have tens of millions of ruptures) and uses a much less aggressive strategy to collapse ruptures, which has the advantage of requiring less RAM.

Equivalent Epicenter Distance Approximation#

The equivalent epicenter distance approximation (reqv for short) was introduced in engine 3.2 to enable the comparison of the OpenQuake engine with time-honored Fortran codes using the same approximation.

You can enable it in the engine by adding a [reqv] section to the job.ini, like in our example in openquake/qa_tests_data/logictree/case_02/job.ini:

reqv_hdf5 = {'active shallow crust': 'lookup_asc.hdf5',
             'stable shallow crust': 'lookup_sta.hdf5'}

For each tectonic region type to which the approximation should be applied, the user must provide a lookup table in .hdf5 format containing arrays mags of shape M, repi of shape N and reqv of shape (M, N).

The examples in openquake/qa_tests_data/classical/case_2 will give you the exact format required. M is the number of magnitudes (in the examples there are 26 magnitudes ranging from 6.05 to 8.55) and N is the number of epicenter distances (in the examples ranging from 1 km to 1000 km).

Depending on the tectonic region type and rupture magnitude, the engine converts the epicentral distance repi into an equivalent distance by looking at the lookup table and use it to determine the rjb and rrup distances, instead of the regular routines. This means that within this approximation ruptures are treated as pointwise and not rectangular as the engine usually does.

Notice that the equivalent epicenter distance approximation only applies to ruptures coming from PointSources/AreaSources/MultiPointSources, fault sources are untouched.

Infer Unknown Basin Parameters#

The z1pt0 (depth to a shear wave velocity of 1000 m/s) and z2pt5 (depth to shear wave velocity of 2500 m/s) provided in the site model CSV are used in GMPEs to compute the basin term (if the GMPE includes one). In engine version 3.24 we introduce the ability to specify a sentinel value of -999 for either parameter in the site model to permit each GMPE considered within a calculation to compute z1pt0 or z2pt5 independelty using the vs30 values in the site model (most GMPEs using basin terms also provide a vs30 to z1pt0 or vs30 to z2pt5 relationship which can be used to estimate the required value for each site). For GMPEs which do use either z1pt0 or z2pt5, but do not provide in their associated publications a vs30 to z1pt0 or vs30 to z2pt5 relationship, we use the most appropriate of those available from other GMPEs (e.g. for the HassaniAtkinson2020 GMPE developed for application to Japan no vs30 to z2pt5 relationship is provided in the journal article, so we assume the CampbellBozorgnia2014 GMPE’s vs30 to z2pt5 relationship for Japan is the most appropriate to use here). It is strongly advisable that if you make use of this feature that you examine this behaviour within any GMPEs included in the GMC logic tree.

This capability was added as required for implementation of the 2023 Conterminous USA model developed by the USGS.

Disaggregation by source disagg_by_src#

Given a system of various sources affecting a specific site, one very common question to ask is: what are the more relevant sources, i.e. which sources contribute the most to the mean hazard curve? The engine is able to answer such question by setting the disagg_by_src flag in the job.ini file. When doing that, the engine saves in the datastore a 4-dimensional ArrayWrapper called mean_rates_by_src with dimensions (site ID, intensity measure type, intensity measure level, source ID). From that it is possible to extract the contribution of each source to the mean hazard curve (interested people should look at the code in the function check_disagg_by_src). The ArrayWrapper mean_rates_by_src can also be converted into a pandas DataFrame, then getting something like the following:

>> dstore['mean_rates_by_src'].to_dframe().set_index('src_id')
               site_id  imt  lvl         value
ASCTRAS407           0  PGA    0  9.703749e-02
IF-CFS-GRID03        0  PGA    0  3.720510e-02
ASCTRAS407           0  PGA    1  6.735009e-02
IF-CFS-GRID03        0  PGA    1  2.851081e-02
ASCTRAS407           0  PGA    2  4.546237e-02
...                ...  ...  ...           ...
IF-CFS-GRID03        0  PGA   17  6.830692e-05
ASCTRAS407           0  PGA   18  1.072884e-06
IF-CFS-GRID03        0  PGA   18  1.275539e-05
ASCTRAS407           0  PGA   19  1.192093e-07
IF-CFS-GRID03        0  PGA   19  5.960464e-07

The value field here is the probability of exceedence in the hazard curve. The lvl field is an integer corresponding to the intensity measure level in the hazard curve.

In engine 3.15 we introduced the so-called “colon convention” on source IDs: if you have many sources that for some reason should be collected together - for instance because they all account for seismicity in the same tectonic region, or because they are components of a same source but are split into separate sources by magnitude - you can tell the engine to collect them into one source in the mean_rates_by_src matrix. The trick is to use IDs with the same prefix, a colon, and then a numeric index. For instance, if you had 3 sources with IDs src_mag_6.65, src_mag_6.75, src_mag_6.85, fragments of the same source with different magnitudes, you could change their IDs to something like src:0, src:1, src:2 and that would reduce the size of the matrix mean_rates_by_src by 3 times by collecting together the contributions of each source. There is no restriction on the numeric indices to start from 0, so using the names src:665, src:675, src:685 would work too and would be clearer: the IDs should be unique, however.

If the IDs are not unique and the engine determines that the underlying sources are different, then an extension “semicolon + incremental index” is automatically added. This is useful when the hazard modeler wants to define a model where the more than one version of the same source appears in one source model, having changed some of the parameters, or when varied versions of a source appear in each branch of a logic tree. In that case, the modeler should use always the exact same ID (i.e. without the colon and numeric index): the engine will automatically distinguish the sources during the calculation of the hazard curves and consider them the same when saving the array mean_rates_by_src: you can see an example in the test qa_tests_data/classical/case_20/job_bis.ini in the engine code base. In that case the source_info dataset will list 7 sources CHAR1;0 CHAR1;1 CHAR1;2 COMFLT1;0 COMFLT1;1 SFLT1;0 SFLT1;1 but the matrix mean_rates_by_src will see only three sources CHAR1 COMFLT1 SFLT1 obtained by composing together the versions of the underlying sources.

In version 3.15 mean_rates_by_src was extended to work with mutually exclusive sources, i.e. for the Japan model. You can see an example in the test qa_tests_data/classical/case_27. However, the case of mutually exclusive ruptures - an example is the New Madrid cluster in the USA model - is not supported yet.

In some cases it is tricky to discern whether use of the colon convention or identical source IDs is appropriate. The following list indicates several possible cases that a user may encounter, and the appropriate approach to assigning source IDs. Note that this list includes the cases that have been tested so far, and is not a comprehensive list of all cases that may arise.

  1. Sources in the same source group/source model are scaled alternatives of each other. For example, this occurs when for a given source, epistemic uncertainties such as occurrence rates or geometries are considered, but the modeller has pre-scaled the rates rather than including the alternative hypothesis in separate logic tree branches.

    Naming approach: identical IDs.

  2. Sources in different files are alternatives of each other, e.g. each is used in a different branch of the source model logic tree.

    Naming approach: identical IDs.

  3. A source is defined in OQ by numerous sources, either in the same file or different ones. For example, one could have a set of non-parametric sources, each with many ruptures, that are grouped together into single files by magnitude. Or, one could have many point sources that together represent the seismicity from one source.

    Naming approach: colon convention

  4. One source consists of many mutually exclusive sources, as in qa_tests_data/classical/case_27.

    Naming approach: colon convention

Cases 1 and 2 could include include more than one source typology, as in qa_tests_data/classical/case_79.

NB: disagg_by_src can be set to true only if the ps_grid_spacing approximation is disabled. The reason is that the ps_grid_spacing approximation builds effective sources which are not in the original source model, thus breaking the connection between the values of the matrix and the original sources.

On this page
Show Source