added min_max normalization (temporary fix)

test/CMakeLists.txt

@@ -131,6 +131,13 @@ else()
 message(WARNING "Catch2 not found. Disabling Unit Tests")
 endif()
 
+if (PCMS_ENABLE_MESHFIELDS)
+  add_executable(thwaites_adapt thwaites_adapt.cpp)
+  target_link_libraries(thwaites_adapt PUBLIC meshfields::meshfields)
+  target_link_libraries(thwaites_adapt PUBLIC pcms::core pcms_interpolator)
+  target_include_directories(thwaites_adapt PUBLIC ${CMAKE_CURRENT_SOURCE_DIR})
+endif()
+
 if (PCMS_ENABLE_C)
     find_package(Kokkos REQUIRED)
     if (PCMS_ENABLE_XGC)

test/test_spr_meshfields.cpp

@@ -11,6 +11,7 @@
 #include <catch2/matchers/catch_matchers_floating_point.hpp>
 #include <vector>
 #include <iostream>
+#include <iomanip>
 
 using namespace std;
 using namespace Omega_h;
@@ -29,6 +30,31 @@ void print_patches(Omega_h::Graph& patches) {
     }
     std::cout << "\n";
   }
+}
+
+void print_patch_vals(Omega_h::Graph& patches_d, Reals vtxCoords_d,
+    Reals elmSrcVals_d, Reals elmCentroids, size_t vtx) {
+  const auto meshDim = 2;
+  const auto offsets = HostRead(patches_d.a2ab);
+  const auto values = HostRead(patches_d.ab2b);
+  const auto vtxCoords = HostRead(vtxCoords_d);
+  const auto elmSrcVals = HostRead(elmSrcVals_d);
+  std::cout << std::setprecision (15);
+  std::cout << "num patches " << patches_d.nnodes() << "\n";
+  for(int patch=0; patch<patches_d.nnodes(); patch++) {
+    if(patch == vtx) {
+      std::cout << "vtxCoords[" << patch << "] " << vtxCoords[patch*meshDim] << " " << vtxCoords[patch*meshDim+1] << "\n";
+      std::cout << "<elementIdx> <centroid x> <centroid y> <source_value>\n";
+      for (auto valIdx = offsets[patch]; valIdx < offsets[patch + 1]; ++valIdx) {
+        auto patchElm = values[valIdx];
+        std::cout <<  patchElm << " " 
+                  << elmCentroids[patchElm*meshDim] << " "
+                  << elmCentroids[patchElm*meshDim+1] << " "
+                  << elmSrcVals[patchElm] << "\n";
+      }
+      std::cout << "\n";
+    }
+  }
 
 }
 
