diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index 5bde113b..3c9785c4 100644 --- a/test/CMakeLists.txt +++ b/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) diff --git a/test/test_spr_meshfields.cpp b/test/test_spr_meshfields.cpp index 4d61b3c0..6fed65bc 100644 --- a/test/test_spr_meshfields.cpp +++ b/test/test_spr_meshfields.cpp @@ -11,6 +11,7 @@ #include #include #include +#include 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 \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 source_values(nfaces, 0, "exact target values"); diff --git a/test/thwaites_adapt.cpp b/test/thwaites_adapt.cpp new file mode 100644 index 00000000..71ecb5fe --- /dev/null +++ b/test/thwaites_adapt.cpp @@ -0,0 +1,159 @@ +#include + +#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::binary +#include //Omega_h::atomic_fetch_add +#include //ostringstream +#include //precision +#include +#include + +#include "thwaites_errorEstimator.hpp" + +//detect floating point exceptions +#include + +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(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 +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(VERT, "recoveredStrain", 1, recoveredStrain); + auto error = computeError(mesh,effectiveStrain,recoveredStrain); + + MeshField::OmegahMeshField omf( + mesh); + + const auto MeshDim = 2; + const auto ShapeOrder = 1; + auto recoveredStrainField = omf.CreateLagrangeField(); + setFieldAtVertices(mesh, recoveredStrain, recoveredStrainField); + + auto coordField = omf.getCoordField(); + const auto [shp, map] = MeshField::Omegah::getTriangleElement(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 tgtLength_oh(tgtLength); + mesh.add_tag(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; +} diff --git a/test/thwaites_errorEstimator.hpp b/test/thwaites_errorEstimator.hpp new file mode 100644 index 00000000..dd8e090c --- /dev/null +++ b/test/thwaites_errorEstimator.hpp @@ -0,0 +1,245 @@ +#ifndef THWAITES_ERROR_ESTIMATOR_H +#define THWAITES_ERROR_ESTIMATOR_H +#include +#include +#include +#include +#include +#include +#include +#include + + +/* useful for initializing values to quickly + detect "uninitialized" value bugs */ +static double getNaN() +{ + return std::numeric_limits::quiet_NaN(); +} + +/* common base for Scalar Integrator. */ +class SInt : public MeshField::Integrator { + public: + SInt(int order): + MeshField::Integrator(order),r(0) + {} + void reset() {r=0;} + MeshField::Real r; +}; + +template +class Estimation { + private: + Estimation() {} + public: + Estimation(Omega_h::Mesh& mesh_in, Omega_h::Reals eps_in, ShapeField eps_star_in, + MeshField::Real tolerance_in) + : mesh(mesh_in), eps(eps_in), eps_star(eps_star_in), + tolerance(tolerance_in), size_factor(getNaN()) + { + /* note that getOrder being used to convey this meaning + is a bit of a hack, but looks decent if you don't + try to define the FieldShape API too rigorously */ + integration_order = 1; //FIXME get from eps_star + /* so far recovery order is directly tied to the + mesh's coordinate field order, coordinate this + with field recovery code */ + recovered_order = 1; //FIXME + } + Omega_h::Mesh& mesh; + /* the maximum polynomial order that can be + integrated exactly using the input integration points. + not necessarily equal to recovered_order, sometimes + users give more points than strictly necessary */ + int integration_order; + /* the polynomial order of the recovered field */ + int recovered_order; + /* the input field consisting of values of a key + quantity at integration points (stress or strain, for example). + this field is related to integration_order */ + Omega_h::Reals eps; + /* the recovered field, consisting of nodal values + of the quantity of interest. + this uses a Lagrange basis of recovered_order */ + ShapeField& eps_star; + /* the acceptable margin of error, expressed as a factor + greater than zero. + setting this equal to zero would request zero element + size everywhere, i.e. infinite refinement. + increasing it scales up the desired element size + throughout. + basically, this is a (nonlinear) scaling factor on the resulting + size field */ + Omega_h::Real tolerance; + /* the uniform linear scaling factor derived from the tolerance + and integrals of the recovered field over the mesh. + desired element size = current size * current error * size_factor */ + Omega_h::Real size_factor; + /* a temporary field storing desired sizes at elements */ + Kokkos::View element_size; + /* the resulting size field, recovered from the element_size field + (using a local average recovery method much weaker than SPR) */ + Omega_h::Reals size; +}; + +/* computes $\|f\|^2$ */ +template +class SelfProduct : public SInt +{ + public: + SelfProduct(EstimationT& estimation_in, OmegahMeshField& omf_in): + SInt(estimation_in.integration_order), + estimation(estimation_in), + omf(omf_in) + { + } + void atPoints(Kokkos::View p, + Kokkos::View w, + Kokkos::View dV) + { + + std::cerr << "SelfProduct::atPoints(...)\n"; + const size_t numPtsPerElem = p.extent(0)/estimation.mesh.nelems(); + auto eps_star_atPts = omf.triangleLocalPointEval(p,numPtsPerElem, + estimation.eps_star); + + r = 0; + Kokkos::parallel_reduce( + "eval", estimation.mesh.nelems(), + KOKKOS_LAMBDA(const int &elm, MeshField::Real &r_local) { + const auto first = elm * numPtsPerElem; + const auto last = first + numPtsPerElem; + for (auto pt = first; pt < last; pt++) { + const auto vPt = eps_star_atPts(pt,0); + const auto wPt = w(pt); + const auto dVPt = dV(pt); + r_local += (vPt * vPt) * wPt * dVPt; + } + }, + r); + + } + private: + EstimationT& estimation; + OmegahMeshField& omf; +}; + +/* computes the integral over the element of the + sum of the squared differences between the + original and recovered fields */ +template +class Error : public SInt +{ + public: + Error(EstimationT& estimation_in, OmegahMeshField& omf_in): + SInt(estimation_in.integration_order), + estimation(estimation_in), + omf(omf_in), + errorNorm("errorNorm", estimation_in.mesh.nelems()) + { + } + void atPoints(Kokkos::View p, + Kokkos::View w, + Kokkos::View dV) + { + std::cerr << "SelfProduct::atPoints(...)\n"; + const size_t numPtsPerElem = p.extent(0)/estimation.mesh.nelems(); + //FIXME eps isn't a ShapeField so we can't call + // omf.triangleLocalPointEval. For now, just get + // the value from the element and assert that there is + // one integration point per element. + assert(numPtsPerElem == 1); + auto eps_star_atPts = omf.triangleLocalPointEval(p,numPtsPerElem, + estimation.eps_star); + double meshDim = estimation.mesh.dim(); + double orderP = estimation.recovered_order; + + const auto& eps = estimation.eps; + r = 0; + Kokkos::parallel_reduce( + "eval", estimation.mesh.nelems(), + KOKKOS_LAMBDA(const int &elm, MeshField::Real &r_local) { + const auto first = elm * numPtsPerElem; + const auto last = first + numPtsPerElem; + const auto eps_elm = eps[elm]; + MeshField::Real sum = 0; + for (auto pt = first; pt < last; pt++) { + const auto eps_star_Pt = eps_star_atPts(pt,0); + const auto diff = eps_elm - eps_star_Pt; + const auto wPt = w(pt); + const auto dVPt = dV(pt); + sum += (diff * diff) * wPt * dVPt; + } + r_local += Kokkos::pow(Kokkos::sqrt(sum), ((2 * meshDim) / (2 * orderP + meshDim))); + errorNorm(elm) = Kokkos::pow(Kokkos::sqrt(sum), -(2 / (2 * orderP + meshDim))); + }, + r); // $\sum_{i=1}^n \|e_\epsilon\|^{\frac{2d}{2p+d}}$ + } + EstimationT& estimation; + OmegahMeshField& omf; + Kokkos::View errorNorm; // (||e_eps||_e)^(-2/(2p+d)) +}; + +//TODO move this into Estimation class +template +void computeSizeFactor(EstimationT& e, OmegahMeshField& omf, FieldElement& coordFe, ErrorT& errorIntegrator) { + SelfProduct sp(e,omf); + sp.process(coordFe); + const double epsStarNorm = Kokkos::sqrt(sp.r); + std::cout << "SelfProduct: " << epsStarNorm << "\n"; + const double a = e.tolerance * e.tolerance * // (n^hat)^2 + epsStarNorm * epsStarNorm; // ||e*||^2 + const double b = a / errorIntegrator.r; // term in parenthesis in section 4 of spr.tex + const double p = e.recovered_order; + e.size_factor = Kokkos::pow(b, 1.0 / (2.0 * p)); + std::cout << "size_factor: " << e.size_factor << "\n"; + +} + +/* computes h_e^new from section 4 of spr.tex */ +template +void getElementSizeField(EstimationT& e, ErrorT& errorIntegrator) { + Kokkos::View eSize("eSize", e.mesh.nelems()); + const auto errorNorm = errorIntegrator.errorNorm; + const auto size_factor = e.size_factor; + const auto currentElmSize = e.mesh.ask_sizes(); + e.mesh.template add_tag(e.mesh.dim(), "curElmSize", 1, currentElmSize); + Kokkos::parallel_for(e.mesh.nelems(), KOKKOS_LAMBDA(const int elm) { + const double h = currentElmSize[elm]; // h_e^current + eSize(elm) = h * errorNorm(elm) * size_factor; + }); + e.element_size = eSize; +} + +Kokkos::View +averageToVertex(Omega_h::Mesh& mesh, const Kokkos::View& elmSize) { + Omega_h::Write elmSize_oh(elmSize); + mesh.add_tag(mesh.dim(), "elmSizeField", 1, elmSize_oh); + Kokkos::View sizeField("sizeField", mesh.nverts()); + const auto v2e = mesh.ask_up(Omega_h::VERT, mesh.dim()); + const auto v2e_offsets = v2e.a2ab; + const auto v2e_values = v2e.ab2b; + Kokkos::parallel_for(mesh.nverts(), KOKKOS_LAMBDA(const int vtx) { + MeshField::Real s = 0; + for (auto idx = v2e_offsets[vtx]; idx < v2e_offsets[vtx + 1]; ++idx) { + const auto elm = v2e_values[idx]; + s += elmSize(elm); + } + const auto numUpElms = v2e_offsets[vtx+1] - v2e_offsets[vtx]; + sizeField(vtx) = s / numUpElms; + }); + return sizeField; +} + +template +Kokkos::View +getSprSizeField(EstimationT& e, OmegahMeshField& omf, FieldElement& coordFe) { + Error errorIntegrator(e,omf); + errorIntegrator.process(coordFe); + std::cout << "Error: " << errorIntegrator.r << "\n"; + computeSizeFactor(e, omf, coordFe, errorIntegrator); + getElementSizeField(e, errorIntegrator); + return averageToVertex(e.mesh, e.element_size); +} + +#endif