Merge pull request #320 from EinsteinToolkit/main

Update with changes in main

CarpetX/src/interpolate.cxx

@@ -20,6 +20,7 @@
 #include <array>
 #include <cassert>
 #include <cmath>
+#include <iostream>
 #include <limits>
 #include <map>
 #include <set>
@@ -33,9 +34,6 @@ namespace {
 
 // Interpolate a grid function at one point, dimensionally recursive
 template <typename T, int order, int centering> struct interpolator {
-  // TODO: `centering` is not used; remove it
-  // TODO: remove `dx` in favour of `x1`?
-
   static constexpr vect<bool, dim> indextype{(centering & 0b100) != 0,
                                              (centering & 0b010) != 0,
                                              (centering & 0b001) != 0};
@@ -73,12 +71,9 @@ template <typename T, int order, int centering> struct interpolator {
         for (int b = -1; b <= +1; ++b)
           for (int a = -1; a <= +1; ++a)
             if (vars.contains(j[0] + a, j[1] + b, j[2] + c))
-              // std::cerr << "  val[" << a << "," << b << "," << c
-              //           << "]=" << vars(j[0] + a, j[1] + b, j[2] + c, vi) <<
-              //           "\n";
-              std::fprintf(stderr, "val[%d,%d,%d]=%.17g\n", j[0] + a, j[1] + b,
-                           j[2] + c,
-                           double(vars(j[0] + a, j[1] + b, j[2] + c, vi)));
+              std::cerr << "  val[" << a << "," << b << "," << c
+                        << "]=" << vars(j[0] + a, j[1] + b, j[2] + c, vi)
+                        << "\n";
     }
     assert(isfinite(val));
 #endif
@@ -229,9 +224,9 @@ template <typename T, int order, int centering> struct interpolator {
     // const auto x1 = x0 + (grid.lsh - 1 - indextype) * grid.dx;
     const auto dx = grid.dx;
 
-    assert(vars.end.x - vars.begin.x == grid.lsh[0]);
-    assert(vars.end.y - vars.begin.y == grid.lsh[1]);
-    assert(vars.end.z - vars.begin.z == grid.lsh[2]);
+    assert(vars.end.x - vars.begin.x == grid.lsh[0] - indextype[0]);
+    assert(vars.end.y - vars.begin.y == grid.lsh[1] - indextype[1]);
+    assert(vars.end.z - vars.begin.z == grid.lsh[2] - indextype[2]);
 
     // We assume that the input is synchronized, i.e. that all ghost
     // zones are valid, but all outer boundaries are invalid.
@@ -445,13 +440,14 @@ extern "C" void CarpetX_Interpolate(const CCTK_POINTER_TO_CONST cctkGH_,
                               patches.data(), localsx.data(), localsy.data(),
                               localsz.data());
   } else {
-    // TODO: Don't copy
-    for (int n = 0; n < npoints; ++n) {
-      patches[n] = 0;
-      localsx[n] = globalsx[n];
-      localsy[n] = globalsy[n];
-      localsz[n] = globalsz[n];
-    }
+    for (int n = 0; n < npoints; ++n)
+      patches.at(n) = 0;
+    for (int n = 0; n < npoints; ++n)
+      localsx.at(n) = globalsx[n];
+    for (int n = 0; n < npoints; ++n)
+      localsy.at(n) = globalsy[n];
+    for (int n = 0; n < npoints; ++n)
+      localsz.at(n) = globalsz[n];
   }
 
   // Apply symmetries to coordinates
@@ -485,7 +481,7 @@ extern "C" void CarpetX_Interpolate(const CCTK_POINTER_TO_CONST cctkGH_,
   std::vector<CCTK_REAL> posy(npoints);
   std::vector<CCTK_REAL> posz(npoints);
   for (int n = 0; n < npoints; ++n) {
-    const int patch = patches[n];
+    const int patch = patches.at(n);
     const amrex::Geometry &geom = ghext->patchdata.at(patch).amrcore->Geom(0);
     const CCTK_REAL *restrict const xmin = geom.ProbLo();
     const CCTK_REAL *restrict const xmax = geom.ProbHi();
@@ -535,8 +531,9 @@ extern "C" void CarpetX_Interpolate(const CCTK_POINTER_TO_CONST cctkGH_,
 
   // Send particles to interpolation points
   for (auto &container : containers) {
-    const int patch = int(&container - containers.data());
 #ifdef CCTK_DEBUG
+    const int patch = int(&container - containers.data());
+
     std::size_t old_nparticles = 0;
     std::set<int> oldids;
     {

Loop/src/loop.hxx

@@ -187,8 +187,9 @@ public:
             // Outward boundary normal (if on outermost interior point), else 0
             const vect<int, dim> BI =
                 vect<int, dim>(I == bnd_max - 1) - vect<int, dim>(I == bnd_min);
+            const vect<bool, dim> CI1{CI, CJ, CK};
             const PointDesc p =
-                point_desc({CI, CJ, CK}, I, iter, NI, I0, BI, bnd_min, bnd_max,
+                point_desc(CI1, I, iter, NI, I0, BI, bnd_min, bnd_max,
                            loop_min, loop_max);
             f(p);
           }

TestInterpolate/param.ccl

@@ -1 +1,9 @@
 # Parameter definitions for thorn TestInterpolate
+
+BOOLEAN test_cell_interpolation "Test cell interpolation"
+{
+} "no"
+
+BOOLEAN test_vertex_interpolation "Test vertex interpolation"
+{
+} "no"

TestInterpolate/schedule.ccl

@@ -4,10 +4,22 @@ STORAGE: results_interpolate
 STORAGE: results_driver_interpolate
 STORAGE: reference
 
-SCHEDULE TestInterpolate_test_interpolation AT poststep
-{
-  LANG: C
-  OPTIONS: global
-  READS: CoordinatesX::vertex_coords(everywhere)
-  WRITES: results_interpolate, results_driver_interpolate, reference
-} "Test interpolation operators"
+if (test_cell_interpolation) {
+  SCHEDULE TestInterpolate_test_vertex_interpolation AT poststep
+  {
+    LANG: C
+    OPTIONS: global
+    READS: CoordinatesX::vertex_coords(everywhere)
+    WRITES: results_interpolate, results_driver_interpolate, reference
+  } "Test vertex interpolation operators"
+}
+
+if (test_vertex_interpolation) {
+  SCHEDULE TestInterpolate_test_cell_interpolation AT poststep
+  {
+    LANG: C
+    OPTIONS: global
+    READS: CoordinatesX::cell_coords(everywhere)
+    WRITES: results_interpolate, results_driver_interpolate, reference
+  } "Test cell interpolation operators"
+}

TestInterpolate/src/make.code.defn

@@ -1,7 +1,7 @@
 # Main make.code.defn file for thorn TestInterpolate
 
 # Source files in this directory
-SRCS = test_interpolation.cxx
+SRCS = test_cell_interpolation.cxx test_vertex_interpolation.cxx
 
 # Subdirectories containing source files
 SUBDIRS =

TestInterpolate/src/test_cell_interpolation.cxx

@@ -0,0 +1,163 @@
+#include <cctk.h>
+#include <cctk_Arguments.h>
+#include <util_Table.h>
+
+namespace TestInterpolate {
+
+#define DIM(v) ((int)(sizeof(v)/sizeof((v)[0])))
+
+extern "C" void TestInterpolate_test_cell_interpolation(CCTK_ARGUMENTS) {
+  DECLARE_CCTK_ARGUMENTS_TestInterpolate_test_cell_interpolation;
+
+  constexpr CCTK_INT N_dims = 3;
+
+  const CCTK_INT all_operations[1 + N_dims + N_dims*(N_dims-1)] = {0, 1, 2, 3, 11, 12, 13, 22, 23, 33};
+
+  const CCTK_INT all_varinds[N_dims] = {
+      CCTK_VarIndex("CoordinatesX::ccoordx"),
+      CCTK_VarIndex("CoordinatesX::ccoordy"),
+      CCTK_VarIndex("CoordinatesX::ccoordz"),
+  };
+
+  constexpr int nvars = DIM(all_varinds) * DIM(all_operations);
+
+  CCTK_INT varinds[DIM(all_operations)][DIM(all_varinds)];
+  CCTK_INT operations[DIM(all_operations)][DIM(all_varinds)];
+  for (int op = 0 ; op < DIM(all_operations) ; op++) {
+    for (int var = 0 ; var < DIM(all_varinds) ; var++) {
+      varinds[op][var] = all_varinds[var];
+      operations[op][var] = all_operations[op];
+    }
+  }
+
+  constexpr int npoints = 4;
+  const CCTK_REAL coords[N_dims][npoints] = {
+      {0.1, 0.2, 0.3, 0.1},
+      {0.1, 0.2, 0.3, 0.2},
+      {0.1, 0.2, 0.3, 0.3},
+  };
+
+  CCTK_REAL *resultptrs[DIM(all_operations)][DIM(all_varinds)] = {
+    interpolate_x_o1_op0,
+    interpolate_y_o1_op0,
+    interpolate_z_o1_op0,
+
+    interpolate_x_o1_op1, interpolate_y_o1_op1, interpolate_z_o1_op1,
+    interpolate_x_o1_op2, interpolate_y_o1_op2, interpolate_z_o1_op2,
+    interpolate_x_o1_op3, interpolate_y_o1_op3, interpolate_z_o1_op3,
+
+    interpolate_x_o1_op11, interpolate_y_o1_op11, interpolate_z_o1_op11,
+    interpolate_x_o1_op12, interpolate_y_o1_op12, interpolate_z_o1_op12,
+    interpolate_x_o1_op13, interpolate_y_o1_op13, interpolate_z_o1_op13,
+    interpolate_x_o1_op22, interpolate_y_o1_op22, interpolate_z_o1_op22,
+    interpolate_x_o1_op23, interpolate_y_o1_op23, interpolate_z_o1_op23,
+    interpolate_x_o1_op33, interpolate_y_o1_op33, interpolate_z_o1_op33,
+  };
+  CCTK_REAL *resultptrs2[DIM(all_operations)][DIM(all_varinds)] = {
+    driver_interpolate_x_o1_op0,
+    driver_interpolate_y_o1_op0,
+    driver_interpolate_z_o1_op0,
+
+    driver_interpolate_x_o1_op1, driver_interpolate_y_o1_op1, driver_interpolate_z_o1_op1,
+    driver_interpolate_x_o1_op2, driver_interpolate_y_o1_op2, driver_interpolate_z_o1_op2,
+    driver_interpolate_x_o1_op3, driver_interpolate_y_o1_op3, driver_interpolate_z_o1_op3,
+
+    driver_interpolate_x_o1_op11, driver_interpolate_y_o1_op11, driver_interpolate_z_o1_op11,
+    driver_interpolate_x_o1_op12, driver_interpolate_y_o1_op12, driver_interpolate_z_o1_op12,
+    driver_interpolate_x_o1_op13, driver_interpolate_y_o1_op13, driver_interpolate_z_o1_op13,
+    driver_interpolate_x_o1_op22, driver_interpolate_y_o1_op22, driver_interpolate_z_o1_op22,
+    driver_interpolate_x_o1_op23, driver_interpolate_y_o1_op23, driver_interpolate_z_o1_op23,
+    driver_interpolate_x_o1_op33, driver_interpolate_y_o1_op33, driver_interpolate_z_o1_op33,
+  };
+  CCTK_REAL *referenceptrs[DIM(all_operations)][DIM(all_varinds)] = {
+    reference_x_o1_op0,
+    reference_y_o1_op0,
+    reference_z_o1_op0,
+
+    reference_x_o1_op1, reference_y_o1_op1, reference_z_o1_op1,
+    reference_x_o1_op2, reference_y_o1_op2, reference_z_o1_op2,
+    reference_x_o1_op3, reference_y_o1_op3, reference_z_o1_op3,
+
+    reference_x_o1_op11, reference_y_o1_op11, reference_z_o1_op11,
+    reference_x_o1_op12, reference_y_o1_op12, reference_z_o1_op12,
+    reference_x_o1_op13, reference_y_o1_op13, reference_z_o1_op13,
+    reference_x_o1_op22, reference_y_o1_op22, reference_z_o1_op22,
+    reference_x_o1_op23, reference_y_o1_op23, reference_z_o1_op23,
+    reference_x_o1_op33, reference_y_o1_op33, reference_z_o1_op33,
+  };
+
+  Interpolate(cctkGH, npoints, coords[0], coords[1], coords[2], nvars,
+      (CCTK_INT const * const)varinds, (CCTK_INT const * const)operations, 1,
+      (CCTK_REAL **)resultptrs);
+
+  const void* interp_coords[N_dims] = {
+    coords[0], coords[1], coords[2]
+  };
+  void *const *output_array = (void *const *)resultptrs2;
+
+  /* DriverInterpolate arguments that aren't currently used */
+  CCTK_INT const coord_system_handle = 0;
+  CCTK_INT const interp_coords_type_code = 0;
+  CCTK_INT const output_array_type_codes[1] = {0};
+  int interp_handle = CCTK_InterpHandle("CarpetX");
+
+  if(interp_handle < 0)
+    CCTK_VERROR("Could not obtain inteprolator handle for built-in 'CarpetX' interpolator: %d", interp_handle);
+
+  /* Table generation */
+  int operands[DIM(all_operations)][DIM(all_varinds)];
+  for (int op = 0 ; op < DIM(all_operations) ; op++) {
+    for (int var = 0 ; var < DIM(all_varinds) ; var++) {
+      operands[op][var] = var;
+    }
+  }
+
+  int ierr;
+  int param_table_handle = Util_TableCreate(UTIL_TABLE_FLAGS_DEFAULT);
+  if (param_table_handle < 0)
+    CCTK_VERROR("Can't create parameter table: %d", param_table_handle);
+  if ((ierr = Util_TableSetInt(param_table_handle, 1, "order")) < 0)
+    CCTK_VERROR("Can't set order in parameter table: %d", ierr);
+  if ((ierr = Util_TableSetIntArray(param_table_handle, nvars, (int const*const)operands,
+                            "operand_indices")) < 0)
+    CCTK_VERROR("Can't set operand_indices array in parameter table: %d", ierr);
+  if ((ierr = Util_TableSetIntArray(param_table_handle, nvars, (int const*const)operations,
+                            "operation_codes")) < 0)
+    CCTK_VERROR("Can't set operation_codes array in parameter table: %d", ierr);
+
+  DriverInterpolate(
+      cctkGH, N_dims, interp_handle, param_table_handle, coord_system_handle,
+      npoints, interp_coords_type_code, interp_coords, nvars, (int const*const)varinds,
+      nvars, output_array_type_codes, output_array);
+
+  for (int op = 0 ; op < DIM(all_operations) ; op++) {
+    for (int var = 0 ; var < DIM(all_varinds) ; var++) {
+      for (int n = 0; n < npoints; ++n) {
+        const int operand = operations[op][var];
+        switch (operand) {
+        case 0:
+          referenceptrs[op][var][n] = coords[var][n];
+          break;
+        case 1:
+        case 2:
+        case 3:
+          referenceptrs[op][var][n] = (operand - 1 == var ? 1. : 0.);
+          break;
+        case 11:
+        case 12:
+        case 13:
+        case 22:
+        case 23:
+        case 33:
+          referenceptrs[op][var][n] = 0.;
+          break;
+        default:
+          CCTK_VERROR("Unsupported derivative operator %d", operand);
+          break;
+        }
+      }
+    }
+  }
+}
+
+} // namespace TestInterpolate

TestInterpolate/src/test_interpolation.cxx → TestInterpolate/src/test_vertex_interpolation.cxx

@@ -6,8 +6,8 @@ namespace TestInterpolate {
 
 #define DIM(v) ((int)(sizeof(v)/sizeof((v)[0])))
 
-extern "C" void TestInterpolate_test_interpolation(CCTK_ARGUMENTS) {
-  DECLARE_CCTK_ARGUMENTS_TestInterpolate_test_interpolation;
+extern "C" void TestInterpolate_test_vertex_interpolation(CCTK_ARGUMENTS) {
+  DECLARE_CCTK_ARGUMENTS_TestInterpolate_test_vertex_interpolation;
 
   constexpr CCTK_INT N_dims = 3;
 

TestInterpolate/test/interpolate.par → TestInterpolate/test/cell_interpolate.par

@@ -35,6 +35,8 @@ Driver::regrid_error_threshold = 5.0
 
 CarpetX::prolongation_type = "ddf"
 
+TestInterpolate::test_cell_interpolation = yes
+
 IO::out_dir = $parfile
 IO::out_fileinfo = "axis labels"
 IO::parfile_write = "no"

TestInterpolate/test/vertex_interpolate.par

@@ -0,0 +1,50 @@
+ActiveThorns = "
+    CarpetX
+    CoordinatesX
+    IOUtil
+    TestInterpolate
+"
+
+$nlevels = 3
+$ncells = 32
+
+Cactus::cctk_show_schedule = no
+Cactus::presync_mode = "mixed-error"
+
+Cactus::terminate = "time"
+Cactus::cctk_final_time = 0
+
+CarpetX::verbose = no
+CarpetX::poison_undefined_values = yes
+
+CarpetX::xmin = -1.0
+CarpetX::ymin = -1.0
+CarpetX::zmin = -1.0
+
+CarpetX::xmax = +1.0
+CarpetX::ymax = +1.0
+CarpetX::zmax = +1.0
+
+CarpetX::ncells_x = $ncells
+CarpetX::ncells_y = $ncells
+CarpetX::ncells_z = $ncells
+
+CarpetX::max_num_levels = $nlevels
+CarpetX::regrid_every = 16
+Driver::regrid_error_threshold = 5.0
+
+CarpetX::prolongation_type = "ddf"
+
+TestInterpolate::test_vertex_interpolation = yes
+
+IO::out_dir = $parfile
+IO::out_fileinfo = "axis labels"
+IO::parfile_write = "no"
+IO::out_every = 1
+
+CarpetX::out_metadata = no
+CarpetX::out_tsv_vars = "
+    TestInterpolate::reference
+    TestInterpolate::results_driver_interpolate
+    TestInterpolate::results_interpolate
+"

TestInterpolate/test/vertex_interpolate/testinterpolate-reference.it000000.tsv

@@ -0,0 +1,2 @@
+# 1:iteration	2:time	3:reference_x_o1_op0[0]	4:reference_x_o1_op0[1]	5:reference_x_o1_op0[2]	6:reference_x_o1_op0[3]	7:reference_y_o1_op0[0]	8:reference_y_o1_op0[1]	9:reference_y_o1_op0[2]	10:reference_y_o1_op0[3]	11:reference_z_o1_op0[0]	12:reference_z_o1_op0[1]	13:reference_z_o1_op0[2]	14:reference_z_o1_op0[3]	15:reference_x_o1_op1[0]	16:reference_x_o1_op1[1]	17:reference_x_o1_op1[2]	18:reference_x_o1_op1[3]	19:reference_y_o1_op1[0]	20:reference_y_o1_op1[1]	21:reference_y_o1_op1[2]	22:reference_y_o1_op1[3]	23:reference_z_o1_op1[0]	24:reference_z_o1_op1[1]	25:reference_z_o1_op1[2]	26:reference_z_o1_op1[3]	27:reference_x_o1_op2[0]	28:reference_x_o1_op2[1]	29:reference_x_o1_op2[2]	30:reference_x_o1_op2[3]	31:reference_y_o1_op2[0]	32:reference_y_o1_op2[1]	33:reference_y_o1_op2[2]	34:reference_y_o1_op2[3]	35:reference_z_o1_op2[0]	36:reference_z_o1_op2[1]	37:reference_z_o1_op2[2]	38:reference_z_o1_op2[3]	39:reference_x_o1_op3[0]	40:reference_x_o1_op3[1]	41:reference_x_o1_op3[2]	42:reference_x_o1_op3[3]	43:reference_y_o1_op3[0]	44:reference_y_o1_op3[1]	45:reference_y_o1_op3[2]	46:reference_y_o1_op3[3]	47:reference_z_o1_op3[0]	48:reference_z_o1_op3[1]	49:reference_z_o1_op3[2]	50:reference_z_o1_op3[3]	51:reference_x_o1_op11[0]	52:reference_x_o1_op11[1]	53:reference_x_o1_op11[2]	54:reference_x_o1_op11[3]	55:reference_y_o1_op11[0]	56:reference_y_o1_op11[1]	57:reference_y_o1_op11[2]	58:reference_y_o1_op11[3]	59:reference_z_o1_op11[0]	60:reference_z_o1_op11[1]	61:reference_z_o1_op11[2]	62:reference_z_o1_op11[3]	63:reference_x_o1_op12[0]	64:reference_x_o1_op12[1]	65:reference_x_o1_op12[2]	66:reference_x_o1_op12[3]	67:reference_y_o1_op12[0]	68:reference_y_o1_op12[1]	69:reference_y_o1_op12[2]	70:reference_y_o1_op12[3]	71:reference_z_o1_op12[0]	72:reference_z_o1_op12[1]	73:reference_z_o1_op12[2]	74:reference_z_o1_op12[3]	75:reference_x_o1_op13[0]	76:reference_x_o1_op13[1]	77:reference_x_o1_op13[2]	78:reference_x_o1_op13[3]	79:reference_y_o1_op13[0]	80:reference_y_o1_op13[1]	81:reference_y_o1_op13[2]	82:reference_y_o1_op13[3]	83:reference_z_o1_op13[0]	84:reference_z_o1_op13[1]	85:reference_z_o1_op13[2]	86:reference_z_o1_op13[3]	87:reference_x_o1_op22[0]	88:reference_x_o1_op22[1]	89:reference_x_o1_op22[2]	90:reference_x_o1_op22[3]	91:reference_y_o1_op22[0]	92:reference_y_o1_op22[1]	93:reference_y_o1_op22[2]	94:reference_y_o1_op22[3]	95:reference_z_o1_op22[0]	96:reference_z_o1_op22[1]	97:reference_z_o1_op22[2]	98:reference_z_o1_op22[3]	99:reference_x_o1_op23[0]	100:reference_x_o1_op23[1]	101:reference_x_o1_op23[2]	102:reference_x_o1_op23[3]	103:reference_y_o1_op23[0]	104:reference_y_o1_op23[1]	105:reference_y_o1_op23[2]	106:reference_y_o1_op23[3]	107:reference_z_o1_op23[0]	108:reference_z_o1_op23[1]	109:reference_z_o1_op23[2]	110:reference_z_o1_op23[3]	111:reference_x_o1_op33[0]	112:reference_x_o1_op33[1]	113:reference_x_o1_op33[2]	114:reference_x_o1_op33[3]	115:reference_y_o1_op33[0]	116:reference_y_o1_op33[1]	117:reference_y_o1_op33[2]	118:reference_y_o1_op33[3]	119:reference_z_o1_op33[0]	120:reference_z_o1_op33[1]	121:reference_z_o1_op33[2]	122:reference_z_o1_op33[3]
+0	0.0000000000000000e+00	1.0000000000000001e-01	2.0000000000000001e-01	2.9999999999999999e-01	1.0000000000000001e-01	1.0000000000000001e-01	2.0000000000000001e-01	2.9999999999999999e-01	2.0000000000000001e-01	1.0000000000000001e-01	2.0000000000000001e-01	2.9999999999999999e-01	2.9999999999999999e-01	1.0000000000000000e+00	1.0000000000000000e+00	1.0000000000000000e+00	1.0000000000000000e+00	0.0000000000000000e+00	0.0000000000000000e+00	0.0000000000000000e+00	0.0000000000000000e+00	0.0000000000000000e+00	0.0000000000000000e+00	0.0000000000000000e+00	0.0000000000000000e+00	0.0000000000000000e+00	0.0000000000000000e+00	0.0000000000000000e+00	0.0000000000000000e+00	1.0000000000000000e+00	1.0000000000000000e+00	1.0000000000000000e+00	1.0000000000000000e+00	0.0000000000000000e+00	0.0000000000000000e+00	0.0000000000000000e+00	0.0000000000000000e+00	0.0000000000000000e+00	0.0000000000000000e+00	0.0000000000000000e+00	0.0000000000000000e+00	0.0000000000000000e+00	0.0000000000000000e+00	0.0000000000000000e+00	0.0000000000000000e+00	1.0000000000000000e+00	1.0000000000000000e+00	1.0000000000000000e+00	1.0000000000000000e+00	0.0000000000000000e+00	0.0000000000000000e+00	0.0000000000000000e+00	0.0000000000000000e+00	0.0000000000000000e+00	0.0000000000000000e+00	0.0000000000000000e+00	0.0000000000000000e+00	0.0000000000000000e+00	0.0000000000000000e+00	0.0000000000000000e+00	0.0000000000000000e+00	0.0000000000000000e+00	0.0000000000000000e+00	0.0000000000000000e+00	0.0000000000000000e+00	0.0000000000000000e+00	0.0000000000000000e+00	0.0000000000000000e+00	0.0000000000000000e+00	0.0000000000000000e+00	0.0000000000000000e+00	0.0000000000000000e+00	0.0000000000000000e+00	0.0000000000000000e+00	0.0000000000000000e+00	0.0000000000000000e+00	0.0000000000000000e+00	0.0000000000000000e+00	0.0000000000000000e+00	0.0000000000000000e+00	0.0000000000000000e+00	0.0000000000000000e+00	0.0000000000000000e+00	0.0000000000000000e+00	0.0000000000000000e+00	0.0000000000000000e+00	0.0000000000000000e+00	0.0000000000000000e+00	0.0000000000000000e+00	0.0000000000000000e+00	0.0000000000000000e+00	0.0000000000000000e+00	0.0000000000000000e+00	0.0000000000000000e+00	0.0000000000000000e+00	0.0000000000000000e+00	0.0000000000000000e+00	0.0000000000000000e+00	0.0000000000000000e+00	0.0000000000000000e+00	0.0000000000000000e+00	0.0000000000000000e+00	0.0000000000000000e+00	0.0000000000000000e+00	0.0000000000000000e+00	0.0000000000000000e+00	0.0000000000000000e+00	0.0000000000000000e+00	0.0000000000000000e+00	0.0000000000000000e+00	0.0000000000000000e+00	0.0000000000000000e+00	0.0000000000000000e+00	0.0000000000000000e+00	0.0000000000000000e+00	0.0000000000000000e+00	0.0000000000000000e+00	0.0000000000000000e+00	0.0000000000000000e+00	0.0000000000000000e+00	0.0000000000000000e+00