@@ -88,7 +114,13 @@ void test(Mesh& mesh, Omega_h::Graph& patches, int degree,
   REQUIRE(m == n);
 
   for (size_t i = 0; i < m; ++i) {
-    std::cout << "vtx " << i << "\n";
+    if( delta_abs[i] > tolerance ) {
+      std::cout << std::setprecision (15);
+      std::cout << "exact_target_values[" << i << "] " << exact_target_values[i] << "\n";
+      std::cout << "approx_target_values[" << i << "] " << approx_target_values[i] << "\n";
+      std::cout << "abs(exact-approx) " << delta_abs[i] << "\n";
+      print_patch_vals(patches, target_coordinates, source_values, source_coordinates, i);
+    }
     CHECK_THAT(
         host_exact_target_values[i],
         Catch::Matchers::WithinAbs(host_approx_target_values[i], tolerance));
@@ -106,7 +138,10 @@ TEST_CASE("meshfields_spr_test")
   auto rank = lib.world()->rank();
   const auto boxSize = 1.0;
   const auto nElms = 6;
-  auto mesh = build_box(world, OMEGA_H_SIMPLEX, boxSize, boxSize, 0, nElms, nElms, 0, false);
+  const auto thwaitesMeshFile = "/lore/smithc11/projects/landice/thwaites_basal/thwaites_basalClass.osh";
+  Omega_h::Mesh mesh(&lib);
+  Omega_h::binary::read(thwaitesMeshFile, lib.world(), &mesh);
+  //auto mesh = build_box(world, OMEGA_H_SIMPLEX, boxSize, boxSize, 0, nElms, nElms, 0, false);
   std::cout << "mesh: elms " << mesh.nelems() << " verts " << mesh.nverts() << "\n";
 
   const auto dim = mesh.dim();
@@ -166,7 +201,7 @@ TEST_CASE("meshfields_spr_test")
         if(interp_degree == 3) minPatchSize = 10; //why so large?
         std::cerr << "minPatchSize " << minPatchSize << "\n";
         auto patches = mesh.get_vtx_patches(minPatchSize);
-        print_patches(patches);
+        //print_patches(patches);
 
         Write<Real> source_values(nfaces, 0, "exact target values");
 

test/thwaites_adapt.cpp

@@ -0,0 +1,159 @@
+#include <iostream>
+
+#include "Omega_h_adapt.hpp"
+#include "Omega_h_array_ops.hpp"
+#include "Omega_h_build.hpp"
+#include "Omega_h_class.hpp"
+#include "Omega_h_compare.hpp"
+#include "Omega_h_for.hpp"
+#include "Omega_h_shape.hpp"
+#include "Omega_h_timer.hpp"
+#include "Omega_h_recover.hpp" //project_by_fit
+#include <Omega_h_file.hpp> //Omega_h::binary
+#include <Omega_h_atomics.hpp> //Omega_h::atomic_fetch_add
+#include <sstream> //ostringstream
+#include <iomanip> //precision
+#include <Omega_h_dbg.hpp>
+#include <Omega_h_for.hpp>
+
+#include "thwaites_errorEstimator.hpp"
+
+//detect floating point exceptions
+#include <fenv.h>
+
+using ExecutionSpace = Kokkos::DefaultExecutionSpace;
+using MemorySpace = Kokkos::DefaultExecutionSpace::memory_space;
+
+using namespace Omega_h;
+
+void setupFieldTransfer(AdaptOpts& opts) {
+  opts.xfer_opts.type_map["solution_1"] = OMEGA_H_LINEAR_INTERP;
+  opts.xfer_opts.type_map["solution_2"] = OMEGA_H_LINEAR_INTERP;
+  opts.xfer_opts.type_map["ice_thickness"] = OMEGA_H_LINEAR_INTERP;
+  opts.xfer_opts.type_map["surface_height"] = OMEGA_H_LINEAR_INTERP;
+  opts.xfer_opts.type_map["bed_topography"] = OMEGA_H_LINEAR_INTERP;
+  opts.xfer_opts.type_map["mu_log"] = OMEGA_H_LINEAR_INTERP;
+  opts.xfer_opts.type_map["prescribed_velocity_1"] = OMEGA_H_LINEAR_INTERP;
+  opts.xfer_opts.type_map["prescribed_velocity_2"] = OMEGA_H_LINEAR_INTERP;
+  opts.xfer_opts.type_map["observed_surface_velocity_rms"] = OMEGA_H_LINEAR_INTERP;
+  opts.xfer_opts.type_map["observed_surface_velocity_1"] = OMEGA_H_LINEAR_INTERP;
+  opts.xfer_opts.type_map["observed_surface_velocity_2"] = OMEGA_H_LINEAR_INTERP;
+  opts.xfer_opts.type_map["stiffening_factor_log"] = OMEGA_H_LINEAR_INTERP;
+  const int numLayers = 11;
+  for(int i=1; i<=numLayers; i++) {
+    std::stringstream ss;
+    ss << "temperature_" << std::setfill('0') << std::setw(2) << i;
+    opts.xfer_opts.type_map[ss.str()] = OMEGA_H_LINEAR_INTERP;
+  }
+}
+
+/**
+ * retrieve the effective strain rate from the mesh
+ *
+ * despite the name being 'effective_stress' it is the effective strain:
+ * the frobenius norm of  the strain tensor
+  */
+Reals getEffectiveStrainRate(Mesh& mesh) {
+  return mesh.get_array<Real>(2, "effective_stress");
+}
+
+Reals recoverLinearStrain(Mesh& mesh, Reals effectiveStrain) {
+  return project_by_fit(&mesh, effectiveStrain);
+}
+
+Reals computeError(Mesh& mesh, Reals effectiveStrain, Reals recoveredStrain) {
+  return Reals(mesh.nelems());
+}
+
+template <typename ShapeField>
+void setFieldAtVertices(Omega_h::Mesh &mesh, Reals recoveredStrain, ShapeField field) {
+  const auto MeshDim = mesh.dim();
+  auto setFieldAtVertices = KOKKOS_LAMBDA(const int &vtx) {
+    field(0, 0, vtx, MeshField::Vertex) = recoveredStrain[vtx];
+  };
+  MeshField::parallel_for(ExecutionSpace(), {0}, {mesh.nverts()},
+                          setFieldAtVertices, "setFieldAtVertices");
+}
+
+void printTriCount(Mesh* mesh) {
+  const auto nTri = mesh->nglobal_ents(2);
+  if (!mesh->comm()->rank())
+    std::cout << "nTri: " << nTri << "\n";
+}
+
+int main(int argc, char** argv) {
+  feenableexcept(FE_ALL_EXCEPT & ~FE_INEXACT);  // Enable all floating point exceptions but FE_INEXACT
+  auto lib = Library(&argc, &argv);
+  if( argc != 6 ) {
+    fprintf(stderr, "Usage: %s inputMesh.osh outputMeshPrefix adaptRatio min_size max_size\n", argv[0]);
+    exit(EXIT_FAILURE);
+  }
+  auto world = lib.world();
+  Omega_h::Mesh mesh(&lib);
+  Omega_h::binary::read(argv[1], world, &mesh);
+  const auto prefix = std::string(argv[2]);
+  //adaptRatio = 0.1 is used in scorec/core:cws/sprThwaites test/spr_test.cc
+  const MeshField::Real adaptRatio = std::stof(argv[3]);
+  const MeshField::Real min_size = std::stof(argv[4]); //the smallest edge length in initial mesh is 1
+  const MeshField::Real max_size = std::stof(argv[5]); //a value of 8 or 16 roughly maintains the boundary shape
+                                                       //in the downstream parts of the domain
+                                                       //where there is significant coarsening
+  std::cout << "input mesh: " << argv[1] << " outputMeshPrefix: "
+            << prefix << " adaptRatio: " << adaptRatio
+            << " max_size: " << max_size
+            << " min_size: " << min_size << "\n";
+  const auto outname = prefix + "_adaptRatio_" + std::string(argv[3]) +
+                       "_maxSz_" + std::to_string(max_size) +
+                       "_minSz_" + std::to_string(min_size);
+
+  auto effectiveStrain = getEffectiveStrainRate(mesh);
+  auto recoveredStrain = recoverLinearStrain(mesh,effectiveStrain);
+  mesh.add_tag<Real>(VERT, "recoveredStrain", 1, recoveredStrain);
+  auto error = computeError(mesh,effectiveStrain,recoveredStrain);
+
+  MeshField::OmegahMeshField<ExecutionSpace, MeshField::KokkosController> omf(
+        mesh);
+
+  const auto MeshDim = 2;
+  const auto ShapeOrder = 1;
+  auto recoveredStrainField = omf.CreateLagrangeField<Real, ShapeOrder, MeshDim>();
+  setFieldAtVertices(mesh, recoveredStrain, recoveredStrainField);
+
+  auto coordField = omf.getCoordField();
+  const auto [shp, map] = MeshField::Omegah::getTriangleElement<ShapeOrder>(mesh);
+  MeshField::FieldElement coordFe(mesh.nelems(), coordField, shp, map);
+
+  auto estimation = Estimation(mesh, effectiveStrain,
+      recoveredStrainField, adaptRatio);
+
+  const auto tgtLength = getSprSizeField(estimation, omf, coordFe);
+  Omega_h::Write<MeshField::Real> tgtLength_oh(tgtLength);
+  mesh.add_tag<Real>(VERT, "tgtLength", 1, tgtLength_oh);
+
+  { //write vtk
+  const std::string vtkFileName = "beforeAdapt" + outname + ".vtk";
+  Omega_h::vtk::write_parallel(vtkFileName, &mesh, 2);
+  const std::string vtkFileName_edges = "beforeAdapt" + outname + "_edges.vtk";
+  Omega_h::vtk::write_parallel(vtkFileName_edges, &mesh, 1);
+  }
+
+  //adapt
+  auto opts = Omega_h::AdaptOpts(&mesh);
+  setupFieldTransfer(opts);
+  printTriCount(&mesh);
+
+  auto verbose = true;
+  const auto isos = Omega_h::isos_from_lengths(tgtLength_oh);
+  auto metric = clamp_metrics(mesh.nverts(), isos, min_size, max_size);
+  Omega_h::grade_fix_adapt(&mesh, opts, metric, verbose);
+
+  { //write vtk and osh for adapted mesh
+  const std::string outfilename = "afterAdapt" + outname;
+  Omega_h::vtk::write_parallel(outfilename + ".vtk", &mesh, 2);
+  Omega_h::binary::write(outfilename + ".osh", &mesh);
+  std::cout << "wrote adapted mesh: " << outfilename + ".osh" << "\n";
+  }
+
+
+  return 0;
+}