From 05cab3926da69962a76cd2fd981e86a5236cbfdf Mon Sep 17 00:00:00 2001 From: Cameron Smith Date: Sat, 15 Feb 2025 13:34:13 -0500 Subject: [PATCH 01/26] copy thwaites/gis driver from omegah source: git@github.com:SCOREC/omega_h cws/thwaitesAdaptDemo@aa45bffc --- test/CMakeLists.txt | 7 + test/thwaites_adapt.cpp | 282 ++++++++++++++++++++++++++++++++++++++++ 2 files changed, 289 insertions(+) create mode 100644 test/thwaites_adapt.cpp 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/thwaites_adapt.cpp b/test/thwaites_adapt.cpp new file mode 100644 index 00000000..487f2b69 --- /dev/null +++ b/test/thwaites_adapt.cpp @@ -0,0 +1,282 @@ +#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" //recover_hessians +#include //Omega_h::binary +#include //Omega_h::atomic_fetch_add +#include //ostringstream +#include //precision +#include +#include + +//detect floating point exceptions +#include + +using namespace Omega_h; + +const Real ICE_DENSITY = 910; //kg/m^3 + +template +Reals get_eigen_values_dim(Reals metrics) { + auto n = divide_no_remainder(metrics.size(), symm_ncomps(dim)); + auto isDegen = Write(n); + Write eigenVals(dim*n); + auto f = OMEGA_H_LAMBDA(LO i) { + auto m = get_symm(metrics, i); + auto m_dc = decompose_eigen(m); + Few m_ew_is_degen; + for (Int j = 0; j < dim; ++j) { + eigenVals[i*2+j] = m_dc.l[j]; + } + }; + parallel_for(n, f, "get_eigen_vals"); + return read(eigenVals); +} + +Reals get_eigen_values(Mesh* mesh, Reals metrics) { + if( mesh->dim() == 1) + return get_eigen_values_dim<1>(metrics); // is this supported? + else if( mesh->dim() == 2) + return get_eigen_values_dim<2>(metrics); + else if( mesh->dim() == 3) + return get_eigen_values_dim<3>(metrics); + else + return Reals(0); +} + +template +void printTagInfo(Omega_h::Mesh& mesh, std::ostringstream& oss, int dim, int tag, std::string type) { + auto tagbase = mesh.get_tag(dim, tag); + auto array = Omega_h::as(tagbase)->array(); + + Omega_h::Real min = get_min(array); + Omega_h::Real max = get_max(array); + + oss << std::setw(18) << std::left << tagbase->name().c_str() + << std::setw(5) << std::left << dim + << std::setw(7) << std::left << type + << std::setw(5) << std::left << tagbase->ncomps() + << std::setw(10) << std::left << min + << std::setw(10) << std::left << max + << "\n"; +} + +Reals getVolumeFromIceThickness(Omega_h::Mesh& mesh) { + auto coords = mesh.coords(); + auto triArea = mesh.ask_sizes(); + auto iceThickness = mesh.get_array(OMEGA_H_VERT, "ice_thickness"); + auto faces2verts = mesh.ask_elem_verts(); + Write vol(mesh.nfaces()); + auto getVol = OMEGA_H_LAMBDA(LO faceIdx) { + auto triVerts = Omega_h::gather_verts<3>(faces2verts, faceIdx); + auto triThickness = Omega_h::gather_scalars<3>(iceThickness, triVerts); + Real avgThickness = 0; + for(LO i=0; i<3; i++) + avgThickness += triThickness[i]; + avgThickness /= 3; + vol[faceIdx] = triArea[faceIdx]*avgThickness; + }; + parallel_for(mesh.nfaces(), getVol); + return vol; +} + +Reals getIceMass(Mesh& mesh, Reals vol) { + return multiply_each_by(vol, ICE_DENSITY); +} + +void printTags(Mesh& mesh) { + std::ostringstream oss; + // always print two places to the right of the decimal + // for floating point types (i.e., imbalance) + oss.precision(2); + oss << std::fixed; + + if (!mesh.comm()->rank()) { + oss << "\nTag Properties by Dimension: (Name, Dim, Type, Number of Components, Min. Value, Max. Value)\n"; + for (int dim=0; dim <= mesh.dim(); dim++) + for (int tag=0; tag < mesh.ntags(dim); tag++) { + auto tagbase = mesh.get_tag(dim, tag); + if (tagbase->type() == OMEGA_H_I8) + printTagInfo(mesh, oss, dim, tag, "I8"); + if (tagbase->type() == OMEGA_H_I32) + printTagInfo(mesh, oss, dim, tag, "I32"); + if (tagbase->type() == OMEGA_H_I64) + printTagInfo(mesh, oss, dim, tag, "I64"); + if (tagbase->type() == OMEGA_H_F64) + printTagInfo(mesh, oss, dim, tag, "F64"); + } + + std::cout << oss.str(); + } + +} + +static void check_total_mass(Mesh& mesh) { + auto vol = getVolumeFromIceThickness(mesh); + auto masses = getIceMass(mesh, vol); + auto owned_masses = mesh.owned_array(mesh.dim(), masses, 1); + auto mass = get_sum(mesh.comm(), owned_masses); + if (!mesh.comm()->rank()) { + std::cout << "mass " << mass << '\n'; + } +} + +void printTriCount(Mesh* mesh) { + const auto nTri = mesh->nglobal_ents(2); + if (!mesh->comm()->rank()) + std::cout << "nTri: " << nTri << "\n"; +} + +void setupFieldTransfer(AdaptOpts& opts) { + opts.xfer_opts.type_map["velocity"] = OMEGA_H_LINEAR_INTERP; + opts.xfer_opts.type_map["ice_thickness"] = 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; + } +} + +Reals exponentialVelocity(Mesh& mesh) { + auto coords = mesh.coords(); + auto coords_h = HostRead(coords); + auto vel_h = HostWrite(mesh.nverts()); + for(int i=0; i(coords_h, i); + if(std::fabs(x[0]) < 100) { + vel_h[i] = 20 * std::exp(-0.02 * std::fabs(x[0])); + } else { + vel_h[i] = 0; + } + } + return Reals(vel_h); +} + +Reals squareVelocity(Mesh& mesh) { + auto coords = mesh.coords(); + auto coords_h = HostRead(coords); + auto vel_h = HostWrite(mesh.nverts()); + for(int i=0; i(coords_h, i); + vel_h[i] = 20 * (-1.0e-6 * std::pow(x[0],2)); + } + return Reals(vel_h); +} + +Reals normSquaredVelocity(Mesh& mesh) { + auto coords = mesh.coords(); + auto coords_h = HostRead(coords); + auto vel_h = HostWrite(mesh.nverts()); + for(int i=0; i(coords_h, i); + vel_h[i] = norm_squared(x); + } + return Reals(vel_h); +} + +Reals magnitudeOfVelocity(Mesh& mesh) { + auto vel_x = mesh.get_array(VERT, "solution_1"); + auto vel_y = mesh.get_array(VERT, "solution_2"); + auto vel_mag = Write(mesh.nverts()); + auto f = OMEGA_H_LAMBDA(LO i) { + Vector<2> v{vel_x[i], vel_y[i]}; + vel_mag[i] = norm(v); + }; + parallel_for(mesh.nverts(), f, "magnitude_of_velocity"); + return Read(vel_mag); +} + +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 != 3 ) { + fprintf(stderr, "Usage: %s inputMesh.osh outputMeshPrefix\n", argv[0]); + exit(EXIT_FAILURE); + } + auto world = lib.world(); + Omega_h::Mesh mesh(&lib); + Omega_h::binary::read(argv[1], world, &mesh); + + Omega_h::vtk::write_parallel("beforeClassFix_edges.vtk", &mesh, 1); + classify_by_angles(&mesh,60*(Omega_h::PI/180)); + + //create analytic velocity field + auto velocity = magnitudeOfVelocity(mesh); + mesh.add_tag(VERT, "velocity", 1, Reals(velocity)); + const auto minSizeThreshold = 0.5; + auto max_velocity = get_max(velocity)*minSizeThreshold; + auto norm_velocity = divide_each_by(velocity, max_velocity); + mesh.add_tag(VERT, "norm_velocity", 1, Reals(norm_velocity)); + const auto minScaleFactor = 0.25; + const auto maxScaleFactor = 4.0; + std::cout << "minScaleFacor " << minScaleFactor << " " << + "maxScaleFacor " << maxScaleFactor << "\n"; + auto scale = Write(mesh.nverts()); + auto set_metric_scaling = OMEGA_H_LAMBDA(LO i) { + const auto a = 3.85939018034835; + const auto b = 20.0; + const auto d = 0.25; + const auto x = norm_velocity[i]; + scale[i] = a * std::log(b*x+1) + d; + if( x > minSizeThreshold ) scale[i] = 12; + }; + parallel_for(mesh.nverts(), set_metric_scaling, "set_metric_scaling"); + auto scale_r = Read(scale); + mesh.add_tag(VERT, "scale", 1, scale_r); + + Omega_h::add_implied_metric_tag(&mesh); + auto metric = mesh.get_array(VERT, "metric"); + auto tgt_metric = multiply_each(metric, scale_r); + mesh.add_tag(VERT, "tgt_metric", 3, Reals(tgt_metric)); + + auto genopts = Omega_h::MetricInput(); + auto metric_scaling = 1.0; //1.0: no scaling + genopts.sources.push_back( + Omega_h::MetricSource{OMEGA_H_GIVEN, metric_scaling, "tgt_metric"}); + genopts.nsmoothing_steps = 2; + Omega_h::generate_target_metric_tag(&mesh, genopts); + + mesh.set_parting(OMEGA_H_GHOSTED); + + //DEBUG { + auto edge_lengths_source = measure_edges_metric(&mesh, mesh.get_array(VERT, "metric")); + mesh.add_tag(EDGE, "lengths_source", 1, edge_lengths_source); + auto edge_lengths_target = measure_edges_metric(&mesh, mesh.get_array(VERT, "target_metric")); + mesh.add_tag(EDGE, "lengths_target", 1, edge_lengths_target); + auto elen = measure_edges_real(&mesh); + mesh.add_tag(EDGE, "lengths_real", 1, elen); + Omega_h::vtk::write_parallel("beforeAdapt_edges.vtk", &mesh, 1); + // } + + Omega_h::AdaptOpts opts(&mesh); + setupFieldTransfer(opts); + + printTriCount(&mesh); + printTags(mesh); + Omega_h::vtk::write_parallel("beforeAdapt.vtk", &mesh, 2); + check_total_mass(mesh); + + while (Omega_h::approach_metric(&mesh, opts)) { + adapt(&mesh, opts); + } + + auto post_elen = measure_edges_real(&mesh); + mesh.add_tag(EDGE, "post_lengths_real", 1, post_elen); + + printTriCount(&mesh); + printTags(mesh); + check_total_mass(mesh); + const std::string vtkFileName = std::string(argv[2]) + ".vtk"; + Omega_h::vtk::write_parallel(vtkFileName, &mesh, 2); + const std::string vtkFileName_edges = std::string(argv[2]) + "_edges.vtk"; + Omega_h::vtk::write_parallel(vtkFileName_edges, &mesh, 1); + return 0; +} From 9bbfbcd7f83d55730a37e2040dff811dd8a90976 Mon Sep 17 00:00:00 2001 From: Cameron Smith Date: Sun, 23 Feb 2025 13:59:20 -0500 Subject: [PATCH 02/26] mls_interpolation on thwaites mesh interpolation degree 2 can't recover a degree 2 function with a minPatchSize of 8 or 16, something is wrong --- test/test_spr_meshfields.cpp | 11 +++++++---- 1 file changed, 7 insertions(+), 4 deletions(-) diff --git a/test/test_spr_meshfields.cpp b/test/test_spr_meshfields.cpp index 4d61b3c0..ce071d57 100644 --- a/test/test_spr_meshfields.cpp +++ b/test/test_spr_meshfields.cpp @@ -88,7 +88,7 @@ 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"; +// std::cout << "vtx " << i << "\n"; CHECK_THAT( host_exact_target_values[i], Catch::Matchers::WithinAbs(host_approx_target_values[i], tolerance)); @@ -106,7 +106,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(); @@ -162,11 +165,11 @@ TEST_CASE("meshfields_spr_test") std::cerr << "start " << interp_degree << ", " << func_degree << " \n"; int minPatchSize; if(interp_degree == 1) minPatchSize = 3; - if(interp_degree == 2) minPatchSize = 8; //why so large? + if(interp_degree == 2) minPatchSize = 16; //why so large? 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"); From bdd35fbabbfccc21b755ba561c0e283f311593e3 Mon Sep 17 00:00:00 2001 From: Cameron Smith Date: Sun, 23 Feb 2025 20:24:44 -0500 Subject: [PATCH 03/26] write bad target point and its patch --- test/test_spr_meshfields.cpp | 31 +++++++++++++++++++++++++++++-- 1 file changed, 29 insertions(+), 2 deletions(-) diff --git a/test/test_spr_meshfields.cpp b/test/test_spr_meshfields.cpp index ce071d57..0ef8f39d 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,26 @@ 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, 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 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, i); + } CHECK_THAT( host_exact_target_values[i], Catch::Matchers::WithinAbs(host_approx_target_values[i], tolerance)); @@ -165,7 +192,7 @@ TEST_CASE("meshfields_spr_test") std::cerr << "start " << interp_degree << ", " << func_degree << " \n"; int minPatchSize; if(interp_degree == 1) minPatchSize = 3; - if(interp_degree == 2) minPatchSize = 16; //why so large? + if(interp_degree == 2) minPatchSize = 8; //why so large? if(interp_degree == 3) minPatchSize = 10; //why so large? std::cerr << "minPatchSize " << minPatchSize << "\n"; auto patches = mesh.get_vtx_patches(minPatchSize); From e5dbd5531261b6a624fa5c198fb8da71b7dabab9 Mon Sep 17 00:00:00 2001 From: Cameron Smith Date: Sun, 23 Feb 2025 21:23:44 -0500 Subject: [PATCH 04/26] more debug --- test/test_spr_meshfields.cpp | 11 ++++++++--- 1 file changed, 8 insertions(+), 3 deletions(-) diff --git a/test/test_spr_meshfields.cpp b/test/test_spr_meshfields.cpp index 0ef8f39d..6fed65bc 100644 --- a/test/test_spr_meshfields.cpp +++ b/test/test_spr_meshfields.cpp @@ -32,7 +32,8 @@ void print_patches(Omega_h::Graph& patches) { } } -void print_patch_vals(Omega_h::Graph& patches_d, Reals vtxCoords_d, Reals elmSrcVals_d, size_t vtx) { +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); @@ -43,9 +44,13 @@ void print_patch_vals(Omega_h::Graph& patches_d, Reals vtxCoords_d, Reals elmSrc for(int patch=0; patch \n"; for (auto valIdx = offsets[patch]; valIdx < offsets[patch + 1]; ++valIdx) { auto patchElm = values[valIdx]; - std::cout << "elmSrcVal[" << patchElm << "] " << elmSrcVals[patchElm] << "\n"; + std::cout << patchElm << " " + << elmCentroids[patchElm*meshDim] << " " + << elmCentroids[patchElm*meshDim+1] << " " + << elmSrcVals[patchElm] << "\n"; } std::cout << "\n"; } @@ -114,7 +119,7 @@ void test(Mesh& mesh, Omega_h::Graph& patches, int degree, 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, i); + print_patch_vals(patches, target_coordinates, source_values, source_coordinates, i); } CHECK_THAT( host_exact_target_values[i], From 52952ffb447ebca9a7ab626592b2eaeaf52e1ea2 Mon Sep 17 00:00:00 2001 From: Cameron Smith Date: Sun, 23 Feb 2025 22:38:47 -0500 Subject: [PATCH 05/26] adapt driver: remove adapt calls, WIP --- test/thwaites_adapt.cpp | 198 +++++----------------------------------- 1 file changed, 21 insertions(+), 177 deletions(-) diff --git a/test/thwaites_adapt.cpp b/test/thwaites_adapt.cpp index 487f2b69..a8828a91 100644 --- a/test/thwaites_adapt.cpp +++ b/test/thwaites_adapt.cpp @@ -16,41 +16,14 @@ #include #include +#include +#include + //detect floating point exceptions #include using namespace Omega_h; -const Real ICE_DENSITY = 910; //kg/m^3 - -template -Reals get_eigen_values_dim(Reals metrics) { - auto n = divide_no_remainder(metrics.size(), symm_ncomps(dim)); - auto isDegen = Write(n); - Write eigenVals(dim*n); - auto f = OMEGA_H_LAMBDA(LO i) { - auto m = get_symm(metrics, i); - auto m_dc = decompose_eigen(m); - Few m_ew_is_degen; - for (Int j = 0; j < dim; ++j) { - eigenVals[i*2+j] = m_dc.l[j]; - } - }; - parallel_for(n, f, "get_eigen_vals"); - return read(eigenVals); -} - -Reals get_eigen_values(Mesh* mesh, Reals metrics) { - if( mesh->dim() == 1) - return get_eigen_values_dim<1>(metrics); // is this supported? - else if( mesh->dim() == 2) - return get_eigen_values_dim<2>(metrics); - else if( mesh->dim() == 3) - return get_eigen_values_dim<3>(metrics); - else - return Reals(0); -} - template void printTagInfo(Omega_h::Mesh& mesh, std::ostringstream& oss, int dim, int tag, std::string type) { auto tagbase = mesh.get_tag(dim, tag); @@ -68,29 +41,6 @@ void printTagInfo(Omega_h::Mesh& mesh, std::ostringstream& oss, int dim, int tag << "\n"; } -Reals getVolumeFromIceThickness(Omega_h::Mesh& mesh) { - auto coords = mesh.coords(); - auto triArea = mesh.ask_sizes(); - auto iceThickness = mesh.get_array(OMEGA_H_VERT, "ice_thickness"); - auto faces2verts = mesh.ask_elem_verts(); - Write vol(mesh.nfaces()); - auto getVol = OMEGA_H_LAMBDA(LO faceIdx) { - auto triVerts = Omega_h::gather_verts<3>(faces2verts, faceIdx); - auto triThickness = Omega_h::gather_scalars<3>(iceThickness, triVerts); - Real avgThickness = 0; - for(LO i=0; i<3; i++) - avgThickness += triThickness[i]; - avgThickness /= 3; - vol[faceIdx] = triArea[faceIdx]*avgThickness; - }; - parallel_for(mesh.nfaces(), getVol); - return vol; -} - -Reals getIceMass(Mesh& mesh, Reals vol) { - return multiply_each_by(vol, ICE_DENSITY); -} - void printTags(Mesh& mesh) { std::ostringstream oss; // always print two places to the right of the decimal @@ -118,16 +68,6 @@ void printTags(Mesh& mesh) { } -static void check_total_mass(Mesh& mesh) { - auto vol = getVolumeFromIceThickness(mesh); - auto masses = getIceMass(mesh, vol); - auto owned_masses = mesh.owned_array(mesh.dim(), masses, 1); - auto mass = get_sum(mesh.comm(), owned_masses); - if (!mesh.comm()->rank()) { - std::cout << "mass " << mass << '\n'; - } -} - void printTriCount(Mesh* mesh) { const auto nTri = mesh->nglobal_ents(2); if (!mesh->comm()->rank()) @@ -145,53 +85,23 @@ void setupFieldTransfer(AdaptOpts& opts) { } } -Reals exponentialVelocity(Mesh& mesh) { - auto coords = mesh.coords(); - auto coords_h = HostRead(coords); - auto vel_h = HostWrite(mesh.nverts()); - for(int i=0; i(coords_h, i); - if(std::fabs(x[0]) < 100) { - vel_h[i] = 20 * std::exp(-0.02 * std::fabs(x[0])); - } else { - vel_h[i] = 0; - } - } - return Reals(vel_h); -} - -Reals squareVelocity(Mesh& mesh) { - auto coords = mesh.coords(); - auto coords_h = HostRead(coords); - auto vel_h = HostWrite(mesh.nverts()); - for(int i=0; i(coords_h, i); - vel_h[i] = 20 * (-1.0e-6 * std::pow(x[0],2)); - } - return Reals(vel_h); -} - -Reals normSquaredVelocity(Mesh& mesh) { - auto coords = mesh.coords(); - auto coords_h = HostRead(coords); - auto vel_h = HostWrite(mesh.nverts()); - for(int i=0; i(coords_h, i); - vel_h[i] = norm_squared(x); - } - return Reals(vel_h); -} - -Reals magnitudeOfVelocity(Mesh& mesh) { - auto vel_x = mesh.get_array(VERT, "solution_1"); - auto vel_y = mesh.get_array(VERT, "solution_2"); - auto vel_mag = Write(mesh.nverts()); - auto f = OMEGA_H_LAMBDA(LO i) { - Vector<2> v{vel_x[i], vel_y[i]}; - vel_mag[i] = norm(v); - }; - parallel_for(mesh.nverts(), f, "magnitude_of_velocity"); - return Read(vel_mag); +//Reals recoverField(Mesh& mesh) { +// const auto interpDegree = 2; +// auto patches = mesh.get_vtx_patches(minPatchSize); +// int dim = mesh.dim(); +// Omega_h::Write ignored(patches.ab2b.size(), 1); +// SupportResults support{patches.a2ab,patches.ab2b,ignored}; +// auto approx_target_values = +// pcms::mls_interpolation(source_values, source_coordinates, target_coordinates, +// support, dim, degree, support.radii2, pcms::RadialBasisFunction::NO_OP); +// return approx_target_values; +//} + +/** + * for now this will be a stub that returns an analytic + * linear field at elm centroids + */ +Reals getEffectiveStrainRate(Mesh& mesh) { } int main(int argc, char** argv) { @@ -206,74 +116,8 @@ int main(int argc, char** argv) { Omega_h::binary::read(argv[1], world, &mesh); Omega_h::vtk::write_parallel("beforeClassFix_edges.vtk", &mesh, 1); - classify_by_angles(&mesh,60*(Omega_h::PI/180)); - - //create analytic velocity field - auto velocity = magnitudeOfVelocity(mesh); - mesh.add_tag(VERT, "velocity", 1, Reals(velocity)); - const auto minSizeThreshold = 0.5; - auto max_velocity = get_max(velocity)*minSizeThreshold; - auto norm_velocity = divide_each_by(velocity, max_velocity); - mesh.add_tag(VERT, "norm_velocity", 1, Reals(norm_velocity)); - const auto minScaleFactor = 0.25; - const auto maxScaleFactor = 4.0; - std::cout << "minScaleFacor " << minScaleFactor << " " << - "maxScaleFacor " << maxScaleFactor << "\n"; - auto scale = Write(mesh.nverts()); - auto set_metric_scaling = OMEGA_H_LAMBDA(LO i) { - const auto a = 3.85939018034835; - const auto b = 20.0; - const auto d = 0.25; - const auto x = norm_velocity[i]; - scale[i] = a * std::log(b*x+1) + d; - if( x > minSizeThreshold ) scale[i] = 12; - }; - parallel_for(mesh.nverts(), set_metric_scaling, "set_metric_scaling"); - auto scale_r = Read(scale); - mesh.add_tag(VERT, "scale", 1, scale_r); - - Omega_h::add_implied_metric_tag(&mesh); - auto metric = mesh.get_array(VERT, "metric"); - auto tgt_metric = multiply_each(metric, scale_r); - mesh.add_tag(VERT, "tgt_metric", 3, Reals(tgt_metric)); - - auto genopts = Omega_h::MetricInput(); - auto metric_scaling = 1.0; //1.0: no scaling - genopts.sources.push_back( - Omega_h::MetricSource{OMEGA_H_GIVEN, metric_scaling, "tgt_metric"}); - genopts.nsmoothing_steps = 2; - Omega_h::generate_target_metric_tag(&mesh, genopts); - - mesh.set_parting(OMEGA_H_GHOSTED); - - //DEBUG { - auto edge_lengths_source = measure_edges_metric(&mesh, mesh.get_array(VERT, "metric")); - mesh.add_tag(EDGE, "lengths_source", 1, edge_lengths_source); - auto edge_lengths_target = measure_edges_metric(&mesh, mesh.get_array(VERT, "target_metric")); - mesh.add_tag(EDGE, "lengths_target", 1, edge_lengths_target); - auto elen = measure_edges_real(&mesh); - mesh.add_tag(EDGE, "lengths_real", 1, elen); - Omega_h::vtk::write_parallel("beforeAdapt_edges.vtk", &mesh, 1); - // } - - Omega_h::AdaptOpts opts(&mesh); - setupFieldTransfer(opts); - - printTriCount(&mesh); - printTags(mesh); - Omega_h::vtk::write_parallel("beforeAdapt.vtk", &mesh, 2); - check_total_mass(mesh); - - while (Omega_h::approach_metric(&mesh, opts)) { - adapt(&mesh, opts); - } - - auto post_elen = measure_edges_real(&mesh); - mesh.add_tag(EDGE, "post_lengths_real", 1, post_elen); + classify_by_angles(&mesh,60*(Omega_h::PI/180)); //NEEDED? - printTriCount(&mesh); - printTags(mesh); - check_total_mass(mesh); const std::string vtkFileName = std::string(argv[2]) + ".vtk"; Omega_h::vtk::write_parallel(vtkFileName, &mesh, 2); const std::string vtkFileName_edges = std::string(argv[2]) + "_edges.vtk"; From 8c779b9534404cad64ac9e0154b5e61e8beda56a Mon Sep 17 00:00:00 2001 From: Cameron Smith Date: Mon, 24 Feb 2025 16:35:33 -0500 Subject: [PATCH 06/26] get effective strain from mesh classify by angles was throwing floating point exceptions... the classification exists from the side sets on the exodus source mesh seems fine --- test/thwaites_adapt.cpp | 10 +++++++--- 1 file changed, 7 insertions(+), 3 deletions(-) diff --git a/test/thwaites_adapt.cpp b/test/thwaites_adapt.cpp index a8828a91..e6c24891 100644 --- a/test/thwaites_adapt.cpp +++ b/test/thwaites_adapt.cpp @@ -98,10 +98,13 @@ void setupFieldTransfer(AdaptOpts& opts) { //} /** - * for now this will be a stub that returns an analytic - * linear field at elm centroids + * 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"); } int main(int argc, char** argv) { @@ -116,7 +119,8 @@ int main(int argc, char** argv) { Omega_h::binary::read(argv[1], world, &mesh); Omega_h::vtk::write_parallel("beforeClassFix_edges.vtk", &mesh, 1); - classify_by_angles(&mesh,60*(Omega_h::PI/180)); //NEEDED? + + auto effectiveStrain = getEffectiveStrainRate(mesh); const std::string vtkFileName = std::string(argv[2]) + ".vtk"; Omega_h::vtk::write_parallel(vtkFileName, &mesh, 2); From 9dc43f0dfdbbfc5566c68482651897be3ca251a0 Mon Sep 17 00:00:00 2001 From: Cameron Smith Date: Mon, 24 Feb 2025 21:40:41 -0500 Subject: [PATCH 07/26] use omega_h project_by_fit should recover a linear field from the element effective strain --- test/thwaites_adapt.cpp | 20 +++++++------------- 1 file changed, 7 insertions(+), 13 deletions(-) diff --git a/test/thwaites_adapt.cpp b/test/thwaites_adapt.cpp index e6c24891..fcfcbeef 100644 --- a/test/thwaites_adapt.cpp +++ b/test/thwaites_adapt.cpp @@ -8,7 +8,7 @@ #include "Omega_h_for.hpp" #include "Omega_h_shape.hpp" #include "Omega_h_timer.hpp" -#include "Omega_h_recover.hpp" //recover_hessians +#include "Omega_h_recover.hpp" //project_by_fit #include //Omega_h::binary #include //Omega_h::atomic_fetch_add #include //ostringstream @@ -85,18 +85,6 @@ void setupFieldTransfer(AdaptOpts& opts) { } } -//Reals recoverField(Mesh& mesh) { -// const auto interpDegree = 2; -// auto patches = mesh.get_vtx_patches(minPatchSize); -// int dim = mesh.dim(); -// Omega_h::Write ignored(patches.ab2b.size(), 1); -// SupportResults support{patches.a2ab,patches.ab2b,ignored}; -// auto approx_target_values = -// pcms::mls_interpolation(source_values, source_coordinates, target_coordinates, -// support, dim, degree, support.radii2, pcms::RadialBasisFunction::NO_OP); -// return approx_target_values; -//} - /** * retrieve the effective strain rate from the mesh * @@ -107,6 +95,10 @@ Reals getEffectiveStrainRate(Mesh& mesh) { return mesh.get_array(2, "effective_stress"); } +Reals recoverLinearStrain(Mesh& mesh, Reals effectiveStrain) { + return project_by_fit(&mesh, effectiveStrain); +} + 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); @@ -121,6 +113,8 @@ int main(int argc, char** argv) { Omega_h::vtk::write_parallel("beforeClassFix_edges.vtk", &mesh, 1); auto effectiveStrain = getEffectiveStrainRate(mesh); + auto recoveredStrain = recoverLinearStrain(mesh,effectiveStrain); + mesh.add_tag(VERT, "recoveredStrain", 1, recoveredStrain); const std::string vtkFileName = std::string(argv[2]) + ".vtk"; Omega_h::vtk::write_parallel(vtkFileName, &mesh, 2); From ca33b7af61b5f27503f32ad6911303fae1d08913 Mon Sep 17 00:00:00 2001 From: Cameron Smith Date: Mon, 24 Feb 2025 22:43:18 -0500 Subject: [PATCH 08/26] WIP: start porting from sprEstimateError, does not compile --- test/thwaites_adapt.cpp | 101 ++++++++++++++++++++++++++++++++++++++-- 1 file changed, 98 insertions(+), 3 deletions(-) diff --git a/test/thwaites_adapt.cpp b/test/thwaites_adapt.cpp index fcfcbeef..4f12d47e 100644 --- a/test/thwaites_adapt.cpp +++ b/test/thwaites_adapt.cpp @@ -16,14 +16,105 @@ #include #include -#include -#include - //detect floating point exceptions #include using namespace Omega_h; +/* 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; +}; + +struct Estimation { + 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 */ + Omega_h::Reals 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 */ + Omega_h::Reals 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; +}; + +/* useful for initializing values to quickly + detect "uninitialized" value bugs */ +static double getNaN() +{ + return std::numeric_limits::quiet_NaN(); +} + +static void setupEstimation(Mesh& mesh, Estimation* e, const Reals eps, Real tolerance) +{ + /* 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 */ + e->integration_order = apf::getShape(eps)->getOrder(); //FIXME + e->mesh = mesh; + /* so far recovery order is directly tied to the + mesh's coordinate field order, coordinate this + with field recovery code */ + e->recovered_order = e->mesh->getShape()->getOrder(); //FIXME + e->eps = eps; + e->tolerance = tolerance; + e->size_factor = getNaN(); +} + +/* computes $\|f\|^2$ */ +class SelfProduct : public SInt +{ + public: + SelfProduct(Estimation* e): + SInt(e->integration_order), estimation(e), element(0) + { + v.setSize(apf::countComponents(e->eps_star)); //FIXME Vector2? + } + void atPoint(apf::Vector3 const& p, double w, double dV) //FIXME bulk operation + { + apf::getComponents(element, p, &v[0]); //FIXME get components of eps_star + r += (v * v) * w * dV; //FIXME vector*vector needed + } + private: + Estimation* estimation; + apf::Element* element; + apf::DynamicVector v; +}; + + + + + template void printTagInfo(Omega_h::Mesh& mesh, std::ostringstream& oss, int dim, int tag, std::string type) { auto tagbase = mesh.get_tag(dim, tag); @@ -99,6 +190,10 @@ Reals recoverLinearStrain(Mesh& mesh, Reals effectiveStrain) { return project_by_fit(&mesh, effectiveStrain); } +Reals computeError(Mesh& mesh, Reals effectiveStrain, Reals recoveredStrain) { + return error; +} + 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); From eda1cefbb30953761224e7db9bd7e8c513d719ba Mon Sep 17 00:00:00 2001 From: Cameron Smith Date: Tue, 25 Feb 2025 15:44:44 -0500 Subject: [PATCH 09/26] compiles, runtime error on parametric coords --- test/thwaites_adapt.cpp | 181 ++++++------------------------- test/thwaites_errorEstimator.hpp | 129 ++++++++++++++++++++++ 2 files changed, 165 insertions(+), 145 deletions(-) create mode 100644 test/thwaites_errorEstimator.hpp diff --git a/test/thwaites_adapt.cpp b/test/thwaites_adapt.cpp index 4f12d47e..161899b3 100644 --- a/test/thwaites_adapt.cpp +++ b/test/thwaites_adapt.cpp @@ -16,154 +16,15 @@ #include #include +#include "thwaites_errorEstimator.hpp" + //detect floating point exceptions #include -using namespace Omega_h; - -/* 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; -}; - -struct Estimation { - 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 */ - Omega_h::Reals 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 */ - Omega_h::Reals 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; -}; - -/* useful for initializing values to quickly - detect "uninitialized" value bugs */ -static double getNaN() -{ - return std::numeric_limits::quiet_NaN(); -} - -static void setupEstimation(Mesh& mesh, Estimation* e, const Reals eps, Real tolerance) -{ - /* 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 */ - e->integration_order = apf::getShape(eps)->getOrder(); //FIXME - e->mesh = mesh; - /* so far recovery order is directly tied to the - mesh's coordinate field order, coordinate this - with field recovery code */ - e->recovered_order = e->mesh->getShape()->getOrder(); //FIXME - e->eps = eps; - e->tolerance = tolerance; - e->size_factor = getNaN(); -} - -/* computes $\|f\|^2$ */ -class SelfProduct : public SInt -{ - public: - SelfProduct(Estimation* e): - SInt(e->integration_order), estimation(e), element(0) - { - v.setSize(apf::countComponents(e->eps_star)); //FIXME Vector2? - } - void atPoint(apf::Vector3 const& p, double w, double dV) //FIXME bulk operation - { - apf::getComponents(element, p, &v[0]); //FIXME get components of eps_star - r += (v * v) * w * dV; //FIXME vector*vector needed - } - private: - Estimation* estimation; - apf::Element* element; - apf::DynamicVector v; -}; - - - - - -template -void printTagInfo(Omega_h::Mesh& mesh, std::ostringstream& oss, int dim, int tag, std::string type) { - auto tagbase = mesh.get_tag(dim, tag); - auto array = Omega_h::as(tagbase)->array(); - - Omega_h::Real min = get_min(array); - Omega_h::Real max = get_max(array); - - oss << std::setw(18) << std::left << tagbase->name().c_str() - << std::setw(5) << std::left << dim - << std::setw(7) << std::left << type - << std::setw(5) << std::left << tagbase->ncomps() - << std::setw(10) << std::left << min - << std::setw(10) << std::left << max - << "\n"; -} - -void printTags(Mesh& mesh) { - std::ostringstream oss; - // always print two places to the right of the decimal - // for floating point types (i.e., imbalance) - oss.precision(2); - oss << std::fixed; - - if (!mesh.comm()->rank()) { - oss << "\nTag Properties by Dimension: (Name, Dim, Type, Number of Components, Min. Value, Max. Value)\n"; - for (int dim=0; dim <= mesh.dim(); dim++) - for (int tag=0; tag < mesh.ntags(dim); tag++) { - auto tagbase = mesh.get_tag(dim, tag); - if (tagbase->type() == OMEGA_H_I8) - printTagInfo(mesh, oss, dim, tag, "I8"); - if (tagbase->type() == OMEGA_H_I32) - printTagInfo(mesh, oss, dim, tag, "I32"); - if (tagbase->type() == OMEGA_H_I64) - printTagInfo(mesh, oss, dim, tag, "I64"); - if (tagbase->type() == OMEGA_H_F64) - printTagInfo(mesh, oss, dim, tag, "F64"); - } - - std::cout << oss.str(); - } - -} +using ExecutionSpace = Kokkos::DefaultExecutionSpace; +using MemorySpace = Kokkos::DefaultExecutionSpace::memory_space; -void printTriCount(Mesh* mesh) { - const auto nTri = mesh->nglobal_ents(2); - if (!mesh->comm()->rank()) - std::cout << "nTri: " << nTri << "\n"; -} +using namespace Omega_h; void setupFieldTransfer(AdaptOpts& opts) { opts.xfer_opts.type_map["velocity"] = OMEGA_H_LINEAR_INTERP; @@ -191,7 +52,17 @@ Reals recoverLinearStrain(Mesh& mesh, Reals effectiveStrain) { } Reals computeError(Mesh& mesh, Reals effectiveStrain, Reals recoveredStrain) { - return error; + 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"); } int main(int argc, char** argv) { @@ -210,6 +81,26 @@ int main(int argc, char** argv) { 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); + + const auto tolerance = 1e-6; //FIXME - set to 'adaptRatio' in apf::spr + auto estimation = Estimation(mesh, effectiveStrain, + recoveredStrainField, tolerance); + + SelfProduct sp(estimation,omf); + sp.process(coordFe); const std::string vtkFileName = std::string(argv[2]) + ".vtk"; Omega_h::vtk::write_parallel(vtkFileName, &mesh, 2); diff --git a/test/thwaites_errorEstimator.hpp b/test/thwaites_errorEstimator.hpp new file mode 100644 index 00000000..07afda56 --- /dev/null +++ b/test/thwaites_errorEstimator.hpp @@ -0,0 +1,129 @@ +#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 */ + Omega_h::Reals 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.size()/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; +}; + + + +#endif From c408ae54fb43d323067a60ed668eb6e8bb98b165 Mon Sep 17 00:00:00 2001 From: Cameron Smith Date: Tue, 25 Feb 2025 22:08:57 -0500 Subject: [PATCH 10/26] executes --- test/thwaites_errorEstimator.hpp | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/test/thwaites_errorEstimator.hpp b/test/thwaites_errorEstimator.hpp index 07afda56..67ad2e7e 100644 --- a/test/thwaites_errorEstimator.hpp +++ b/test/thwaites_errorEstimator.hpp @@ -40,11 +40,11 @@ class Estimation { /* 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 + 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 + recovered_order = 1; //FIXME } Omega_h::Mesh& mesh; /* the maximum polynomial order that can be @@ -99,7 +99,7 @@ class SelfProduct : public SInt { std::cerr << "SelfProduct::atPoints(...)\n"; - const size_t numPtsPerElem = p.size()/estimation.mesh.nelems(); + const size_t numPtsPerElem = p.extent(0)/estimation.mesh.nelems(); auto eps_star_atPts = omf.triangleLocalPointEval(p,numPtsPerElem, estimation.eps_star); From c9702301f7589da629afd5651fc37ff7bc05f947 Mon Sep 17 00:00:00 2001 From: Cameron Smith Date: Wed, 26 Feb 2025 14:53:50 -0500 Subject: [PATCH 11/26] matches pumi result pumi: 7.70212*10^6 meshfields: 7.70811*10^6 difference: 5.990*10^3 percent_difference vs pumi: 7.78*10-4 --- test/thwaites_adapt.cpp | 1 + 1 file changed, 1 insertion(+) diff --git a/test/thwaites_adapt.cpp b/test/thwaites_adapt.cpp index 161899b3..2b21522d 100644 --- a/test/thwaites_adapt.cpp +++ b/test/thwaites_adapt.cpp @@ -101,6 +101,7 @@ int main(int argc, char** argv) { SelfProduct sp(estimation,omf); sp.process(coordFe); + std::cout << "SelfProduct: " << Kokkos::sqrt(sp.r) << "\n"; const std::string vtkFileName = std::string(argv[2]) + ".vtk"; Omega_h::vtk::write_parallel(vtkFileName, &mesh, 2); From 7eb02d79b4cd4392ee935d648a1cbe4be44a2870 Mon Sep 17 00:00:00 2001 From: Cameron Smith Date: Wed, 26 Feb 2025 15:01:47 -0500 Subject: [PATCH 12/26] WIP computeSizeFactor --- test/thwaites_adapt.cpp | 4 +--- test/thwaites_errorEstimator.hpp | 7 ++++++- 2 files changed, 7 insertions(+), 4 deletions(-) diff --git a/test/thwaites_adapt.cpp b/test/thwaites_adapt.cpp index 2b21522d..cd02ee7d 100644 --- a/test/thwaites_adapt.cpp +++ b/test/thwaites_adapt.cpp @@ -99,9 +99,7 @@ int main(int argc, char** argv) { auto estimation = Estimation(mesh, effectiveStrain, recoveredStrainField, tolerance); - SelfProduct sp(estimation,omf); - sp.process(coordFe); - std::cout << "SelfProduct: " << Kokkos::sqrt(sp.r) << "\n"; + computeSizeFactor(estimation, omf, coordFe); const std::string vtkFileName = std::string(argv[2]) + ".vtk"; Omega_h::vtk::write_parallel(vtkFileName, &mesh, 2); diff --git a/test/thwaites_errorEstimator.hpp b/test/thwaites_errorEstimator.hpp index 67ad2e7e..e57ff7c6 100644 --- a/test/thwaites_errorEstimator.hpp +++ b/test/thwaites_errorEstimator.hpp @@ -124,6 +124,11 @@ class SelfProduct : public SInt OmegahMeshField& omf; }; - +template +void computeSizeFactor(EstimationT& e, OmegahMeshField& omf, FieldElement& coordFe) { + SelfProduct sp(e,omf); + sp.process(coordFe); + std::cout << "SelfProduct: " << Kokkos::sqrt(sp.r) << "\n"; +} #endif From 7838f5bce4431d61bc5c4c16dc49f9365f9126e5 Mon Sep 17 00:00:00 2001 From: Cameron Smith Date: Wed, 26 Feb 2025 16:04:34 -0500 Subject: [PATCH 13/26] results are close to pumi pumi: scorec/core cws/sprThwaites @ e1f0c79d, from modified test/spr_test.cc epsStarNorm 7.70212e+06 Error: 5.66159e+07 size_factor: 102.363 meshfields: SelfProduct: 7.70811e+06 Error: 5.63518e+07 size_factor: 102.682 --- test/thwaites_adapt.cpp | 4 +- test/thwaites_errorEstimator.hpp | 69 +++++++++++++++++++++++++++++++- 2 files changed, 70 insertions(+), 3 deletions(-) diff --git a/test/thwaites_adapt.cpp b/test/thwaites_adapt.cpp index cd02ee7d..9e91e06c 100644 --- a/test/thwaites_adapt.cpp +++ b/test/thwaites_adapt.cpp @@ -95,9 +95,9 @@ int main(int argc, char** argv) { const auto [shp, map] = MeshField::Omegah::getTriangleElement(mesh); MeshField::FieldElement coordFe(mesh.nelems(), coordField, shp, map); - const auto tolerance = 1e-6; //FIXME - set to 'adaptRatio' in apf::spr + const auto adaptRatio = 0.1; //same value used in scorec/core:cws/sprThwaites test/spr_test.cc auto estimation = Estimation(mesh, effectiveStrain, - recoveredStrainField, tolerance); + recoveredStrainField, adaptRatio); computeSizeFactor(estimation, omf, coordFe); diff --git a/test/thwaites_errorEstimator.hpp b/test/thwaites_errorEstimator.hpp index e57ff7c6..9b11c032 100644 --- a/test/thwaites_errorEstimator.hpp +++ b/test/thwaites_errorEstimator.hpp @@ -124,11 +124,78 @@ class SelfProduct : public SInt 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) + { + } + 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 += pow(sqrt(sum), ((2 * meshDim) / (2 * orderP + meshDim))); + }, + r); + + } + EstimationT& estimation; + OmegahMeshField& omf; +}; + template void computeSizeFactor(EstimationT& e, OmegahMeshField& omf, FieldElement& coordFe) { SelfProduct sp(e,omf); sp.process(coordFe); - std::cout << "SelfProduct: " << Kokkos::sqrt(sp.r) << "\n"; + const double epsStarNorm = Kokkos::sqrt(sp.r); + std::cout << "SelfProduct: " << epsStarNorm << "\n"; + + Error errorIntegrator(e,omf); + errorIntegrator.process(coordFe); + std::cout << "Error: " << errorIntegrator.r << "\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 = pow(b, 1.0 / (2.0 * p)); + std::cout << "size_factor: " << e.size_factor << "\n"; + } #endif From 1ddf3956a849e56dd1ad7b30c07519283ce5237d Mon Sep 17 00:00:00 2001 From: Cameron Smith Date: Thu, 27 Feb 2025 12:34:25 -0500 Subject: [PATCH 14/26] compute other error term --- test/thwaites_errorEstimator.hpp | 18 +++++++++++++----- 1 file changed, 13 insertions(+), 5 deletions(-) diff --git a/test/thwaites_errorEstimator.hpp b/test/thwaites_errorEstimator.hpp index 9b11c032..c01a4dbb 100644 --- a/test/thwaites_errorEstimator.hpp +++ b/test/thwaites_errorEstimator.hpp @@ -134,7 +134,8 @@ class Error : public SInt Error(EstimationT& estimation_in, OmegahMeshField& omf_in): SInt(estimation_in.integration_order), estimation(estimation_in), - omf(omf_in) + omf(omf_in), + errorNorm("errorNorm", estimation_in.mesh.nelems()) { } void atPoints(Kokkos::View p, @@ -169,15 +170,17 @@ class Error : public SInt const auto dVPt = dV(pt); sum += (diff * diff) * wPt * dVPt; } - r_local += pow(sqrt(sum), ((2 * meshDim) / (2 * orderP + meshDim))); + r_local += Kokkos::pow(Kokkos::sqrt(sum), ((2 * meshDim) / (2 * orderP + meshDim))); + errorNorm(elm) = Kokkos::pow(Kokkos::sqrt(sum), -(2 / (2 * orderP + meshDim))); }, - r); - + r); // $\sum_{i=1}^n \|e_\epsilon\|^{\frac{2d}{2p+d}}$ } EstimationT& estimation; OmegahMeshField& omf; + Kokkos::View errorNorm; // $\|e_\epsilon\|^{-\frac{2}{2p+d}}_e$ }; +//TODO move this into Estimation class template void computeSizeFactor(EstimationT& e, OmegahMeshField& omf, FieldElement& coordFe) { SelfProduct sp(e,omf); @@ -193,9 +196,14 @@ void computeSizeFactor(EstimationT& e, OmegahMeshField& omf, FieldElement& coord 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 = pow(b, 1.0 / (2.0 * p)); + e.size_factor = Kokkos::pow(b, 1.0 / (2.0 * p)); std::cout << "size_factor: " << e.size_factor << "\n"; } +template +void estimateError(EstimationT& e, OmegahMeshField& omf, FieldElement& coordFe) { + computeSizeFactor(e, omf, coordFe); +} + #endif From 240fba77468769d38c72af075e14cb56c1373fda Mon Sep 17 00:00:00 2001 From: Cameron Smith Date: Thu, 27 Feb 2025 12:50:50 -0500 Subject: [PATCH 15/26] compute size field --- test/thwaites_adapt.cpp | 2 +- test/thwaites_errorEstimator.hpp | 34 ++++++++++++++++++++++---------- 2 files changed, 25 insertions(+), 11 deletions(-) diff --git a/test/thwaites_adapt.cpp b/test/thwaites_adapt.cpp index 9e91e06c..e5b13242 100644 --- a/test/thwaites_adapt.cpp +++ b/test/thwaites_adapt.cpp @@ -99,7 +99,7 @@ int main(int argc, char** argv) { auto estimation = Estimation(mesh, effectiveStrain, recoveredStrainField, adaptRatio); - computeSizeFactor(estimation, omf, coordFe); + estimateError(estimation, omf, coordFe); const std::string vtkFileName = std::string(argv[2]) + ".vtk"; Omega_h::vtk::write_parallel(vtkFileName, &mesh, 2); diff --git a/test/thwaites_errorEstimator.hpp b/test/thwaites_errorEstimator.hpp index c01a4dbb..b9dc9d04 100644 --- a/test/thwaites_errorEstimator.hpp +++ b/test/thwaites_errorEstimator.hpp @@ -76,7 +76,7 @@ class Estimation { desired element size = current size * current error * size_factor */ Omega_h::Real size_factor; /* a temporary field storing desired sizes at elements */ - Omega_h::Reals element_size; + 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; @@ -177,21 +177,16 @@ class Error : public SInt } EstimationT& estimation; OmegahMeshField& omf; - Kokkos::View errorNorm; // $\|e_\epsilon\|^{-\frac{2}{2p+d}}_e$ + Kokkos::View errorNorm; // (||e_eps||_e)^(-2/(2p+d)) }; //TODO move this into Estimation class -template -void computeSizeFactor(EstimationT& e, OmegahMeshField& omf, FieldElement& coordFe) { +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"; - - Error errorIntegrator(e,omf); - errorIntegrator.process(coordFe); - std::cout << "Error: " << errorIntegrator.r << "\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 @@ -201,9 +196,28 @@ void computeSizeFactor(EstimationT& e, OmegahMeshField& omf, FieldElement& coord } +/* 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(); + Kokkos::parallel_for(e.mesh.nelems(), KOKKOS_LAMBDA(const int elm) { + const double h = currentElmSize[elm]; // h_e^current //FIXME + eSize(elm) = h * errorNorm(elm) * size_factor; + }); + e.element_size = eSize; +} + template void estimateError(EstimationT& e, OmegahMeshField& omf, FieldElement& coordFe) { - computeSizeFactor(e, omf, coordFe); + Error errorIntegrator(e,omf); + errorIntegrator.process(coordFe); + std::cout << "Error: " << errorIntegrator.r << "\n"; + computeSizeFactor(e, omf, coordFe, errorIntegrator); + getElementSizeField(e, errorIntegrator); +// averageSizeField(e); } #endif From e43d2e133ad097fa83044cc49fe90668d2f683de Mon Sep 17 00:00:00 2001 From: Cameron Smith Date: Thu, 27 Feb 2025 13:32:02 -0500 Subject: [PATCH 16/26] average elm size field to vertices --- test/thwaites_adapt.cpp | 2 +- test/thwaites_errorEstimator.hpp | 22 ++++++++++++++++++++-- 2 files changed, 21 insertions(+), 3 deletions(-) diff --git a/test/thwaites_adapt.cpp b/test/thwaites_adapt.cpp index e5b13242..e03c3af6 100644 --- a/test/thwaites_adapt.cpp +++ b/test/thwaites_adapt.cpp @@ -99,7 +99,7 @@ int main(int argc, char** argv) { auto estimation = Estimation(mesh, effectiveStrain, recoveredStrainField, adaptRatio); - estimateError(estimation, omf, coordFe); + const auto sizeField = getSprSizeField(estimation, omf, coordFe); const std::string vtkFileName = std::string(argv[2]) + ".vtk"; Omega_h::vtk::write_parallel(vtkFileName, &mesh, 2); diff --git a/test/thwaites_errorEstimator.hpp b/test/thwaites_errorEstimator.hpp index b9dc9d04..5bcc4299 100644 --- a/test/thwaites_errorEstimator.hpp +++ b/test/thwaites_errorEstimator.hpp @@ -210,14 +210,32 @@ void getElementSizeField(EstimationT& e, ErrorT& errorIntegrator) { e.element_size = eSize; } +Kokkos::View +averageToVertex(Omega_h::Mesh& mesh, const Kokkos::View& elmSize) { + 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 -void estimateError(EstimationT& e, OmegahMeshField& omf, FieldElement& coordFe) { +void 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); -// averageSizeField(e); + return averageToVertex(e.mesh, e.element_size); } #endif From 7ec2a4c10286f766c08fc60e730c56d032f82c05 Mon Sep 17 00:00:00 2001 From: Cameron Smith Date: Thu, 27 Feb 2025 13:43:18 -0500 Subject: [PATCH 17/26] write size field to vtk --- test/thwaites_adapt.cpp | 2 ++ test/thwaites_errorEstimator.hpp | 3 ++- 2 files changed, 4 insertions(+), 1 deletion(-) diff --git a/test/thwaites_adapt.cpp b/test/thwaites_adapt.cpp index e03c3af6..87112f87 100644 --- a/test/thwaites_adapt.cpp +++ b/test/thwaites_adapt.cpp @@ -100,6 +100,8 @@ int main(int argc, char** argv) { recoveredStrainField, adaptRatio); const auto sizeField = getSprSizeField(estimation, omf, coordFe); + Omega_h::Write sizeField_oh(sizeField); + mesh.add_tag(VERT, "sizeField", 1, sizeField_oh); const std::string vtkFileName = std::string(argv[2]) + ".vtk"; Omega_h::vtk::write_parallel(vtkFileName, &mesh, 2); diff --git a/test/thwaites_errorEstimator.hpp b/test/thwaites_errorEstimator.hpp index 5bcc4299..2e182695 100644 --- a/test/thwaites_errorEstimator.hpp +++ b/test/thwaites_errorEstimator.hpp @@ -229,7 +229,8 @@ averageToVertex(Omega_h::Mesh& mesh, const Kokkos::View& elmSi } template -void getSprSizeField(EstimationT& e, OmegahMeshField& omf, FieldElement& coordFe) { +Kokkos::View +getSprSizeField(EstimationT& e, OmegahMeshField& omf, FieldElement& coordFe) { Error errorIntegrator(e,omf); errorIntegrator.process(coordFe); std::cout << "Error: " << errorIntegrator.r << "\n"; From 6857b6afd1299827b10112b5a3536ac43261788a Mon Sep 17 00:00:00 2001 From: Cameron Smith Date: Thu, 27 Feb 2025 13:51:02 -0500 Subject: [PATCH 18/26] write elm size field --- test/thwaites_errorEstimator.hpp | 2 ++ 1 file changed, 2 insertions(+) diff --git a/test/thwaites_errorEstimator.hpp b/test/thwaites_errorEstimator.hpp index 2e182695..ae5a9dfd 100644 --- a/test/thwaites_errorEstimator.hpp +++ b/test/thwaites_errorEstimator.hpp @@ -212,6 +212,8 @@ void getElementSizeField(EstimationT& e, ErrorT& errorIntegrator) { 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; From 7e908b4088e97a82e628cf7eace369d73b422389 Mon Sep 17 00:00:00 2001 From: Cameron Smith Date: Thu, 27 Feb 2025 16:14:55 -0500 Subject: [PATCH 19/26] remove FIXME --- test/thwaites_errorEstimator.hpp | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/test/thwaites_errorEstimator.hpp b/test/thwaites_errorEstimator.hpp index ae5a9dfd..dd8e090c 100644 --- a/test/thwaites_errorEstimator.hpp +++ b/test/thwaites_errorEstimator.hpp @@ -203,8 +203,9 @@ void getElementSizeField(EstimationT& e, ErrorT& errorIntegrator) { 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 //FIXME + const double h = currentElmSize[elm]; // h_e^current eSize(elm) = h * errorNorm(elm) * size_factor; }); e.element_size = eSize; From 4b4f4ac6c226bafaaa3622f1f6f9cdeb569e6f06 Mon Sep 17 00:00:00 2001 From: Cameron Smith Date: Thu, 27 Feb 2025 16:16:15 -0500 Subject: [PATCH 20/26] coarsens the 106k tri mesh down to 21k tri --- test/thwaites_adapt.cpp | 33 ++++++++++++++++++++++++++++----- 1 file changed, 28 insertions(+), 5 deletions(-) diff --git a/test/thwaites_adapt.cpp b/test/thwaites_adapt.cpp index 87112f87..9250a9d2 100644 --- a/test/thwaites_adapt.cpp +++ b/test/thwaites_adapt.cpp @@ -65,6 +65,12 @@ void setFieldAtVertices(Omega_h::Mesh &mesh, Reals recoveredStrain, ShapeField f 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); @@ -99,13 +105,30 @@ int main(int argc, char** argv) { auto estimation = Estimation(mesh, effectiveStrain, recoveredStrainField, adaptRatio); - const auto sizeField = getSprSizeField(estimation, omf, coordFe); - Omega_h::Write sizeField_oh(sizeField); - mesh.add_tag(VERT, "sizeField", 1, sizeField_oh); + const auto tgtLength = getSprSizeField(estimation, omf, coordFe); + Omega_h::Write tgtLength_oh(tgtLength); + mesh.add_tag(VERT, "tgtLength", 1, tgtLength_oh); - const std::string vtkFileName = std::string(argv[2]) + ".vtk"; + { //write vtk + const std::string vtkFileName = "beforeAdapt" + std::string(argv[2]) + ".vtk"; Omega_h::vtk::write_parallel(vtkFileName, &mesh, 2); - const std::string vtkFileName_edges = std::string(argv[2]) + "_edges.vtk"; + const std::string vtkFileName_edges = "beforeAdapt" + std::string(argv[2]) + "_edges.vtk"; Omega_h::vtk::write_parallel(vtkFileName_edges, &mesh, 1); + } + + //adapt + auto opts = Omega_h::AdaptOpts(&mesh); + //setupFieldTransfer(opts); //FIXME - add this back + printTriCount(&mesh); + + auto verbose = true; + const auto isos = Omega_h::isos_from_lengths(tgtLength_oh); + Omega_h::grade_fix_adapt(&mesh, opts, isos, verbose); + + { //write vtk + const std::string vtkFileName = "afterAdapt" + std::string(argv[2]) + ".vtk"; + Omega_h::vtk::write_parallel(vtkFileName, &mesh, 2); + } + return 0; } From 69bc43877fb387e29775297c04f72841600f2b02 Mon Sep 17 00:00:00 2001 From: Cameron Smith Date: Thu, 27 Feb 2025 16:36:01 -0500 Subject: [PATCH 21/26] adaptRatio arg and adaptRatio in output files --- test/thwaites_adapt.cpp | 18 ++++++++++++------ 1 file changed, 12 insertions(+), 6 deletions(-) diff --git a/test/thwaites_adapt.cpp b/test/thwaites_adapt.cpp index 9250a9d2..31820776 100644 --- a/test/thwaites_adapt.cpp +++ b/test/thwaites_adapt.cpp @@ -74,13 +74,19 @@ void printTriCount(Mesh* mesh) { 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 != 3 ) { - fprintf(stderr, "Usage: %s inputMesh.osh outputMeshPrefix\n", argv[0]); + if( argc != 4 ) { + fprintf(stderr, "Usage: %s inputMesh.osh outputMeshPrefix adaptRatio\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]); + std::cout << "input mesh: " << argv[1] << " outputMeshPrefix: " + << prefix << " adaptRatio: " << adaptRatio << "\n"; + const auto outname = prefix + "adaptRatio_" + std::string(argv[3]); Omega_h::vtk::write_parallel("beforeClassFix_edges.vtk", &mesh, 1); @@ -101,7 +107,6 @@ int main(int argc, char** argv) { const auto [shp, map] = MeshField::Omegah::getTriangleElement(mesh); MeshField::FieldElement coordFe(mesh.nelems(), coordField, shp, map); - const auto adaptRatio = 0.1; //same value used in scorec/core:cws/sprThwaites test/spr_test.cc auto estimation = Estimation(mesh, effectiveStrain, recoveredStrainField, adaptRatio); @@ -110,14 +115,15 @@ int main(int argc, char** argv) { mesh.add_tag(VERT, "tgtLength", 1, tgtLength_oh); { //write vtk - const std::string vtkFileName = "beforeAdapt" + std::string(argv[2]) + ".vtk"; + const std::string vtkFileName = "beforeAdapt" + outname + ".vtk"; Omega_h::vtk::write_parallel(vtkFileName, &mesh, 2); - const std::string vtkFileName_edges = "beforeAdapt" + std::string(argv[2]) + "_edges.vtk"; + 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); + //opts.max_length_allowed = 64; //made the mesh even coarser.... //setupFieldTransfer(opts); //FIXME - add this back printTriCount(&mesh); @@ -126,7 +132,7 @@ int main(int argc, char** argv) { Omega_h::grade_fix_adapt(&mesh, opts, isos, verbose); { //write vtk - const std::string vtkFileName = "afterAdapt" + std::string(argv[2]) + ".vtk"; + const std::string vtkFileName = "afterAdapt" + outname + ".vtk"; Omega_h::vtk::write_parallel(vtkFileName, &mesh, 2); } From 575dae4b2053f457c3f74d1edaf24cf272d2fd57 Mon Sep 17 00:00:00 2001 From: Cameron Smith Date: Thu, 27 Feb 2025 21:50:28 -0500 Subject: [PATCH 22/26] clamp_metric --- test/thwaites_adapt.cpp | 17 +++++++++++++---- 1 file changed, 13 insertions(+), 4 deletions(-) diff --git a/test/thwaites_adapt.cpp b/test/thwaites_adapt.cpp index 31820776..4f137393 100644 --- a/test/thwaites_adapt.cpp +++ b/test/thwaites_adapt.cpp @@ -84,9 +84,18 @@ int main(int argc, char** argv) { 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]); + auto max_size = Real(16); //roughly maintains the boundary shape + //in the downstream parts of the domain + //where there is significant coarsening + auto min_size = Real(1); //roughly the smallest edge length in the + //source mesh std::cout << "input mesh: " << argv[1] << " outputMeshPrefix: " - << prefix << " adaptRatio: " << adaptRatio << "\n"; - const auto outname = prefix + "adaptRatio_" + std::string(argv[3]); + << 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); Omega_h::vtk::write_parallel("beforeClassFix_edges.vtk", &mesh, 1); @@ -123,13 +132,13 @@ int main(int argc, char** argv) { //adapt auto opts = Omega_h::AdaptOpts(&mesh); - //opts.max_length_allowed = 64; //made the mesh even coarser.... //setupFieldTransfer(opts); //FIXME - add this back printTriCount(&mesh); auto verbose = true; const auto isos = Omega_h::isos_from_lengths(tgtLength_oh); - Omega_h::grade_fix_adapt(&mesh, opts, isos, verbose); + auto metric = clamp_metrics(mesh.nverts(), isos, min_size, max_size); + Omega_h::grade_fix_adapt(&mesh, opts, metric, verbose); { //write vtk const std::string vtkFileName = "afterAdapt" + outname + ".vtk"; From d800d43997062d91359ec3ba76b1c63d6f344852 Mon Sep 17 00:00:00 2001 From: Cameron Smith Date: Fri, 28 Feb 2025 14:29:31 -0500 Subject: [PATCH 23/26] transfer fields, write osh --- test/thwaites_adapt.cpp | 14 +++++++++----- 1 file changed, 9 insertions(+), 5 deletions(-) diff --git a/test/thwaites_adapt.cpp b/test/thwaites_adapt.cpp index 4f137393..0ee3814d 100644 --- a/test/thwaites_adapt.cpp +++ b/test/thwaites_adapt.cpp @@ -27,7 +27,8 @@ using MemorySpace = Kokkos::DefaultExecutionSpace::memory_space; using namespace Omega_h; void setupFieldTransfer(AdaptOpts& opts) { - opts.xfer_opts.type_map["velocity"] = OMEGA_H_LINEAR_INTERP; + 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; const int numLayers = 11; for(int i=1; i<=numLayers; i++) { @@ -132,7 +133,7 @@ int main(int argc, char** argv) { //adapt auto opts = Omega_h::AdaptOpts(&mesh); - //setupFieldTransfer(opts); //FIXME - add this back + setupFieldTransfer(opts); printTriCount(&mesh); auto verbose = true; @@ -140,10 +141,13 @@ int main(int argc, char** argv) { auto metric = clamp_metrics(mesh.nverts(), isos, min_size, max_size); Omega_h::grade_fix_adapt(&mesh, opts, metric, verbose); - { //write vtk - const std::string vtkFileName = "afterAdapt" + outname + ".vtk"; - Omega_h::vtk::write_parallel(vtkFileName, &mesh, 2); + { //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; } From 48bfc43c54423fed6e505cea3faa258f31dede7e Mon Sep 17 00:00:00 2001 From: Cameron Smith Date: Fri, 28 Feb 2025 16:11:12 -0500 Subject: [PATCH 24/26] transfer remaining nodal fields --- test/thwaites_adapt.cpp | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/test/thwaites_adapt.cpp b/test/thwaites_adapt.cpp index 0ee3814d..9d3edb02 100644 --- a/test/thwaites_adapt.cpp +++ b/test/thwaites_adapt.cpp @@ -30,6 +30,12 @@ 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; const int numLayers = 11; for(int i=1; i<=numLayers; i++) { std::stringstream ss; From 1f47fb019ebbb6962f098b27740ff3f84ae1cb9b Mon Sep 17 00:00:00 2001 From: Cameron Smith Date: Fri, 28 Feb 2025 20:52:00 -0500 Subject: [PATCH 25/26] transfer more fields --- test/thwaites_adapt.cpp | 3 +++ 1 file changed, 3 insertions(+) diff --git a/test/thwaites_adapt.cpp b/test/thwaites_adapt.cpp index 9d3edb02..36427f01 100644 --- a/test/thwaites_adapt.cpp +++ b/test/thwaites_adapt.cpp @@ -36,6 +36,9 @@ void setupFieldTransfer(AdaptOpts& opts) { 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; From 76382a3f89a8fa532f158f761b93b41b18278d18 Mon Sep 17 00:00:00 2001 From: Cameron Smith Date: Sat, 1 Mar 2025 13:33:21 -0500 Subject: [PATCH 26/26] min and max size input args --- test/thwaites_adapt.cpp | 15 ++++++--------- 1 file changed, 6 insertions(+), 9 deletions(-) diff --git a/test/thwaites_adapt.cpp b/test/thwaites_adapt.cpp index 36427f01..71ecb5fe 100644 --- a/test/thwaites_adapt.cpp +++ b/test/thwaites_adapt.cpp @@ -84,8 +84,8 @@ void printTriCount(Mesh* mesh) { 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 != 4 ) { - fprintf(stderr, "Usage: %s inputMesh.osh outputMeshPrefix adaptRatio\n", argv[0]); + 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(); @@ -94,11 +94,10 @@ int main(int argc, char** argv) { 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]); - auto max_size = Real(16); //roughly maintains the boundary shape - //in the downstream parts of the domain - //where there is significant coarsening - auto min_size = Real(1); //roughly the smallest edge length in the - //source mesh + 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 @@ -107,8 +106,6 @@ int main(int argc, char** argv) { "_maxSz_" + std::to_string(max_size) + "_minSz_" + std::to_string(min_size); - Omega_h::vtk::write_parallel("beforeClassFix_edges.vtk", &mesh, 1); - auto effectiveStrain = getEffectiveStrainRate(mesh); auto recoveredStrain = recoverLinearStrain(mesh,effectiveStrain); mesh.add_tag(VERT, "recoveredStrain", 1, recoveredStrain);