Added benchmark script

examples/f_120_sumatra.yaml

@@ -21,4 +21,4 @@
         return r/(0.7*Vs)+(0.081.0*r_crit/(0.7.0*Vs))*(1.0/(1.0-pow(r/r_crit, 2.0))-1.0);
       }
       return 1000000000.0;
-[s_xx, s_yy, s_zz, s_xy, s_yz, s_xz]: !Include ../examples/120_initial_stress.yaml
+[s_xx, s_yy, s_zz, s_xy, s_yz, s_xz, plastCo]: !Include ../examples/120_initial_stress.yaml

examples/f_16_scec.yaml

@@ -19,3 +19,4 @@ components: !SCECFile
         s_xy: return s_strike;
         s_yz: return s_dip;
     [mu_s, mu_d, d_c, cohesion, forced_rupture_time]: !IdentityMap {}
+    [s_normal, s_strike, s_dip]: !IdentityMap {}

tools/benchmark.py

@@ -0,0 +1,40 @@
+#!/usr/bin/env python3
+
+import argparse
+import os
+import subprocess
+
+cmdLineParser = argparse.ArgumentParser()
+cmdLineParser.add_argument('-N', type=int, required=True)
+cmdLineParser.add_argument('--easicube', type=str, default='../build/easicube')
+cmdLineParser.add_argument('--example_dir', type=str, default='../examples')
+args = cmdLineParser.parse_args()
+
+class Example(object):
+    def __init__(self, name, frm, to):
+        self.name = name
+        self.frm = frm
+        self.to = to
+
+examples = [
+    Example('0_constant', [0.0, 0.0, 0.0], [1.0, 1.0, 1.0]),
+    Example('2_prem', [0.0, 0.0, 0.0], [3678298.565, 3678298.565, 3678298.565]),
+    Example('3_layered_linear', [0.0, 0.0, -4000.0], [1.0, 1.0, 0.0]),
+    Example('5_function', [0.0, 0.0, 0.0], [1.0, 25.0, 1.0]),
+    Example('26_function', [-10.0, -1000.0, -100000.0], [10.0, 0.0, 100000.0]),
+    Example('33_layered_constant', [0.0, -1000.0, 0.0], [1.0, 1000.0, 1.0]),
+    Example('f_16_scec', [-5.0, 0.0, -1000.0], [1000.0, 1.0, 5.0]),
+    Example('f_120_sumatra', [0.0, 0.0, -50000.0], [50000.0, 100000.0, 0.0])
+]
+
+print('example,npoints,tsetup,teval')
+
+for example in examples:
+    print(example.name, end=',', flush=True)
+    spec = '{{from: [{}], to: [{}], N: [{}]}}'.format(
+                ','.join(map(str, example.frm)),
+                ','.join(map(str, example.to)),
+                ','.join(map(str, [args.N]*len(example.frm)))
+            )
+    fileName = os.path.join(args.example_dir, example.name + '.yaml')
+    subprocess.run([args.easicube, spec, fileName])

tools/easicube.cpp

@@ -4,6 +4,7 @@
 #include <fstream>
 #include <vector>
 #include <string>
+#include <chrono>
 #include <yaml-cpp/yaml.h>
 #include <netcdf.h>
 #include "easi/YAMLParser.h"
@@ -23,28 +24,36 @@ int main(int argc, char** argv)
     std::cerr << "Usage: easicube \"{from: [...], to: [...], N: [...], parameters: [...], group: ...}\" <model_file> [<output_basename>]" << std::endl;
     std::cerr << "  from, to: corner points of cube" << std::endl;
     std::cerr << "  N: number of grid points in each dimension" << std::endl;
-    std::cerr << "  parameters: parameters that shall be queried from the model" << std::endl;
-    std::cerr << "  group: group number" << std::endl;
+    std::cerr << "  parameters: parameters that shall be queried from the model (optional)" << std::endl;
+    std::cerr << "  group: group number (optional)" << std::endl;
     return -1;
   }
   std::string modelFile(argv[2]);
   std::string basename;
   if (argc >= 4) {
     basename = std::string(argv[3]);
-  } else {
-    basename = modelFile;
   }
-
+
+  std::vector<double> from, to;
+  std::vector<unsigned> N;
   YAML::Node config = YAML::Load(argv[1]);
-  std::vector<double> from = config["from"].as<std::vector<double>>();
-  std::vector<double> to = config["to"].as<std::vector<double>>();
-  std::vector<unsigned> N = config["N"].as<std::vector<unsigned>>();
+  try {
+    from = config["from"].as<std::vector<double>>();
+    to = config["to"].as<std::vector<double>>();
+    N = config["N"].as<std::vector<unsigned>>();
+  } catch (std::exception& e) {
+    std::cerr << "Error while parsing required arguments from, to, and N: " << e.what() << std::endl;
+    return -1;
+  }
   unsigned group = 0;
   if (config["group"]) {
     group = config["group"].as<unsigned>();
   }
