From 0cd999ef63139ce51e893fc41eb2621c8600a871 Mon Sep 17 00:00:00 2001 From: Cameron Smith Date: Thu, 16 Oct 2025 14:41:03 -0400 Subject: [PATCH 01/35] build meshfields as TPL via FetchContent --- CMakeLists.txt | 12 ++++ cmake/GetOrInstallMeshFields.cmake | 96 ++++++++++++++++++++++++++++++ 2 files changed, 108 insertions(+) create mode 100644 cmake/GetOrInstallMeshFields.cmake diff --git a/CMakeLists.txt b/CMakeLists.txt index 55e084356c..aaffbb2d95 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -715,6 +715,18 @@ else() MESSAGE("-- Omega_h lib interface NOT Enabled.") ENDIF() +# MeshFields +OPTION(ENABLE_MESHFIELDS "Flag to turn on the MeshFields field interface for +computing error indicators (via super convergent patch recovery) to +drive mesh adaptation" OFF) +IF (ENABLE_MESHFIELDS) + MESSAGE("-- MeshFields field interface Enabled.") + SET(ALBANY_MESHFIELDS TRUE) + include (GetOrInstallMeshFields) +else() + MESSAGE("-- MeshFields interface NOT Enabled.") +ENDIF() + # add a target to generate API documentation with Doxygen option (DISABLE_DOXYGEN "Whether to disable doxygen" ON) IF( NOT DISABLE_DOXYGEN ) diff --git a/cmake/GetOrInstallMeshFields.cmake b/cmake/GetOrInstallMeshFields.cmake new file mode 100644 index 0000000000..fbb91ed767 --- /dev/null +++ b/cmake/GetOrInstallMeshFields.cmake @@ -0,0 +1,96 @@ +set (tmpStr "Looking for valid MeshFields installation ...") +message (STATUS ${tmpStr}) + +# Get all propreties that cmake supports +if(NOT CMAKE_PROPERTY_LIST) + execute_process(COMMAND cmake --help-property-list OUTPUT_VARIABLE CMAKE_PROPERTY_LIST) + + # Convert command output into a CMake list + string(REGEX REPLACE ";" "\\\\;" CMAKE_PROPERTY_LIST "${CMAKE_PROPERTY_LIST}") + string(REGEX REPLACE "\n" ";" CMAKE_PROPERTY_LIST "${CMAKE_PROPERTY_LIST}") + list(REMOVE_DUPLICATES CMAKE_PROPERTY_LIST) +endif() + +function(print_properties) + message("CMAKE_PROPERTY_LIST = ${CMAKE_PROPERTY_LIST}") +endfunction() + +function(print_target_properties target) + if(NOT TARGET ${target}) + message(STATUS "There is no target named '${target}'") + return() + endif() + + foreach(property ${CMAKE_PROPERTY_LIST}) + string(REPLACE "" "${CMAKE_BUILD_TYPE}" property ${property}) + + # Fix https://stackoverflow.com/questions/32197663/how-can-i-remove-the-the-location-property-may-not-be-read-from-target-error-i + if(property STREQUAL "LOCATION" OR property MATCHES "^LOCATION_" OR property MATCHES "_LOCATION$") + continue() + endif() + + get_property(was_set TARGET ${target} PROPERTY ${property} SET) + if(was_set) + get_target_property(value ${target} ${property}) + message("${target} ${property} = ${value}") + endif() + endforeach() +endfunction() + +find_package(MeshFields CONFIG QUIET + # Avoid all defaults. Only check env/CMake var MeshFields_ROOT + NO_CMAKE_PATH + NO_CMAKE_ENVIRONMENT_PATH + NO_SYSTEM_ENVIRONMENT_PATH + NO_CMAKE_PACKAGE_REGISTRY + NO_CMAKE_SYSTEM_PATH + NO_CMAKE_INSTALL_PREFIX + NO_CMAKE_SYSTEM_PACKAGE_REGISTRY +) + +# If we're building shared libs in Albany, we need MeshFields to be built with shared libs +# TODO: if CMake adds a "DYNAMYC=" to find_package (to specify what variant +# of libs we want), then we can use it, and remove (some of) this snippet +if (MeshFields_FOUND) + message (STATUS "${tmpStr} Found.") + message (" -- MeshFields_DIR: ${MeshFields_DIR}") + message (" -- MeshFields_VERSION: ${MeshFields_VERSION}") + get_target_property(MeshFields_LIBTYPE MeshFields::meshfields TYPE) + if (BUILD_SHARED_LIBS AND NOT MeshFields_LIBTYPE STREQUAL SHARED_LIBRARY) + message (" -> Wrong lib variant: ${MeshFields_LIBTYPE} instead of SHARED_LIBRARY.") + message (" Please, point to a compatible MeshFields installation (or none at all to force a new install)") + message (FATAL_ERROR "Aborting...") + else() + foreach (opt IN ITEMS MeshFields_USE_Kokkos MeshFields_USE_MPI) + if (NOT ${opt}) + message (" -> Wrong configuration: ${opt} is OFF (should be ON)") + message (" Please, point to a compatible MeshFields installation (or none at all to force local install)") + message (FATAL_ERROR "Aborting...") + endif() + endforeach() + endif() +else () + message (STATUS "${tmpStr} NOT Found.") + message (STATUS " -> Downloading and building locally in ${CMAKE_BINARY_DIR}/tpls/meshfields") + + include (FetchContent) + + # Fetch and populate the external project + set (FETCHCONTENT_BASE_DIR ${CMAKE_BINARY_DIR}/tpls/meshfields) + + FetchContent_Declare ( + MeshFields + GIT_REPOSITORY git@github.com:SCOREC/meshFields + GIT_TAG origin/cws/supportFetchContent + OVERRIDE_FIND_PACKAGE + ) + + # Set options for MeshFields before adding the subdirectory + get_target_property(Kokkos_INCLUDE_DIR Kokkos::kokkos INTERFACE_INCLUDE_DIRECTORIES) + string(REPLACE "include/kokkos" "" Kokkos_INSTALL_DIR ${Kokkos_INCLUDE_DIR}) + set(Kokkos_PREFIX ${Kokkos_INSTALL_DIR} PATH "Path to Kokkos install") + + message (STATUS " *** Begin of MeshFields configuration ***") + FetchContent_MakeAvailable (MeshFields) + message (STATUS " *** End of MeshFields configuration ***") +endif() From a04dca03bd4ef57d680b7e25e805ecb54c461195 Mon Sep 17 00:00:00 2001 From: Cameron Smith Date: Thu, 16 Oct 2025 14:47:55 -0400 Subject: [PATCH 02/35] link against meshfields --- src/landIce/CMakeLists.txt | 3 +++ 1 file changed, 3 insertions(+) diff --git a/src/landIce/CMakeLists.txt b/src/landIce/CMakeLists.txt index dccedb0b47..f12eed7c69 100644 --- a/src/landIce/CMakeLists.txt +++ b/src/landIce/CMakeLists.txt @@ -221,6 +221,9 @@ target_link_libraries(landIce PUBLIC albanyLib) if (ALBANY_OMEGAH) target_link_libraries(landIce PUBLIC Omega_h::omega_h) endif() +if (ALBANY_MESHFIELDS) + target_link_libraries(landIce PUBLIC MeshFields::meshfields) +endif() if (ALBANY_SUPPRESS_TRILINOS_WARNINGS) target_include_directories(landIce SYSTEM PUBLIC "${Trilinos_INCLUDE_DIRS};${Trilinos_TPL_INCLUDE_DIRS}") From bf44b690082a35d9970f75e16774b43cb631df0f Mon Sep 17 00:00:00 2001 From: Cameron Smith Date: Thu, 16 Oct 2025 15:01:59 -0400 Subject: [PATCH 03/35] meshfields linking --- src/landIce/CMakeLists.txt | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/landIce/CMakeLists.txt b/src/landIce/CMakeLists.txt index f12eed7c69..a01d3acb30 100644 --- a/src/landIce/CMakeLists.txt +++ b/src/landIce/CMakeLists.txt @@ -222,7 +222,7 @@ if (ALBANY_OMEGAH) target_link_libraries(landIce PUBLIC Omega_h::omega_h) endif() if (ALBANY_MESHFIELDS) - target_link_libraries(landIce PUBLIC MeshFields::meshfields) + target_link_libraries(landIce PUBLIC meshfields::meshfields) endif() if (ALBANY_SUPPRESS_TRILINOS_WARNINGS) target_include_directories(landIce SYSTEM PUBLIC From 24dbecefd5b02d6696434189a93673690bc478ab Mon Sep 17 00:00:00 2001 From: Cameron Smith Date: Thu, 16 Oct 2025 15:03:38 -0400 Subject: [PATCH 04/35] no scope? --- src/landIce/CMakeLists.txt | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/landIce/CMakeLists.txt b/src/landIce/CMakeLists.txt index a01d3acb30..b799f6b5de 100644 --- a/src/landIce/CMakeLists.txt +++ b/src/landIce/CMakeLists.txt @@ -222,7 +222,7 @@ if (ALBANY_OMEGAH) target_link_libraries(landIce PUBLIC Omega_h::omega_h) endif() if (ALBANY_MESHFIELDS) - target_link_libraries(landIce PUBLIC meshfields::meshfields) + target_link_libraries(landIce PUBLIC meshfields) endif() if (ALBANY_SUPPRESS_TRILINOS_WARNINGS) target_include_directories(landIce SYSTEM PUBLIC From 6fdb6979f99c1f0d35edc2f0f1c0c538b66342ea Mon Sep 17 00:00:00 2001 From: Cameron Smith Date: Thu, 16 Oct 2025 16:05:32 -0400 Subject: [PATCH 05/35] link albanyLib against meshfields --- src/CMakeLists.txt | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 979984c5fe..eeda3960c8 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -526,6 +526,10 @@ if (ALBANY_OMEGAH) target_link_libraries(albanyLib Omega_h::omega_h) endif() +if (ALBANY_MESHFIELDS) + target_link_libraries(albanyLib meshfields) +endif() + install(TARGETS albanyLib EXPORT albany-export LIBRARY DESTINATION ${CMAKE_INSTALL_LIBDIR} ARCHIVE DESTINATION ${CMAKE_INSTALL_LIBDIR} From c6bdc4bb1f60972c89e0f0400024967d1e7b36cf Mon Sep 17 00:00:00 2001 From: Cameron Smith Date: Thu, 16 Oct 2025 16:15:41 -0400 Subject: [PATCH 06/35] preprocessor symbol for meshfields --- src/Albany_config.h.in | 3 +++ 1 file changed, 3 insertions(+) diff --git a/src/Albany_config.h.in b/src/Albany_config.h.in index bc7b7a1582..d0167a3d19 100644 --- a/src/Albany_config.h.in +++ b/src/Albany_config.h.in @@ -30,6 +30,9 @@ // Whether the Omega_h interface is enabled #cmakedefine ALBANY_OMEGAH +// Whether the MeshFields interface is enabled +#cmakedefine ALBANY_MESHFIELDS + // Whether this is a nightly test build #cmakedefine ALBANY_NIGHTLY_TEST From 054ee66bfd81bbc391d0d149987f411462532abc Mon Sep 17 00:00:00 2001 From: Cameron Smith Date: Thu, 16 Oct 2025 16:15:54 -0400 Subject: [PATCH 07/35] wrap in preprocessor macro --- src/disc/omegah/Albany_OmegahDiscretization.cpp | 9 +++++++++ 1 file changed, 9 insertions(+) diff --git a/src/disc/omegah/Albany_OmegahDiscretization.cpp b/src/disc/omegah/Albany_OmegahDiscretization.cpp index 6917ec9b1f..fe1361e737 100644 --- a/src/disc/omegah/Albany_OmegahDiscretization.cpp +++ b/src/disc/omegah/Albany_OmegahDiscretization.cpp @@ -8,6 +8,11 @@ #include "Omega_h_array_ops.hpp" #include // for Omega_h::binary::write +#ifdef ALBANY_MESHFIELDS +#include +#include +#endif + #include namespace debug { @@ -435,6 +440,10 @@ checkForAdaptation (const Teuchos::RCP& solution , } } +#ifdef ALBANY_MESHFIELDS + MeshField::OmegahMeshField omf(*mesh); +#endif + return adapt_data; } From 32ae5cdda7252147946e063e758145905df7f35e Mon Sep 17 00:00:00 2001 From: Cameron Smith Date: Thu, 16 Oct 2025 16:23:26 -0400 Subject: [PATCH 08/35] switch back to main --- cmake/GetOrInstallMeshFields.cmake | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/cmake/GetOrInstallMeshFields.cmake b/cmake/GetOrInstallMeshFields.cmake index fbb91ed767..5172459e8b 100644 --- a/cmake/GetOrInstallMeshFields.cmake +++ b/cmake/GetOrInstallMeshFields.cmake @@ -81,7 +81,7 @@ else () FetchContent_Declare ( MeshFields GIT_REPOSITORY git@github.com:SCOREC/meshFields - GIT_TAG origin/cws/supportFetchContent + GIT_TAG origin/main OVERRIDE_FIND_PACKAGE ) From 6c2e22ee9fa46424f70a717ec5c7f1cc821f80ad Mon Sep 17 00:00:00 2001 From: Cameron Smith Date: Tue, 28 Oct 2025 21:53:14 -0400 Subject: [PATCH 09/35] start adding spr calls porting from meshfields/test/testSprThwaitesAdapt.cpp --- .../omegah/Albany_OmegahDiscretization.cpp | 154 ++++++++++++++---- 1 file changed, 119 insertions(+), 35 deletions(-) diff --git a/src/disc/omegah/Albany_OmegahDiscretization.cpp b/src/disc/omegah/Albany_OmegahDiscretization.cpp index fe1361e737..f79225ea15 100644 --- a/src/disc/omegah/Albany_OmegahDiscretization.cpp +++ b/src/disc/omegah/Albany_OmegahDiscretization.cpp @@ -7,14 +7,54 @@ #include "Omega_h_adapt.hpp" #include "Omega_h_array_ops.hpp" #include // for Omega_h::binary::write +#include "Omega_h_recover.hpp" //project_by_fit #ifdef ALBANY_MESHFIELDS #include #include +#include #endif #include +//FIXME! +using ExecutionSpace = Kokkos::DefaultExecutionSpace; +using MemorySpace = Kokkos::DefaultExecutionSpace::memory_space; + +namespace { + constexpr bool isMeshfieldsEnabled() { +#ifdef ALBANY_MESHFIELDS + return true; +#else + return false; +#endif + } + + Omega_h::Reals getEffectiveStrainRate(Omega_h::Mesh &mesh) { + return mesh.get_array(2, "solution_grad_norm"); + } + + Omega_h::Reals recoverLinearStrain(Omega_h::Mesh &mesh, Omega_h::Reals effectiveStrain) { + return Omega_h::project_by_fit(&mesh, effectiveStrain); + } + + template + void setFieldAtVertices(Omega_h::Mesh &mesh, Omega_h::Reals recoveredStrain, + ShapeField field) { + auto setFieldAtVertices = KOKKOS_LAMBDA(const int &vtx) { + field(vtx, 0, 0, MeshField::Vertex) = recoveredStrain[vtx]; + }; + MeshField::parallel_for(ExecutionSpace(), {0}, {mesh.nverts()}, + setFieldAtVertices, "setFieldAtVertices"); + } + + void printTriCount(Omega_h::Mesh &mesh, std::string_view prefix) { + const auto nTri = mesh.nglobal_ents(2); + if (!mesh.comm()->rank()) + std::cout << prefix << " nTri: " << nTri << "\n"; + } +} + namespace debug { template void printTagInfo(Omega_h::Mesh mesh, std::ostringstream& oss, int dim, int tag, std::string type) { @@ -390,19 +430,14 @@ checkForAdaptation (const Teuchos::RCP& solution , const Teuchos::RCP& dxdp) { auto adapt_data = Teuchos::rcp(new AdaptationData()); - - // Only do adaptation for simple 1d problems auto mesh = m_mesh_struct->getOmegahMesh(); - if (mesh->dim() != 1) { - std::cout << "NOT a 1D Omega_h mesh...we will not adapt.\n"; - return adapt_data; - } auto& adapt_params = m_disc_params->sublist("Mesh Adaptivity"); auto adapt_type = adapt_params.get("Type","None"); if (adapt_type=="None") { return adapt_data; } - TEUCHOS_TEST_FOR_EXCEPTION (adapt_type!="Minimally-Oscillatory", std::runtime_error, + + TEUCHOS_TEST_FOR_EXCEPTION (adapt_type!="Minimally-Oscillatory", std::runtime_error, //FIXME - need an SPR option "Error! Adaptation type '" << adapt_type << "' not supported.\n" " - valid choices: None, Minimally-Oscillatory\n"); @@ -410,42 +445,91 @@ checkForAdaptation (const Teuchos::RCP& solution , "Error! the dxdp Thyra_MultiVector is expected to be null\n"); if(solution_dot != Teuchos::null and solution_dotdot != Teuchos::null) { - writeSolutionToMeshDatabase(*solution, dxdp, *solution_dot, *solution_dotdot, false); + writeSolutionToMeshDatabase(*solution, dxdp, *solution_dot, *solution_dotdot, false); } else if(solution_dot != Teuchos::null) { - writeSolutionToMeshDatabase(*solution, dxdp, *solution_dot, false); + writeSolutionToMeshDatabase(*solution, dxdp, *solution_dot, false); } else { - writeSolutionToMeshDatabase(*solution, dxdp, false); + writeSolutionToMeshDatabase(*solution, dxdp, false); } - double tol = adapt_params.get("Max Hessian"); - auto data = getLocalData(solution); - // Simple check: refine if a proxy of the hessian of x is larger than a tolerance - // TODO: replace with - // 1. if |C_i| > threshold, mark for refinement the whole mesh - // 2. Interpolate solution (and all elem/node fields if possible, but not necessary for adv-diff example) - int num_nodes = data.size(); - adapt_data->x = solution; - adapt_data->x_dot = solution_dot; - adapt_data->x_dotdot = solution_dotdot; - adapt_data->dxdp = dxdp; - for (int i=1; itol and grad_prev*grad_next<0) { - adapt_data->type = AdaptationType::Topology; - break; + if (mesh->dim() == 1) { + double tol = adapt_params.get("Max Hessian"); + auto data = getLocalData(solution); + // Simple check: refine if a proxy of the hessian of x is larger than a tolerance + // TODO: replace with + // 1. if |C_i| > threshold, mark for refinement the whole mesh + // 2. Interpolate solution (and all elem/node fields if possible, but not necessary for adv-diff example) + int num_nodes = data.size(); + adapt_data->x = solution; + adapt_data->x_dot = solution_dot; + adapt_data->x_dotdot = solution_dotdot; + adapt_data->dxdp = dxdp; + for (int i=1; itol and grad_prev*grad_next<0) { + adapt_data->type = AdaptationType::Topology; + break; + } + } + return adapt_data; + } else if (mesh->dim() == 2) { + if (!isMeshfieldsEnabled()) { + std::cout << "2D Omega_h mesh adaptation requires Meshfields " + << "(configure Albany with ENABLE_MESHFIELDS=ON to enable it) " + << "... we will not adapt.\n"; + return adapt_data; } - } -#ifdef ALBANY_MESHFIELDS - MeshField::OmegahMeshField omf(*mesh); -#endif + #ifdef ALBANY_MESHFIELDS + auto effectiveStrain = getEffectiveStrainRate(*mesh); + auto recoveredStrain = recoverLinearStrain(*mesh, effectiveStrain); + mesh->add_tag(Omega_h::VERT, "recoveredStrain", 1, recoveredStrain, + false, Omega_h::ArrayType::VectorND); + + const auto MeshDim = 2; + const auto ShapeOrder = 1; + MeshField::OmegahMeshField omf(*mesh); + 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 adaptRatio = 0.1; + auto estimation = + Estimation(*mesh, effectiveStrain, recoveredStrainField, adaptRatio); //FIXME - move to meshfields namespace + + const auto tgtLength = getSprSizeField(estimation, omf, coordFe); + Omega_h::Write tgtLength_oh(tgtLength); + mesh->add_tag(Omega_h::VERT, "tgtLength", 1, tgtLength_oh, false, + Omega_h::ArrayType::VectorND); + + { // write vtk + const auto outname = std::string("foo"); + 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); + } - return adapt_data; + printTriCount(*mesh, "beforeAdapt"); + #endif + return adapt_data; + } else { //meshdim != 1 && meshdim != 2 + std::cout << "Only 1D and 2D (with Meshfields enabled) Omega_h mesh " + << "adaptation is supported ... we will not adapt.\n"; + return adapt_data; + } } void OmegahDiscretization:: From 232a0e6b5d87ee666bfad6e803e2fdff6390a196 Mon Sep 17 00:00:00 2001 From: Cameron Smith Date: Tue, 28 Oct 2025 22:18:09 -0400 Subject: [PATCH 10/35] use namespace --- src/disc/omegah/Albany_OmegahDiscretization.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/disc/omegah/Albany_OmegahDiscretization.cpp b/src/disc/omegah/Albany_OmegahDiscretization.cpp index f79225ea15..5824d52f7f 100644 --- a/src/disc/omegah/Albany_OmegahDiscretization.cpp +++ b/src/disc/omegah/Albany_OmegahDiscretization.cpp @@ -505,9 +505,9 @@ checkForAdaptation (const Teuchos::RCP& solution , const auto adaptRatio = 0.1; auto estimation = - Estimation(*mesh, effectiveStrain, recoveredStrainField, adaptRatio); //FIXME - move to meshfields namespace + MeshField::SPR::Estimation(*mesh, effectiveStrain, recoveredStrainField, adaptRatio); - const auto tgtLength = getSprSizeField(estimation, omf, coordFe); + const auto tgtLength = MeshField::SPR::getSprSizeField(estimation, omf, coordFe); Omega_h::Write tgtLength_oh(tgtLength); mesh->add_tag(Omega_h::VERT, "tgtLength", 1, tgtLength_oh, false, Omega_h::ArrayType::VectorND); From ca31e0d87a1ecb056009922f35620c002fafbc0d Mon Sep 17 00:00:00 2001 From: Cameron Smith Date: Wed, 29 Oct 2025 21:58:00 -0400 Subject: [PATCH 11/35] spr: adapt calls --- .../omegah/Albany_OmegahDiscretization.cpp | 61 ++++++++++++------- 1 file changed, 40 insertions(+), 21 deletions(-) diff --git a/src/disc/omegah/Albany_OmegahDiscretization.cpp b/src/disc/omegah/Albany_OmegahDiscretization.cpp index 5824d52f7f..43fa215648 100644 --- a/src/disc/omegah/Albany_OmegahDiscretization.cpp +++ b/src/disc/omegah/Albany_OmegahDiscretization.cpp @@ -5,6 +5,7 @@ #include "OmegahConnManager.hpp" #include "Omega_h_adapt.hpp" +#include "Omega_h_metric.hpp" // isos_from_lengths, clamp_metrics #include "Omega_h_array_ops.hpp" #include // for Omega_h::binary::write #include "Omega_h_recover.hpp" //project_by_fit @@ -523,6 +524,7 @@ checkForAdaptation (const Teuchos::RCP& solution , printTriCount(*mesh, "beforeAdapt"); #endif + adapt_data->type = AdaptationType::Topology; return adapt_data; } else { //meshdim != 1 && meshdim != 2 @@ -550,30 +552,47 @@ adapt (const Teuchos::RCP& adaptData) std::string beforeAdaptName = "before_adapt" + std::to_string(adaptCount) + ".vtk"; Omega_h::vtk::write_parallel(beforeAdaptName, ohMesh.get()); - // Note: the code below is hard-coding a simple adaptation for a 1d mesh, - // where the number of elements is doubled. - auto nelems = ohMesh->nglobal_ents(ohMesh->dim()); - const auto desired_nelems = nelems*2; - - Omega_h::AdaptOpts opts(&(*ohMesh)); - opts.xfer_opts.type_map[solution_dof_name()] = OMEGA_H_LINEAR_INTERP; - opts.xfer_opts.type_map[std::string(solution_dof_name())+"_dot"] = OMEGA_H_LINEAR_INTERP; - while (double(nelems) < desired_nelems) { - if (!ohMesh->has_tag(0, "metric")) { - std::cout << "mesh had no metric, adding implied and adapting to it\n"; - Omega_h::add_implied_metric_tag(ohMesh.get()); + if (ohMesh->dim() == 1) { + // Note: the code below is hard-coding a simple adaptation for a 1d mesh, + // where the number of elements is doubled. + auto nelems = ohMesh->nglobal_ents(ohMesh->dim()); + const auto desired_nelems = nelems*2; + + Omega_h::AdaptOpts opts(&(*ohMesh)); + opts.xfer_opts.type_map[solution_dof_name()] = OMEGA_H_LINEAR_INTERP; + opts.xfer_opts.type_map[std::string(solution_dof_name())+"_dot"] = OMEGA_H_LINEAR_INTERP; + while (double(nelems) < desired_nelems) { + if (!ohMesh->has_tag(0, "metric")) { + std::cout << "mesh had no metric, adding implied and adapting to it\n"; + Omega_h::add_implied_metric_tag(ohMesh.get()); + Omega_h::adapt(ohMesh.get(), opts); + nelems = ohMesh->nglobal_ents(ohMesh->dim()); + } + auto metrics = ohMesh->get_array(0, "metric"); + metrics = Omega_h::multiply_each_by(metrics, 1.2); + auto const metric_ncomps = + Omega_h::divide_no_remainder(metrics.size(), ohMesh->nverts()); + ohMesh->add_tag(0, "metric", metric_ncomps, metrics); + std::cout << "adapting to scaled metric\n"; Omega_h::adapt(ohMesh.get(), opts); nelems = ohMesh->nglobal_ents(ohMesh->dim()); + std::cout << "mesh now has " << nelems << " total elements\n"; } - auto metrics = ohMesh->get_array(0, "metric"); - metrics = Omega_h::multiply_each_by(metrics, 1.2); - auto const metric_ncomps = - Omega_h::divide_no_remainder(metrics.size(), ohMesh->nverts()); - ohMesh->add_tag(0, "metric", metric_ncomps, metrics); - std::cout << "adapting to scaled metric\n"; - Omega_h::adapt(ohMesh.get(), opts); - nelems = ohMesh->nglobal_ents(ohMesh->dim()); - std::cout << "mesh now has " << nelems << " total elements\n"; + } else if (ohMesh->dim() == 2) { + // adapt + Omega_h::AdaptOpts opts(&(*ohMesh)); + opts.xfer_opts.type_map[solution_dof_name()] = OMEGA_H_LINEAR_INTERP; + opts.xfer_opts.type_map[std::string(solution_dof_name())+"_dot"] = OMEGA_H_LINEAR_INTERP; + + auto verbose = false; + const auto tgtLength_oh = ohMesh->get_array(Omega_h::VERT, "tgtLength"); + const auto isos = Omega_h::isos_from_lengths(tgtLength_oh); + const auto min_size = 0.1; + const auto max_size = 1.0; + auto metric = Omega_h::clamp_metrics(ohMesh->nverts(), isos, min_size, max_size); + Omega_h::grade_fix_adapt(&(*ohMesh), opts, metric, verbose); + + printTriCount(*ohMesh, "afterAdapt"); } std::string afterAdaptName = "after_adapt" + std::to_string(adaptCount) + ".vtk"; From 6cd9770e2b7e7aade088a3e89425c609fa01f1ae Mon Sep 17 00:00:00 2001 From: Cameron Smith Date: Thu, 30 Oct 2025 13:13:46 -0400 Subject: [PATCH 12/35] check for meshfields --- src/disc/omegah/Albany_OmegahDiscretization.cpp | 13 ++++++------- 1 file changed, 6 insertions(+), 7 deletions(-) diff --git a/src/disc/omegah/Albany_OmegahDiscretization.cpp b/src/disc/omegah/Albany_OmegahDiscretization.cpp index 43fa215648..74c8b904a8 100644 --- a/src/disc/omegah/Albany_OmegahDiscretization.cpp +++ b/src/disc/omegah/Albany_OmegahDiscretization.cpp @@ -479,8 +479,8 @@ checkForAdaptation (const Teuchos::RCP& solution , return adapt_data; } else if (mesh->dim() == 2) { if (!isMeshfieldsEnabled()) { - std::cout << "2D Omega_h mesh adaptation requires Meshfields " - << "(configure Albany with ENABLE_MESHFIELDS=ON to enable it) " + std::cout << "Warning: 2D Omega_h mesh adaptation requires Meshfields. " + << "Configure Albany with ENABLE_MESHFIELDS=ON to enable it. " << "... we will not adapt.\n"; return adapt_data; } @@ -514,11 +514,10 @@ checkForAdaptation (const Teuchos::RCP& solution , Omega_h::ArrayType::VectorND); { // write vtk - const auto outname = std::string("foo"); - const std::string vtkFileName = "beforeAdapt" + outname + ".vtk"; + const auto outname = std::string("beforeAdapt2d"); + const std::string vtkFileName = outname + ".vtk"; Omega_h::vtk::write_parallel(vtkFileName, &(*mesh), 2); - const std::string vtkFileName_edges = - "beforeAdapt" + outname + "_edges.vtk"; + const std::string vtkFileName_edges = outname + "_edges.vtk"; Omega_h::vtk::write_parallel(vtkFileName_edges, &(*mesh), 1); } @@ -578,7 +577,7 @@ adapt (const Teuchos::RCP& adaptData) nelems = ohMesh->nglobal_ents(ohMesh->dim()); std::cout << "mesh now has " << nelems << " total elements\n"; } - } else if (ohMesh->dim() == 2) { + } else if (ohMesh->dim() == 2 && isMeshfieldsEnabled()) { // adapt Omega_h::AdaptOpts opts(&(*ohMesh)); opts.xfer_opts.type_map[solution_dof_name()] = OMEGA_H_LINEAR_INTERP; From 0c53c3608f571085033aecde787a702eb0550a33 Mon Sep 17 00:00:00 2001 From: Cameron Smith Date: Thu, 30 Oct 2025 13:22:54 -0400 Subject: [PATCH 13/35] 2d omegah adapt test --- tests/corePDEs/AdvDiff/CMakeLists.txt | 8 +- .../AdvDiff/input_2d_adapt_omegah.yaml | 130 ++++++++++++++++++ 2 files changed, 137 insertions(+), 1 deletion(-) create mode 100644 tests/corePDEs/AdvDiff/input_2d_adapt_omegah.yaml diff --git a/tests/corePDEs/AdvDiff/CMakeLists.txt b/tests/corePDEs/AdvDiff/CMakeLists.txt index 064d55662f..9a1c13cfa7 100644 --- a/tests/corePDEs/AdvDiff/CMakeLists.txt +++ b/tests/corePDEs/AdvDiff/CMakeLists.txt @@ -40,9 +40,15 @@ IF(NOT ALBANY_PARALLEL_ONLY AND ALBANY_IFPACK2) # 2d test configure_file(${CMAKE_CURRENT_SOURCE_DIR}/input_2d_omegah.yaml ${CMAKE_CURRENT_BINARY_DIR}/input_2d_omegah.yaml COPYONLY) - add_test(${testName}_2d_omegah ${SerialAlbany.exe} input_2d_stk.yaml) + add_test(${testName}_2d_omegah ${SerialAlbany.exe} input_2d_omegah.yaml) set_tests_properties(${testName}_2d_omegah PROPERTIES LABELS "Demo;Forward;Serial") + # 2d test with adaptivity + configure_file(${CMAKE_CURRENT_SOURCE_DIR}/input_2d_adapt_omegah.yaml + ${CMAKE_CURRENT_BINARY_DIR}/input_2d_adapt_omegah.yaml COPYONLY) + add_test(${testName}_2d_adapt_omegah ${SerialAlbany.exe} input_2d_adapt_omegah.yaml) + set_tests_properties(${testName}_2d_adapt_omegah PROPERTIES LABELS "Demo;Forward;Serial;Adapt") + ENDIF() ENDIF () diff --git a/tests/corePDEs/AdvDiff/input_2d_adapt_omegah.yaml b/tests/corePDEs/AdvDiff/input_2d_adapt_omegah.yaml new file mode 100644 index 0000000000..f2731dfb2e --- /dev/null +++ b/tests/corePDEs/AdvDiff/input_2d_adapt_omegah.yaml @@ -0,0 +1,130 @@ +%YAML 1.1 +--- +ANONYMOUS: + Debug Output: + Write Solution to MatrixMarket: 0 + Report Timers: false + Problem: + Name: AdvDiff 2D + Solution Method: Transient + Number of PDE Equations: 1 + Dirichlet BCs: + DBC on NS NodeSet0 for DOF U0: 0.0 + DBC on NS NodeSet2 for DOF U0: 0.0 + Initial Condition: + Function: Circle + Options: + Advection Coefficient: [1.0,1.0] + Viscosity mu: 0.01 + Parameters: + Number Of Parameters: 0 + Response Functions: + Number Of Responses: 3 + Response 0: + Type: Scalar Response + Name: Solution Max Value + Response 1: + Type: Scalar Response + Name: Solution Min Value + Response 2: + Type: Scalar Response + Name: Solution Average + Discretization: + Method: Omegah + Mesh Creation Method: Box2D + Mesh Adaptivity: + Type: Minimally-Oscillatory + Max Hessian: 10.0 + Refining Factor: 2 + Piro: + Tempus: + Integrator Name: Tempus Integrator + Tempus Integrator: + Integrator Type: Integrator Basic + Stepper Name: Tempus Stepper + Solution History: + Storage Type: Unlimited + Storage Limit: 20 + Time Step Control: + Initial Time: 0.0 + Initial Time Index: 0 + Initial Time Step: 0.05 + Final Time: 1.0 + Maximum Absolute Error: 1.0e-08 + Maximum Relative Error: 1.0e-08 + Maximum Number of Stepper Failures: 10 + Maximum Number of Consecutive Stepper Failures: 5 + Tempus Stepper: + Stepper Type: Backward Euler + Solver Name: Demo Solver + Demo Solver: + NOX: + Direction: + Method: Newton + Newton: + Forcing Term Method: Constant + Rescue Bad Newton Solve: true + Linear Solver: + Tolerance: 1.0e-05 + Line Search: + Full Step: + Full Step: 1.0 + Method: Full Step + Nonlinear Solver: Line Search Based + Printing: + Output Precision: 3 + Output Processor: 0 + Output Information: + Error: true + Warning: true + Outer Iteration: false + Parameters: true + Details: false + Linear Solver Details: true + Stepper Iteration: true + Stepper Details: true + Stepper Parameters: true + Solver Options: + Status Test Check Type: Minimal + Status Tests: + Test Type: Combo + Combo Type: OR + Number of Tests: 2 + Test 0: + Test Type: NormF + Tolerance: 1.0e-08 + Test 1: + Test Type: MaxIters + Maximum Iterations: 10 + Stratimikos: + Linear Solver Type: Belos + Linear Solver Types: + Belos: + Solver Type: Block GMRES + Solver Types: + Block GMRES: + Convergence Tolerance: 1.0e-05 + Output Frequency: 10 + Output Style: 1 + Verbosity: 0 + Maximum Iterations: 10 + Block Size: 1 + Num Blocks: 100 + Flexible Gmres: false + Preconditioner Type: Ifpack2 + Preconditioner Types: + Ifpack2: + Prec Type: ILUT + Overlap: 1 + Ifpack2 Settings: + 'fact: ilut level-of-fill': 1.0e+00 + Regression For Response 0: + Test Value: 6.825290995201e-02 + Relative Tolerance: 1.0e-04 + Regression For Response 1: + Test Value: -8.259050370473e-04 + Relative Tolerance: 1.0e-04 + Regression For Response 2: + Test Value: 2.170508705497e-03 + Relative Tolerance: 1.0e-04 +... From 44f703c6c0983712d81c998917513e840aadd358 Mon Sep 17 00:00:00 2001 From: Cameron Smith Date: Thu, 30 Oct 2025 13:23:13 -0400 Subject: [PATCH 14/35] consistent names --- tests/corePDEs/AdvDiff/CMakeLists.txt | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/corePDEs/AdvDiff/CMakeLists.txt b/tests/corePDEs/AdvDiff/CMakeLists.txt index 9a1c13cfa7..de15ef4323 100644 --- a/tests/corePDEs/AdvDiff/CMakeLists.txt +++ b/tests/corePDEs/AdvDiff/CMakeLists.txt @@ -35,7 +35,7 @@ IF(NOT ALBANY_PARALLEL_ONLY AND ALBANY_IFPACK2) configure_file(${CMAKE_CURRENT_SOURCE_DIR}/input_1d_adapt_omegah.yaml ${CMAKE_CURRENT_BINARY_DIR}/input_1d_adapt_omegah.yaml COPYONLY) add_test(${testName}_1d_adapt_Omega_h ${SerialAlbany.exe} input_1d_adapt_omegah.yaml) - set_tests_properties(${testName}_1d_adapt_Omega_h PROPERTIES LABELS "Demo;Forward;Serial;Adapt") + set_tests_properties(${testName}_1d_adapt_omega_h PROPERTIES LABELS "Demo;Forward;Serial;Adapt") # 2d test configure_file(${CMAKE_CURRENT_SOURCE_DIR}/input_2d_omegah.yaml From 4d7f3b8d45ce067159aeea40c61d7995d4197c18 Mon Sep 17 00:00:00 2001 From: Cameron Smith Date: Thu, 30 Oct 2025 13:24:07 -0400 Subject: [PATCH 15/35] missed an Omegah --- tests/corePDEs/AdvDiff/CMakeLists.txt | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/corePDEs/AdvDiff/CMakeLists.txt b/tests/corePDEs/AdvDiff/CMakeLists.txt index de15ef4323..e11785e447 100644 --- a/tests/corePDEs/AdvDiff/CMakeLists.txt +++ b/tests/corePDEs/AdvDiff/CMakeLists.txt @@ -34,7 +34,7 @@ IF(NOT ALBANY_PARALLEL_ONLY AND ALBANY_IFPACK2) # 1d test with adaptivity configure_file(${CMAKE_CURRENT_SOURCE_DIR}/input_1d_adapt_omegah.yaml ${CMAKE_CURRENT_BINARY_DIR}/input_1d_adapt_omegah.yaml COPYONLY) - add_test(${testName}_1d_adapt_Omega_h ${SerialAlbany.exe} input_1d_adapt_omegah.yaml) + add_test(${testName}_1d_adapt_omega_h ${SerialAlbany.exe} input_1d_adapt_omegah.yaml) set_tests_properties(${testName}_1d_adapt_omega_h PROPERTIES LABELS "Demo;Forward;Serial;Adapt") # 2d test From 16a1a9f8e328c4ca80ffb05c14f9483ba88f1ec9 Mon Sep 17 00:00:00 2001 From: Cameron Smith Date: Thu, 30 Oct 2025 13:25:42 -0400 Subject: [PATCH 16/35] 2d omegah adapt requires meshfields for spr --- tests/corePDEs/AdvDiff/CMakeLists.txt | 12 +++++++----- 1 file changed, 7 insertions(+), 5 deletions(-) diff --git a/tests/corePDEs/AdvDiff/CMakeLists.txt b/tests/corePDEs/AdvDiff/CMakeLists.txt index e11785e447..39ce1f716c 100644 --- a/tests/corePDEs/AdvDiff/CMakeLists.txt +++ b/tests/corePDEs/AdvDiff/CMakeLists.txt @@ -43,11 +43,13 @@ IF(NOT ALBANY_PARALLEL_ONLY AND ALBANY_IFPACK2) add_test(${testName}_2d_omegah ${SerialAlbany.exe} input_2d_omegah.yaml) set_tests_properties(${testName}_2d_omegah PROPERTIES LABELS "Demo;Forward;Serial") - # 2d test with adaptivity - configure_file(${CMAKE_CURRENT_SOURCE_DIR}/input_2d_adapt_omegah.yaml - ${CMAKE_CURRENT_BINARY_DIR}/input_2d_adapt_omegah.yaml COPYONLY) - add_test(${testName}_2d_adapt_omegah ${SerialAlbany.exe} input_2d_adapt_omegah.yaml) - set_tests_properties(${testName}_2d_adapt_omegah PROPERTIES LABELS "Demo;Forward;Serial;Adapt") + if (ENABLE_MESHFIELDS) + # 2d test with adaptivity + configure_file(${CMAKE_CURRENT_SOURCE_DIR}/input_2d_adapt_omegah.yaml + ${CMAKE_CURRENT_BINARY_DIR}/input_2d_adapt_omegah.yaml COPYONLY) + add_test(${testName}_2d_adapt_omegah ${SerialAlbany.exe} input_2d_adapt_omegah.yaml) + set_tests_properties(${testName}_2d_adapt_omegah PROPERTIES LABELS "Demo;Forward;Serial;Adapt") + ENDIF() ENDIF() From c3c56ea0c8eaa12fabe3d1c084dbc9811a400061 Mon Sep 17 00:00:00 2001 From: Cameron Smith Date: Thu, 30 Oct 2025 13:41:28 -0400 Subject: [PATCH 17/35] need to specify initial elm count --- tests/corePDEs/AdvDiff/input_2d_adapt_omegah.yaml | 1 + 1 file changed, 1 insertion(+) diff --git a/tests/corePDEs/AdvDiff/input_2d_adapt_omegah.yaml b/tests/corePDEs/AdvDiff/input_2d_adapt_omegah.yaml index f2731dfb2e..1f74bc6403 100644 --- a/tests/corePDEs/AdvDiff/input_2d_adapt_omegah.yaml +++ b/tests/corePDEs/AdvDiff/input_2d_adapt_omegah.yaml @@ -30,6 +30,7 @@ ANONYMOUS: Type: Scalar Response Name: Solution Average Discretization: + Number of Elements: [10,10] Method: Omegah Mesh Creation Method: Box2D Mesh Adaptivity: From b860884a6acb7753fec863c036f67dd54f005aef Mon Sep 17 00:00:00 2001 From: Cameron Smith Date: Thu, 30 Oct 2025 13:46:27 -0400 Subject: [PATCH 18/35] 2d supported --- src/disc/omegah/Albany_OmegahDiscretization.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/disc/omegah/Albany_OmegahDiscretization.cpp b/src/disc/omegah/Albany_OmegahDiscretization.cpp index 74c8b904a8..0daca621bf 100644 --- a/src/disc/omegah/Albany_OmegahDiscretization.cpp +++ b/src/disc/omegah/Albany_OmegahDiscretization.cpp @@ -545,8 +545,8 @@ adapt (const Teuchos::RCP& adaptData) auto ohMesh = m_mesh_struct->getOmegahMesh(); TEUCHOS_TEST_FOR_EXCEPTION (adaptData->type!=AdaptationType::Topology, std::runtime_error, "Error! Adaptation type not supported. Only 'None' and 'Topology' are currently supported.\n"); - TEUCHOS_TEST_FOR_EXCEPTION (ohMesh->dim()!=1, std::runtime_error, - "Error! Adaptation not supported for this mesh. We only implemented a simple 1d case.\n"); + TEUCHOS_TEST_FOR_EXCEPTION (ohMesh->dim()!=1 && ohMesh->dim()!=2, std::runtime_error, + "Error! Adaptation not supported for this mesh. We only implemented simple 1d and 2d cases.\n"); std::string beforeAdaptName = "before_adapt" + std::to_string(adaptCount) + ".vtk"; Omega_h::vtk::write_parallel(beforeAdaptName, ohMesh.get()); From 1ff67efb4be6723851966bda7b1b6abd8fbf4766 Mon Sep 17 00:00:00 2001 From: Cameron Smith Date: Thu, 30 Oct 2025 13:47:30 -0400 Subject: [PATCH 19/35] match function name --- src/disc/omegah/Albany_OmegahGenericMesh.cpp | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/disc/omegah/Albany_OmegahGenericMesh.cpp b/src/disc/omegah/Albany_OmegahGenericMesh.cpp index e27c2fb247..6bd79259c6 100644 --- a/src/disc/omegah/Albany_OmegahGenericMesh.cpp +++ b/src/disc/omegah/Albany_OmegahGenericMesh.cpp @@ -147,19 +147,19 @@ mark_part_entities (const std::string& name, const bool markDownward) { TEUCHOS_TEST_FOR_EXCEPTION (m_part_topo.find(name)==m_part_topo.end(), std::runtime_error, - "[OmegahGenericMesh::set_part_entities] Error! Part not found.\n" + "[OmegahGenericMesh::mark_part_entities] Error! Part not found.\n" " - part name: " << name << "\n"); auto dim = topo_dim(m_part_topo[name]); TEUCHOS_TEST_FOR_EXCEPTION (is_entity_in_part.size()!=m_mesh->nents(dim), std::logic_error, - "[OmegahGenericMesh::set_part_entities] Error! Input array has the wrong dimensions.\n" + "[OmegahGenericMesh::mark_part_entities] Error! Input array has the wrong dimensions.\n" " - part name: " << name << "\n" " - part dim : " << dim << "\n" " - num ents : " << m_mesh->nents(dim) << "\n" " - array dim: " << is_entity_in_part.size() << "\n"); TEUCHOS_TEST_FOR_EXCEPTION (m_mesh->has_tag(dim,name), std::runtime_error, - "[OmegahGenericMesh::set_part_entities] Error! A tag with this name was already set.\n" + "[OmegahGenericMesh::mark_part_entities] Error! A tag with this name was already set.\n" " - part name: " << name << "\n" " - part dim : " << dim << "\n"); @@ -167,7 +167,7 @@ mark_part_entities (const std::string& name, if (markDownward) { TEUCHOS_TEST_FOR_EXCEPTION (dim==0, std::logic_error, - "[OmegahGenericMesh::set_part_entities] Error! Cannot mark downward if the part dimension is 0.\n" + "[OmegahGenericMesh::mark_part_entities] Error! Cannot mark downward if the part dimension is 0.\n" " - part name: " << name << "\n" " - part dim : " << dim << "\n") Omega_h::Write downMarked; From 45366519479f97fbbcd8877363b12b613895a948 Mon Sep 17 00:00:00 2001 From: Cameron Smith Date: Thu, 30 Oct 2025 13:48:08 -0400 Subject: [PATCH 20/35] missed one name mismatch --- src/disc/omegah/Albany_OmegahGenericMesh.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/disc/omegah/Albany_OmegahGenericMesh.cpp b/src/disc/omegah/Albany_OmegahGenericMesh.cpp index 6bd79259c6..bb87818a32 100644 --- a/src/disc/omegah/Albany_OmegahGenericMesh.cpp +++ b/src/disc/omegah/Albany_OmegahGenericMesh.cpp @@ -197,7 +197,7 @@ mark_part_entities (const std::string& name, } } }; - Kokkos::parallel_for ("OmegahGenericMesh::set_part_entities::markDownward",m_mesh->nents(topo_dim(topo)),f); + Kokkos::parallel_for ("OmegahGenericMesh::mark_part_entities::markDownward",m_mesh->nents(topo_dim(topo)),f); m_mesh->add_tag(topo_dim(down_topo),name,1,read(downMarked)); upMarked = downMarked; From a17e616394dbc4cb031be8c85b52c9a2d26474d7 Mon Sep 17 00:00:00 2001 From: Cameron Smith Date: Thu, 30 Oct 2025 14:44:40 -0400 Subject: [PATCH 21/35] if adapt runs, but doesn't modify the mesh, the nodeset tag will exist --- src/disc/omegah/Albany_OmegahGenericMesh.cpp | 9 +++------ 1 file changed, 3 insertions(+), 6 deletions(-) diff --git a/src/disc/omegah/Albany_OmegahGenericMesh.cpp b/src/disc/omegah/Albany_OmegahGenericMesh.cpp index bb87818a32..89e2a11827 100644 --- a/src/disc/omegah/Albany_OmegahGenericMesh.cpp +++ b/src/disc/omegah/Albany_OmegahGenericMesh.cpp @@ -158,12 +158,9 @@ mark_part_entities (const std::string& name, " - part dim : " << dim << "\n" " - num ents : " << m_mesh->nents(dim) << "\n" " - array dim: " << is_entity_in_part.size() << "\n"); - TEUCHOS_TEST_FOR_EXCEPTION (m_mesh->has_tag(dim,name), std::runtime_error, - "[OmegahGenericMesh::mark_part_entities] Error! A tag with this name was already set.\n" - " - part name: " << name << "\n" - " - part dim : " << dim << "\n"); - - m_mesh->add_tag(dim,name,1,is_entity_in_part); + if( ! m_mesh->has_tag(dim,name) ) { + m_mesh->add_tag(dim,name,1,is_entity_in_part); + } if (markDownward) { TEUCHOS_TEST_FOR_EXCEPTION (dim==0, std::logic_error, From 03d96612a3fcd2aafbbf304be5ccf03654b68f4a Mon Sep 17 00:00:00 2001 From: Cameron Smith Date: Thu, 30 Oct 2025 15:32:59 -0400 Subject: [PATCH 22/35] modify 1d and 2d yaml settings and use --- .../omegah/Albany_OmegahDiscretization.cpp | 20 +++++++++++-------- .../AdvDiff/input_2d_adapt_omegah.yaml | 7 ++++--- 2 files changed, 16 insertions(+), 11 deletions(-) diff --git a/src/disc/omegah/Albany_OmegahDiscretization.cpp b/src/disc/omegah/Albany_OmegahDiscretization.cpp index 0daca621bf..aafd649db0 100644 --- a/src/disc/omegah/Albany_OmegahDiscretization.cpp +++ b/src/disc/omegah/Albany_OmegahDiscretization.cpp @@ -438,10 +438,6 @@ checkForAdaptation (const Teuchos::RCP& solution , return adapt_data; } - TEUCHOS_TEST_FOR_EXCEPTION (adapt_type!="Minimally-Oscillatory", std::runtime_error, //FIXME - need an SPR option - "Error! Adaptation type '" << adapt_type << "' not supported.\n" - " - valid choices: None, Minimally-Oscillatory\n"); - TEUCHOS_TEST_FOR_EXCEPTION (dxdp != Teuchos::null, std::runtime_error, "Error! the dxdp Thyra_MultiVector is expected to be null\n"); @@ -454,6 +450,9 @@ checkForAdaptation (const Teuchos::RCP& solution , } if (mesh->dim() == 1) { + TEUCHOS_TEST_FOR_EXCEPTION (adapt_type!="Minimally-Oscillatory", std::runtime_error, + "Error! Adaptation type '" << adapt_type << "' not supported.\n" + " - valid choices for 1D: None, Minimally-Oscillatory\n"); double tol = adapt_params.get("Max Hessian"); auto data = getLocalData(solution); // Simple check: refine if a proxy of the hessian of x is larger than a tolerance @@ -484,6 +483,9 @@ checkForAdaptation (const Teuchos::RCP& solution , << "... we will not adapt.\n"; return adapt_data; } + TEUCHOS_TEST_FOR_EXCEPTION (adapt_type!="SPR", std::runtime_error, + "Error! Adaptation type '" << adapt_type << "' not supported.\n" + " - valid choices for 2D: None, SPR\n"); #ifdef ALBANY_MESHFIELDS auto effectiveStrain = getEffectiveStrainRate(*mesh); @@ -504,7 +506,7 @@ checkForAdaptation (const Teuchos::RCP& solution , MeshField::Omegah::getTriangleElement(*mesh); MeshField::FieldElement coordFe(mesh->nelems(), coordField, shp, map); - const auto adaptRatio = 0.1; + const auto adaptRatio = adapt_params.get("Adapt Ratio",0.1); auto estimation = MeshField::SPR::Estimation(*mesh, effectiveStrain, recoveredStrainField, adaptRatio); @@ -578,7 +580,7 @@ adapt (const Teuchos::RCP& adaptData) std::cout << "mesh now has " << nelems << " total elements\n"; } } else if (ohMesh->dim() == 2 && isMeshfieldsEnabled()) { - // adapt + Omega_h::AdaptOpts opts(&(*ohMesh)); opts.xfer_opts.type_map[solution_dof_name()] = OMEGA_H_LINEAR_INTERP; opts.xfer_opts.type_map[std::string(solution_dof_name())+"_dot"] = OMEGA_H_LINEAR_INTERP; @@ -586,8 +588,10 @@ adapt (const Teuchos::RCP& adaptData) auto verbose = false; const auto tgtLength_oh = ohMesh->get_array(Omega_h::VERT, "tgtLength"); const auto isos = Omega_h::isos_from_lengths(tgtLength_oh); - const auto min_size = 0.1; - const auto max_size = 1.0; + + auto& adapt_params = m_disc_params->sublist("Mesh Adaptivity"); + const auto min_size = adapt_params.get("Minimum Size",0.08); + const auto max_size = adapt_params.get("Maximum Size",1.0); auto metric = Omega_h::clamp_metrics(ohMesh->nverts(), isos, min_size, max_size); Omega_h::grade_fix_adapt(&(*ohMesh), opts, metric, verbose); diff --git a/tests/corePDEs/AdvDiff/input_2d_adapt_omegah.yaml b/tests/corePDEs/AdvDiff/input_2d_adapt_omegah.yaml index 1f74bc6403..877f0a0267 100644 --- a/tests/corePDEs/AdvDiff/input_2d_adapt_omegah.yaml +++ b/tests/corePDEs/AdvDiff/input_2d_adapt_omegah.yaml @@ -34,9 +34,10 @@ ANONYMOUS: Method: Omegah Mesh Creation Method: Box2D Mesh Adaptivity: - Type: Minimally-Oscillatory - Max Hessian: 10.0 - Refining Factor: 2 + Type: SPR + Adapt Ratio: 0.1 + Minimum Size: 0.81 + Maximum Size: 1.0 Piro: Tempus: Integrator Name: Tempus Integrator From c42a72938b2077fe17c76e91d4da9fe85bb5df60 Mon Sep 17 00:00:00 2001 From: Cameron Smith Date: Thu, 30 Oct 2025 15:57:32 -0400 Subject: [PATCH 23/35] min size for ~300 tris --- tests/corePDEs/AdvDiff/input_2d_adapt_omegah.yaml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/corePDEs/AdvDiff/input_2d_adapt_omegah.yaml b/tests/corePDEs/AdvDiff/input_2d_adapt_omegah.yaml index 877f0a0267..71242a07d0 100644 --- a/tests/corePDEs/AdvDiff/input_2d_adapt_omegah.yaml +++ b/tests/corePDEs/AdvDiff/input_2d_adapt_omegah.yaml @@ -36,7 +36,7 @@ ANONYMOUS: Mesh Adaptivity: Type: SPR Adapt Ratio: 0.1 - Minimum Size: 0.81 + Minimum Size: 0.08 Maximum Size: 1.0 Piro: Tempus: From f4fc48ea71eb9ec3cfbdaabdfdb54f47c0ba157e Mon Sep 17 00:00:00 2001 From: Cameron Smith Date: Fri, 31 Oct 2025 14:58:37 -0400 Subject: [PATCH 24/35] better adapt parameter names --- src/disc/omegah/Albany_OmegahDiscretization.cpp | 4 ++-- tests/corePDEs/AdvDiff/input_2d_adapt_omegah.yaml | 6 +++--- 2 files changed, 5 insertions(+), 5 deletions(-) diff --git a/src/disc/omegah/Albany_OmegahDiscretization.cpp b/src/disc/omegah/Albany_OmegahDiscretization.cpp index aafd649db0..09d69f602b 100644 --- a/src/disc/omegah/Albany_OmegahDiscretization.cpp +++ b/src/disc/omegah/Albany_OmegahDiscretization.cpp @@ -590,8 +590,8 @@ adapt (const Teuchos::RCP& adaptData) const auto isos = Omega_h::isos_from_lengths(tgtLength_oh); auto& adapt_params = m_disc_params->sublist("Mesh Adaptivity"); - const auto min_size = adapt_params.get("Minimum Size",0.08); - const auto max_size = adapt_params.get("Maximum Size",1.0); + const auto min_size = adapt_params.get("Minimum Edge Length",0.08); + const auto max_size = adapt_params.get("Maximum Edge Length",1.0); auto metric = Omega_h::clamp_metrics(ohMesh->nverts(), isos, min_size, max_size); Omega_h::grade_fix_adapt(&(*ohMesh), opts, metric, verbose); diff --git a/tests/corePDEs/AdvDiff/input_2d_adapt_omegah.yaml b/tests/corePDEs/AdvDiff/input_2d_adapt_omegah.yaml index 71242a07d0..af8960cdb2 100644 --- a/tests/corePDEs/AdvDiff/input_2d_adapt_omegah.yaml +++ b/tests/corePDEs/AdvDiff/input_2d_adapt_omegah.yaml @@ -35,9 +35,9 @@ ANONYMOUS: Mesh Creation Method: Box2D Mesh Adaptivity: Type: SPR - Adapt Ratio: 0.1 - Minimum Size: 0.08 - Maximum Size: 1.0 + Adapt Ratio: 1.2 + Minimum Edge Length: 0.02 + Maximum Edge Length: 1.0 Piro: Tempus: Integrator Name: Tempus Integrator From 4ca39adfbaede042531e00229ec25e2b4cd88c31 Mon Sep 17 00:00:00 2001 From: Cameron Smith Date: Fri, 31 Oct 2025 15:56:21 -0400 Subject: [PATCH 25/35] trigger adaptation when error is large --- .../omegah/Albany_OmegahDiscretization.cpp | 35 ++++++++++--------- .../AdvDiff/input_2d_adapt_omegah.yaml | 1 + 2 files changed, 20 insertions(+), 16 deletions(-) diff --git a/src/disc/omegah/Albany_OmegahDiscretization.cpp b/src/disc/omegah/Albany_OmegahDiscretization.cpp index 09d69f602b..2a886ecf0b 100644 --- a/src/disc/omegah/Albany_OmegahDiscretization.cpp +++ b/src/disc/omegah/Albany_OmegahDiscretization.cpp @@ -510,24 +510,27 @@ checkForAdaptation (const Teuchos::RCP& solution , auto estimation = MeshField::SPR::Estimation(*mesh, effectiveStrain, recoveredStrainField, adaptRatio); - const auto tgtLength = MeshField::SPR::getSprSizeField(estimation, omf, coordFe); - Omega_h::Write tgtLength_oh(tgtLength); - mesh->add_tag(Omega_h::VERT, "tgtLength", 1, tgtLength_oh, false, - Omega_h::ArrayType::VectorND); - - { // write vtk - const auto outname = std::string("beforeAdapt2d"); - const std::string vtkFileName = outname + ".vtk"; - Omega_h::vtk::write_parallel(vtkFileName, &(*mesh), 2); - const std::string vtkFileName_edges = outname + "_edges.vtk"; - Omega_h::vtk::write_parallel(vtkFileName_edges, &(*mesh), 1); - } - - printTriCount(*mesh, "beforeAdapt"); - #endif - adapt_data->type = AdaptationType::Topology; + const auto [tgtLength, error] = MeshField::SPR::getSprSizeField(estimation, omf, coordFe); + const auto errorThreshold = adapt_params.get("Error Threshold",0.5); + std::cout << "Error: " << error << " Error Threshold: " << errorThreshold << '\n'; + if( error > errorThreshold ) { //trigger adaptation + Omega_h::Write tgtLength_oh(tgtLength); + mesh->add_tag(Omega_h::VERT, "tgtLength", 1, tgtLength_oh, false, + Omega_h::ArrayType::VectorND); + + { // write vtk + const auto outname = std::string("beforeAdapt2d"); + const std::string vtkFileName = outname + ".vtk"; + Omega_h::vtk::write_parallel(vtkFileName, &(*mesh), 2); + const std::string vtkFileName_edges = outname + "_edges.vtk"; + Omega_h::vtk::write_parallel(vtkFileName_edges, &(*mesh), 1); + } + printTriCount(*mesh, "beforeAdapt"); + adapt_data->type = AdaptationType::Topology; + } return adapt_data; + #endif //ALBANY_MESHFIELDS } else { //meshdim != 1 && meshdim != 2 std::cout << "Only 1D and 2D (with Meshfields enabled) Omega_h mesh " << "adaptation is supported ... we will not adapt.\n"; diff --git a/tests/corePDEs/AdvDiff/input_2d_adapt_omegah.yaml b/tests/corePDEs/AdvDiff/input_2d_adapt_omegah.yaml index af8960cdb2..d5a6012ab0 100644 --- a/tests/corePDEs/AdvDiff/input_2d_adapt_omegah.yaml +++ b/tests/corePDEs/AdvDiff/input_2d_adapt_omegah.yaml @@ -36,6 +36,7 @@ ANONYMOUS: Mesh Adaptivity: Type: SPR Adapt Ratio: 1.2 + Error Threshold: 0.3 Minimum Edge Length: 0.02 Maximum Edge Length: 1.0 Piro: From 238f9edb7733b3784b6773544f73cf1673765d0a Mon Sep 17 00:00:00 2001 From: Cameron Smith Date: Mon, 3 Nov 2025 11:39:34 -0500 Subject: [PATCH 26/35] spr: ctest passes --- tests/corePDEs/AdvDiff/input_2d_adapt_omegah.yaml | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/tests/corePDEs/AdvDiff/input_2d_adapt_omegah.yaml b/tests/corePDEs/AdvDiff/input_2d_adapt_omegah.yaml index d5a6012ab0..37fc2cec43 100644 --- a/tests/corePDEs/AdvDiff/input_2d_adapt_omegah.yaml +++ b/tests/corePDEs/AdvDiff/input_2d_adapt_omegah.yaml @@ -122,12 +122,12 @@ ANONYMOUS: Ifpack2 Settings: 'fact: ilut level-of-fill': 1.0e+00 Regression For Response 0: - Test Value: 6.825290995201e-02 + Test Value: 6.661786260006e-02 Relative Tolerance: 1.0e-04 Regression For Response 1: - Test Value: -8.259050370473e-04 + Test Value: -8.182824963026e-05 Relative Tolerance: 1.0e-04 Regression For Response 2: - Test Value: 2.170508705497e-03 + Test Value: 4.462238299510e-03 Relative Tolerance: 1.0e-04 ... From e7393b6b86835b3e65704206b74c333935b502ca Mon Sep 17 00:00:00 2001 From: Cameron Smith Date: Mon, 3 Nov 2025 12:26:36 -0500 Subject: [PATCH 27/35] spr: less scary if grepping for error --- src/disc/omegah/Albany_OmegahDiscretization.cpp | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/disc/omegah/Albany_OmegahDiscretization.cpp b/src/disc/omegah/Albany_OmegahDiscretization.cpp index 2a886ecf0b..59af18e515 100644 --- a/src/disc/omegah/Albany_OmegahDiscretization.cpp +++ b/src/disc/omegah/Albany_OmegahDiscretization.cpp @@ -512,7 +512,8 @@ checkForAdaptation (const Teuchos::RCP& solution , const auto [tgtLength, error] = MeshField::SPR::getSprSizeField(estimation, omf, coordFe); const auto errorThreshold = adapt_params.get("Error Threshold",0.5); - std::cout << "Error: " << error << " Error Threshold: " << errorThreshold << '\n'; + std::cout << "SPR Computed Error: " << error + << " Error Threshold: " << errorThreshold << '\n'; if( error > errorThreshold ) { //trigger adaptation Omega_h::Write tgtLength_oh(tgtLength); mesh->add_tag(Omega_h::VERT, "tgtLength", 1, tgtLength_oh, false, From 114a8e902299346c37fd4b0c03f41421fd7db17d Mon Sep 17 00:00:00 2001 From: Cameron Smith Date: Mon, 3 Nov 2025 12:26:52 -0500 Subject: [PATCH 28/35] spr: doc strings for yaml --- .../AdvDiff/input_2d_adapt_omegah.yaml | 21 ++++++++++++++++--- 1 file changed, 18 insertions(+), 3 deletions(-) diff --git a/tests/corePDEs/AdvDiff/input_2d_adapt_omegah.yaml b/tests/corePDEs/AdvDiff/input_2d_adapt_omegah.yaml index 37fc2cec43..420cc38269 100644 --- a/tests/corePDEs/AdvDiff/input_2d_adapt_omegah.yaml +++ b/tests/corePDEs/AdvDiff/input_2d_adapt_omegah.yaml @@ -34,11 +34,26 @@ ANONYMOUS: Method: Omegah Mesh Creation Method: Box2D Mesh Adaptivity: - Type: SPR + Type: SPR # superconvergent patch recovery + # 'Adapt Ratio' defines 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 results in larger mesh elements; + # this is a (nonlinear) scaling factor on the resulting + # size field. Adapt Ratio: 1.2 + # Error Thresold: + # SPR computes the sum of error norms raised to the power + # (2d/(2p+d)) over all elements, where d is the mesh + # dimension and p is the recovered field polynomial order. + # For this example d=2 and p=2. + # The 'Error Threshold' specifies the maxium allowed error + # before mesh adaptation is triggered; + # i.e., `if (error > threshold) runMeshAdaptation(...)` . Error Threshold: 0.3 - Minimum Edge Length: 0.02 - Maximum Edge Length: 1.0 + Minimum Edge Length: 0.02 # min desired edge length in physical space + Maximum Edge Length: 1.0 # max desired edge length in physical space Piro: Tempus: Integrator Name: Tempus Integrator From 6f532e22f3d22ac93ab373ba8d0d50203fafebec Mon Sep 17 00:00:00 2001 From: Cameron Smith Date: Mon, 3 Nov 2025 12:32:10 -0500 Subject: [PATCH 29/35] spr: meshfields requires omegah --- CMakeLists.txt | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/CMakeLists.txt b/CMakeLists.txt index aaffbb2d95..e6c8732d33 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -720,6 +720,11 @@ OPTION(ENABLE_MESHFIELDS "Flag to turn on the MeshFields field interface for computing error indicators (via super convergent patch recovery) to drive mesh adaptation" OFF) IF (ENABLE_MESHFIELDS) + IF (NOT ENABLE_OMEGAH) + MESSAGE(FATAL_ERROR + "The MeshFields field interface requires Omega_h to be enabled. \ + Reconfigure with -DENABLE_MESHFIELDS=ON and -DENABLE_OMEGAH=ON.") + ENDIF() MESSAGE("-- MeshFields field interface Enabled.") SET(ALBANY_MESHFIELDS TRUE) include (GetOrInstallMeshFields) From 57de97103b13e6dbe4be11ee15c6002b4938400d Mon Sep 17 00:00:00 2001 From: Cameron Smith Date: Mon, 3 Nov 2025 12:57:01 -0500 Subject: [PATCH 30/35] spr: requires meshfields --- src/disc/omegah/Albany_OmegahDiscretization.cpp | 2 ++ 1 file changed, 2 insertions(+) diff --git a/src/disc/omegah/Albany_OmegahDiscretization.cpp b/src/disc/omegah/Albany_OmegahDiscretization.cpp index 59af18e515..c043cff65e 100644 --- a/src/disc/omegah/Albany_OmegahDiscretization.cpp +++ b/src/disc/omegah/Albany_OmegahDiscretization.cpp @@ -39,6 +39,7 @@ namespace { return Omega_h::project_by_fit(&mesh, effectiveStrain); } + #ifdef ALBANY_MESHFIELDS template void setFieldAtVertices(Omega_h::Mesh &mesh, Omega_h::Reals recoveredStrain, ShapeField field) { @@ -48,6 +49,7 @@ namespace { MeshField::parallel_for(ExecutionSpace(), {0}, {mesh.nverts()}, setFieldAtVertices, "setFieldAtVertices"); } + #endif void printTriCount(Omega_h::Mesh &mesh, std::string_view prefix) { const auto nTri = mesh.nglobal_ents(2); From 8b5ca71c2eb0d819737b0dfb79edcecb467ef6ee Mon Sep 17 00:00:00 2001 From: Cameron Smith Date: Tue, 4 Nov 2025 12:13:50 -0500 Subject: [PATCH 31/35] controls to silence output --- .../omegah/Albany_OmegahDiscretization.cpp | 52 ++++++++++++------- .../AdvDiff/input_2d_adapt_omegah.yaml | 2 + 2 files changed, 34 insertions(+), 20 deletions(-) diff --git a/src/disc/omegah/Albany_OmegahDiscretization.cpp b/src/disc/omegah/Albany_OmegahDiscretization.cpp index c043cff65e..52e9634136 100644 --- a/src/disc/omegah/Albany_OmegahDiscretization.cpp +++ b/src/disc/omegah/Albany_OmegahDiscretization.cpp @@ -439,6 +439,7 @@ checkForAdaptation (const Teuchos::RCP& solution , if (adapt_type=="None") { return adapt_data; } + const auto verbose = adapt_params.get("Verbose",false); TEUCHOS_TEST_FOR_EXCEPTION (dxdp != Teuchos::null, std::runtime_error, "Error! the dxdp Thyra_MultiVector is expected to be null\n"); @@ -480,9 +481,11 @@ checkForAdaptation (const Teuchos::RCP& solution , return adapt_data; } else if (mesh->dim() == 2) { if (!isMeshfieldsEnabled()) { - std::cout << "Warning: 2D Omega_h mesh adaptation requires Meshfields. " - << "Configure Albany with ENABLE_MESHFIELDS=ON to enable it. " - << "... we will not adapt.\n"; + if (!mesh->comm()->rank()) { + std::cout << "Warning: 2D Omega_h mesh adaptation requires Meshfields. " + << "Configure Albany with ENABLE_MESHFIELDS=ON to enable it. " + << "... we will not adapt.\n"; + } return adapt_data; } TEUCHOS_TEST_FOR_EXCEPTION (adapt_type!="SPR", std::runtime_error, @@ -514,14 +517,19 @@ checkForAdaptation (const Teuchos::RCP& solution , const auto [tgtLength, error] = MeshField::SPR::getSprSizeField(estimation, omf, coordFe); const auto errorThreshold = adapt_params.get("Error Threshold",0.5); - std::cout << "SPR Computed Error: " << error - << " Error Threshold: " << errorThreshold << '\n'; + if(verbose) { + //FIXME - should this be a per-rank output? + // - does getSprSizeField have a reduction? + std::cout << "SPR Computed Error: " << error + << " Error Threshold: " << errorThreshold << '\n'; + } if( error > errorThreshold ) { //trigger adaptation Omega_h::Write tgtLength_oh(tgtLength); mesh->add_tag(Omega_h::VERT, "tgtLength", 1, tgtLength_oh, false, Omega_h::ArrayType::VectorND); - { // write vtk + const auto writeVtk = adapt_params.get("Write VTK Files",false); + if( writeVtk ) { const auto outname = std::string("beforeAdapt2d"); const std::string vtkFileName = outname + ".vtk"; Omega_h::vtk::write_parallel(vtkFileName, &(*mesh), 2); @@ -529,14 +537,16 @@ checkForAdaptation (const Teuchos::RCP& solution , Omega_h::vtk::write_parallel(vtkFileName_edges, &(*mesh), 1); } - printTriCount(*mesh, "beforeAdapt"); + if(verbose) printTriCount(*mesh, "beforeAdapt"); adapt_data->type = AdaptationType::Topology; } return adapt_data; #endif //ALBANY_MESHFIELDS } else { //meshdim != 1 && meshdim != 2 - std::cout << "Only 1D and 2D (with Meshfields enabled) Omega_h mesh " - << "adaptation is supported ... we will not adapt.\n"; + if (!mesh->comm()->rank()) { + std::cout << "Only 1D and 2D (with Meshfields enabled) Omega_h mesh " + << "adaptation is supported ... we will not adapt.\n"; + } return adapt_data; } } @@ -556,8 +566,8 @@ adapt (const Teuchos::RCP& adaptData) TEUCHOS_TEST_FOR_EXCEPTION (ohMesh->dim()!=1 && ohMesh->dim()!=2, std::runtime_error, "Error! Adaptation not supported for this mesh. We only implemented simple 1d and 2d cases.\n"); - std::string beforeAdaptName = "before_adapt" + std::to_string(adaptCount) + ".vtk"; - Omega_h::vtk::write_parallel(beforeAdaptName, ohMesh.get()); + auto& adapt_params = m_disc_params->sublist("Mesh Adaptivity"); + const auto verbose = adapt_params.get("Verbose",false); if (ohMesh->dim() == 1) { // Note: the code below is hard-coding a simple adaptation for a 1d mesh, @@ -566,11 +576,12 @@ adapt (const Teuchos::RCP& adaptData) const auto desired_nelems = nelems*2; Omega_h::AdaptOpts opts(&(*ohMesh)); + opts.verbosity = (verbose ? Omega_h::EACH_ADAPT : Omega_h::SILENT); opts.xfer_opts.type_map[solution_dof_name()] = OMEGA_H_LINEAR_INTERP; opts.xfer_opts.type_map[std::string(solution_dof_name())+"_dot"] = OMEGA_H_LINEAR_INTERP; while (double(nelems) < desired_nelems) { if (!ohMesh->has_tag(0, "metric")) { - std::cout << "mesh had no metric, adding implied and adapting to it\n"; + if(verbose) std::cout << "mesh had no metric, adding implied and adapting to it\n"; Omega_h::add_implied_metric_tag(ohMesh.get()); Omega_h::adapt(ohMesh.get(), opts); nelems = ohMesh->nglobal_ents(ohMesh->dim()); @@ -580,32 +591,33 @@ adapt (const Teuchos::RCP& adaptData) auto const metric_ncomps = Omega_h::divide_no_remainder(metrics.size(), ohMesh->nverts()); ohMesh->add_tag(0, "metric", metric_ncomps, metrics); - std::cout << "adapting to scaled metric\n"; + if(verbose) std::cout << "adapting to scaled metric\n"; Omega_h::adapt(ohMesh.get(), opts); nelems = ohMesh->nglobal_ents(ohMesh->dim()); - std::cout << "mesh now has " << nelems << " total elements\n"; + if(verbose) std::cout << "mesh now has " << nelems << " total elements\n"; } } else if (ohMesh->dim() == 2 && isMeshfieldsEnabled()) { Omega_h::AdaptOpts opts(&(*ohMesh)); + opts.verbosity = (verbose ? Omega_h::EACH_ADAPT : Omega_h::SILENT); opts.xfer_opts.type_map[solution_dof_name()] = OMEGA_H_LINEAR_INTERP; opts.xfer_opts.type_map[std::string(solution_dof_name())+"_dot"] = OMEGA_H_LINEAR_INTERP; - auto verbose = false; const auto tgtLength_oh = ohMesh->get_array(Omega_h::VERT, "tgtLength"); const auto isos = Omega_h::isos_from_lengths(tgtLength_oh); - - auto& adapt_params = m_disc_params->sublist("Mesh Adaptivity"); const auto min_size = adapt_params.get("Minimum Edge Length",0.08); const auto max_size = adapt_params.get("Maximum Edge Length",1.0); auto metric = Omega_h::clamp_metrics(ohMesh->nverts(), isos, min_size, max_size); Omega_h::grade_fix_adapt(&(*ohMesh), opts, metric, verbose); - printTriCount(*ohMesh, "afterAdapt"); + if(verbose) printTriCount(*ohMesh, "afterAdapt"); } - std::string afterAdaptName = "after_adapt" + std::to_string(adaptCount) + ".vtk"; - Omega_h::vtk::write_parallel(afterAdaptName, ohMesh.get()); + const auto writeVtk = adapt_params.get("Write VTK Files",false); + if( writeVtk ) { + std::string afterAdaptName = "after_adapt" + std::to_string(adaptCount) + ".vtk"; + Omega_h::vtk::write_parallel(afterAdaptName, ohMesh.get()); + } //create node and side set tags m_mesh_struct->createNodeSets(); diff --git a/tests/corePDEs/AdvDiff/input_2d_adapt_omegah.yaml b/tests/corePDEs/AdvDiff/input_2d_adapt_omegah.yaml index 420cc38269..afc22f9062 100644 --- a/tests/corePDEs/AdvDiff/input_2d_adapt_omegah.yaml +++ b/tests/corePDEs/AdvDiff/input_2d_adapt_omegah.yaml @@ -54,6 +54,8 @@ ANONYMOUS: Error Threshold: 0.3 Minimum Edge Length: 0.02 # min desired edge length in physical space Maximum Edge Length: 1.0 # max desired edge length in physical space + Write VTK Files: false # write VTK files before and after adaptation + Verbose: false # write additional mesh adaptation info to stdout Piro: Tempus: Integrator Name: Tempus Integrator From 29779ad01aa04dd9fc82194d6c28e51be751605b Mon Sep 17 00:00:00 2001 From: Cameron Smith Date: Wed, 5 Nov 2025 20:23:26 -0800 Subject: [PATCH 32/35] need uvm to run on perlmutter gpu --- cmake/GetOrInstallOmegah.cmake | 3 +++ 1 file changed, 3 insertions(+) diff --git a/cmake/GetOrInstallOmegah.cmake b/cmake/GetOrInstallOmegah.cmake index 2e261238c0..6b8c312f60 100644 --- a/cmake/GetOrInstallOmegah.cmake +++ b/cmake/GetOrInstallOmegah.cmake @@ -57,6 +57,9 @@ else () option (Omega_h_USE_Kokkos "Use Kokkos as a backend" ON) option (Omega_h_USE_MPI "Use MPI for parallelism" ON) set (MPIEXEC_EXECUTABLE ${Albany_CXX_COMPILER}) + if (Kokkos_ENABLE_CUDA_UVM) + option (Omega_h_MEM_SPACE_SHARED "enabled shared memory space" ON) + endif() message (STATUS " *** Begin of Omega_h configuration ***") FetchContent_MakeAvailable (Omega_h) From 7612a3fa56dfb013396368ee9b220b7cd217faad Mon Sep 17 00:00:00 2001 From: Cameron Smith Date: Thu, 6 Nov 2025 10:11:06 -0800 Subject: [PATCH 33/35] spr: append adapt count to pre-adapt vtk filenames --- .../omegah/Albany_OmegahDiscretization.cpp | 20 ++++++++++--------- 1 file changed, 11 insertions(+), 9 deletions(-) diff --git a/src/disc/omegah/Albany_OmegahDiscretization.cpp b/src/disc/omegah/Albany_OmegahDiscretization.cpp index 52e9634136..31f449dc19 100644 --- a/src/disc/omegah/Albany_OmegahDiscretization.cpp +++ b/src/disc/omegah/Albany_OmegahDiscretization.cpp @@ -528,14 +528,6 @@ checkForAdaptation (const Teuchos::RCP& solution , mesh->add_tag(Omega_h::VERT, "tgtLength", 1, tgtLength_oh, false, Omega_h::ArrayType::VectorND); - const auto writeVtk = adapt_params.get("Write VTK Files",false); - if( writeVtk ) { - const auto outname = std::string("beforeAdapt2d"); - const std::string vtkFileName = outname + ".vtk"; - Omega_h::vtk::write_parallel(vtkFileName, &(*mesh), 2); - const std::string vtkFileName_edges = outname + "_edges.vtk"; - Omega_h::vtk::write_parallel(vtkFileName_edges, &(*mesh), 1); - } if(verbose) printTriCount(*mesh, "beforeAdapt"); adapt_data->type = AdaptationType::Topology; @@ -568,6 +560,17 @@ adapt (const Teuchos::RCP& adaptData) auto& adapt_params = m_disc_params->sublist("Mesh Adaptivity"); const auto verbose = adapt_params.get("Verbose",false); + const auto writeVtk = adapt_params.get("Write VTK Files",false); + + if( writeVtk ) { + const auto outname = std::string("before_adapt") + std::to_string(adaptCount); + const std::string vtkFileName = outname + ".vtk"; + Omega_h::vtk::write_parallel(vtkFileName, &(*ohMesh), ohMesh->dim()); + if (ohMesh->dim() == 2) { + const std::string vtkFileName_edges = outname + "_edges.vtk"; + Omega_h::vtk::write_parallel(vtkFileName_edges, &(*ohMesh), Omega_h::EDGE); + } + } if (ohMesh->dim() == 1) { // Note: the code below is hard-coding a simple adaptation for a 1d mesh, @@ -613,7 +616,6 @@ adapt (const Teuchos::RCP& adaptData) if(verbose) printTriCount(*ohMesh, "afterAdapt"); } - const auto writeVtk = adapt_params.get("Write VTK Files",false); if( writeVtk ) { std::string afterAdaptName = "after_adapt" + std::to_string(adaptCount) + ".vtk"; Omega_h::vtk::write_parallel(afterAdaptName, ohMesh.get()); From 15ef0cc7124b5e38d10d644094e02e7d38214576 Mon Sep 17 00:00:00 2001 From: Cameron Smith Date: Thu, 6 Nov 2025 10:51:24 -0800 Subject: [PATCH 34/35] spr: loosen response tolerances for cuda uvm runs final mesh looks good --- tests/corePDEs/AdvDiff/input_2d_adapt_omegah.yaml | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/tests/corePDEs/AdvDiff/input_2d_adapt_omegah.yaml b/tests/corePDEs/AdvDiff/input_2d_adapt_omegah.yaml index afc22f9062..e7f599aae2 100644 --- a/tests/corePDEs/AdvDiff/input_2d_adapt_omegah.yaml +++ b/tests/corePDEs/AdvDiff/input_2d_adapt_omegah.yaml @@ -140,11 +140,11 @@ ANONYMOUS: 'fact: ilut level-of-fill': 1.0e+00 Regression For Response 0: Test Value: 6.661786260006e-02 - Relative Tolerance: 1.0e-04 + Absolute Tolerance: 2.0e-05 Regression For Response 1: Test Value: -8.182824963026e-05 - Relative Tolerance: 1.0e-04 + Absolute Tolerance: 3.0e-05 Regression For Response 2: Test Value: 4.462238299510e-03 - Relative Tolerance: 1.0e-04 + Absolute Tolerance: 2.0e-04 ... From 36780166f7302cf78e005de611d3d08cfe87c400 Mon Sep 17 00:00:00 2001 From: Cameron Smith Date: Fri, 7 Nov 2025 10:42:15 -0500 Subject: [PATCH 35/35] already linked to albanyLib thanks @bartgol --- src/landIce/CMakeLists.txt | 6 ------ 1 file changed, 6 deletions(-) diff --git a/src/landIce/CMakeLists.txt b/src/landIce/CMakeLists.txt index b799f6b5de..c84771a98b 100644 --- a/src/landIce/CMakeLists.txt +++ b/src/landIce/CMakeLists.txt @@ -218,12 +218,6 @@ endif() add_library(landIce ${Albany_LIBRARY_TYPE} ${SOURCES}) set_target_properties(landIce PROPERTIES PUBLIC_HEADER "${HEADERS}") target_link_libraries(landIce PUBLIC albanyLib) -if (ALBANY_OMEGAH) - target_link_libraries(landIce PUBLIC Omega_h::omega_h) -endif() -if (ALBANY_MESHFIELDS) - target_link_libraries(landIce PUBLIC meshfields) -endif() if (ALBANY_SUPPRESS_TRILINOS_WARNINGS) target_include_directories(landIce SYSTEM PUBLIC "${Trilinos_INCLUDE_DIRS};${Trilinos_TPL_INCLUDE_DIRS}")