-  std::vector<std::string> parameters = config["parameters"].as<std::vector<std::string>>();  
-
+  std::vector<std::string> parameters;
+  if (config["parameters"]) {
+    parameters = config["parameters"].as<std::vector<std::string>>();
+  }
+
   if (from.size() != 3 || to.size() != 3 || N.size() != 3) {
     std::cerr << "from, to, and N must have 3 dimensions" << std::endl;
     return -1;
@@ -60,7 +69,12 @@ int main(int argc, char** argv)
   for (unsigned d = 0; d < dim; ++d) {
     h.push_back( (to[d]-from[d]) / (N[d]-1) );
   }
-
+
+  auto parseStart = std::chrono::high_resolution_clock::now();
+  easi::YAMLParser parser(dim);
+  easi::Component* model = parser.parse(modelFile);
+  auto parseEnd = std::chrono::high_resolution_clock::now();
+
   easi::Query query(nGP, dim);
   for (unsigned k = 0; k < N[2]; ++k) {
     for (unsigned j = 0; j < N[1]; ++j) {
@@ -74,9 +88,10 @@ int main(int argc, char** argv)
     }
   }
 
-  easi::YAMLParser parser(dim);
-  easi::Component* model = parser.parse(modelFile);
-
+  if (parameters.empty()) {
+    auto supplied = model->suppliedParameters();
+    parameters = std::vector<std::string>(supplied.begin(), supplied.end());
+  }
   easi::ArraysAdapter adapter;
   std::vector<std::vector<double>> material(parameters.size());
   auto it = material.begin();
@@ -85,93 +100,102 @@ int main(int argc, char** argv)
     adapter.addBindingPoint(p, it->data());
     ++it;
   }
-
-  model->evaluate(query, adapter);  
+
+  auto evalStart = std::chrono::high_resolution_clock::now();
+  model->evaluate(query, adapter);
+  auto evalEnd = std::chrono::high_resolution_clock::now();
   delete model;
+
+  std::cout << parameters.size() << "," <<
+               nGP << "," <<
+               std::chrono::duration_cast<std::chrono::nanoseconds>(parseEnd-parseStart).count() * 1.0e-9 << "," <<
+               std::chrono::duration_cast<std::chrono::nanoseconds>(evalEnd-evalStart).count()   * 1.0e-9 << std::endl;
 
   /// Write xdmf
-
-  std::ofstream ofs(basename + ".xdmf", std::ofstream::out);
-
-  ofs << "<?xml version=\"1.0\" ?>\n"
-      << "<!DOCTYPE Xdmf SYSTEM \"Xdmf.dtd\" []>\n"
-      << "<Xdmf xmlns:xi=\"http://www.w3.org/2003/XInclude\" Version=\"2.2\">\n"
-      << "  <Domain>\n"
-      << "    <Grid Name=\"Model\" GridType=\"Uniform\">\n"
-      << "      <Topology TopologyType=\"3DCORECTMesh\" Dimensions=\"" << N[2] << " " << N[1] << " " << N[0] << "\"/>\n"
-      << "      <Geometry GeometryType=\"ORIGIN_DXDYDZ\">\n"
-      << "        <DataItem Name=\"Origin\" Dimensions=\"3\" NumberType=\"Float\" Precision=\"8\" Format=\"XML\">\n"
-      << "          " << from[2] << " " << from[1] << " " << from[0] << "\n"
-      << "        </DataItem>\n"
-      << "        <DataItem Name=\"Spacing\" Dimensions=\"3\" NumberType=\"Float\" Precision=\"8\" Format=\"XML\">\n"
-      << "          " << h[2] << " " << h[1] << " " << h[0] << "\n"
-      << "        </DataItem>\n"
-      << "      </Geometry>\n";
-  for (auto const& p : parameters) {
-    ofs << "      <Attribute Name=\"" << p << "\" AttributeType=\"Scalar\" Center=\"Node\">\n"
-        << "        <DataItem Format=\"HDF\" NumberType=\"Float\" Precision=\"8\" Dimensions=\"" << N[2] << " " << N[1] << " " << N[0] << "\">" << basename << ".nc:/" << p << "</DataItem>\n"
-        << "      </Attribute>\n";
-  }
-  ofs << "    </Grid>\n"
-      << "  </Domain>\n"
-      << "</Xdmf>\n";
 
-  ofs.close();
-
-  /// Write netcdf
-
-  int stat;  /* return status */
-  int ncid;  /* netCDF id */
-
-  /* dimension ids */
-  int x_dim;
-  int y_dim;
-  int z_dim;
-
-  /* dimension lengths */
-  size_t x_len = N[0];
-  size_t y_len = N[1];
-  size_t z_len = N[2];
-
-  /* enter define mode */
-  stat = nc_create(std::string(basename + ".nc").c_str(), NC_NETCDF4, &ncid);
-  check_err(stat,__LINE__,__FILE__);
-
-  /* define dimensions */
-  stat = nc_def_dim(ncid, "x", x_len, &x_dim);
-  check_err(stat,__LINE__,__FILE__);
-  stat = nc_def_dim(ncid, "y", y_len, &y_dim);
-  check_err(stat,__LINE__,__FILE__);
-  stat = nc_def_dim(ncid, "z", z_len, &z_dim);
-  check_err(stat,__LINE__,__FILE__);
-
-  /* define variables */
-
-  std::unordered_map<std::string, int> ids;
-  for (auto const& p : parameters) {
-    int dims[3];
-    dims[0] = z_dim;
-    dims[1] = y_dim;
-    dims[2] = x_dim;
+  if (!basename.empty()) {
+    std::ofstream ofs(basename + ".xdmf", std::ofstream::out);
+
+    ofs << "<?xml version=\"1.0\" ?>\n"
+        << "<!DOCTYPE Xdmf SYSTEM \"Xdmf.dtd\" []>\n"
+        << "<Xdmf xmlns:xi=\"http://www.w3.org/2003/XInclude\" Version=\"2.2\">\n"
+        << "  <Domain>\n"
+        << "    <Grid Name=\"Model\" GridType=\"Uniform\">\n"
+        << "      <Topology TopologyType=\"3DCORECTMesh\" Dimensions=\"" << N[2] << " " << N[1] << " " << N[0] << "\"/>\n"
+        << "      <Geometry GeometryType=\"ORIGIN_DXDYDZ\">\n"
+        << "        <DataItem Name=\"Origin\" Dimensions=\"3\" NumberType=\"Float\" Precision=\"8\" Format=\"XML\">\n"
+        << "          " << from[2] << " " << from[1] << " " << from[0] << "\n"
+        << "        </DataItem>\n"
+        << "        <DataItem Name=\"Spacing\" Dimensions=\"3\" NumberType=\"Float\" Precision=\"8\" Format=\"XML\">\n"
+        << "          " << h[2] << " " << h[1] << " " << h[0] << "\n"
+        << "        </DataItem>\n"
+        << "      </Geometry>\n";
+    for (auto const& p : parameters) {
+      ofs << "      <Attribute Name=\"" << p << "\" AttributeType=\"Scalar\" Center=\"Node\">\n"
+          << "        <DataItem Format=\"HDF\" NumberType=\"Float\" Precision=\"8\" Dimensions=\"" << N[2] << " " << N[1] << " " << N[0] << "\">" << basename << ".nc:/" << p << "</DataItem>\n"
+          << "      </Attribute>\n";
+    }
+    ofs << "    </Grid>\n"
+        << "  </Domain>\n"
+        << "</Xdmf>\n";
+
+    ofs.close();
 
-    stat = nc_def_var(ncid, p.c_str(), NC_DOUBLE, 3, dims, &ids[p]);
+    /// Write netcdf
+
+    int stat;  /* return status */
+    int ncid;  /* netCDF id */
+
+    /* dimension ids */
+    int x_dim;
+    int y_dim;
+    int z_dim;
+
+    /* dimension lengths */
+    size_t x_len = N[0];
+    size_t y_len = N[1];
+    size_t z_len = N[2];
+
+    /* enter define mode */
+    stat = nc_create(std::string(basename + ".nc").c_str(), NC_NETCDF4, &ncid);
     check_err(stat,__LINE__,__FILE__);
-  }
 
-  /* leave define mode */
-  stat = nc_enddef (ncid);
-  check_err(stat,__LINE__,__FILE__);
+    /* define dimensions */
+    stat = nc_def_dim(ncid, "x", x_len, &x_dim);
+    check_err(stat,__LINE__,__FILE__);
+    stat = nc_def_dim(ncid, "y", y_len, &y_dim);
+    check_err(stat,__LINE__,__FILE__);
+    stat = nc_def_dim(ncid, "z", z_len, &z_dim);
+    check_err(stat,__LINE__,__FILE__);
 
-  /* assign variable data */
-
-  it = material.begin();
-  for (auto const& p : parameters) {
-    stat = nc_put_var_double(ncid, ids[p], it->data());
-    ++it;
-  }
+    /* define variables */
+
+    std::unordered_map<std::string, int> ids;
+    for (auto const& p : parameters) {
+      int dims[3];
+      dims[0] = z_dim;
+      dims[1] = y_dim;
+      dims[2] = x_dim;
+
+      stat = nc_def_var(ncid, p.c_str(), NC_DOUBLE, 3, dims, &ids[p]);
+      check_err(stat,__LINE__,__FILE__);
+    }
+
+    /* leave define mode */
+    stat = nc_enddef (ncid);
+    check_err(stat,__LINE__,__FILE__);
+
+    /* assign variable data */
+
+    it = material.begin();
+    for (auto const& p : parameters) {
+      stat = nc_put_var_double(ncid, ids[p], it->data());
+      ++it;
+    }
 
-  stat = nc_close(ncid);
-  check_err(stat,__LINE__,__FILE__);
+    stat = nc_close(ncid);
+    check_err(stat,__LINE__,__FILE__);
+  }
 
   return 0;
 }