diff --git a/CMakeLists.txt b/CMakeLists.txt index 55e084356c..e6c8732d33 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -715,6 +715,23 @@ 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) + 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) +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..5172459e8b --- /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/main + 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() 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) 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 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} diff --git a/src/disc/omegah/Albany_OmegahDiscretization.cpp b/src/disc/omegah/Albany_OmegahDiscretization.cpp index 6917ec9b1f..66ed26f15d 100644 --- a/src/disc/omegah/Albany_OmegahDiscretization.cpp +++ b/src/disc/omegah/Albany_OmegahDiscretization.cpp @@ -2,14 +2,62 @@ #include "Albany_OmegahUtils.hpp" #include "Albany_StringUtils.hpp" #include "Albany_ThyraUtils.hpp" +#include "Albany_KokkosTypes.hpp" // PHX::Device:: #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 + +#ifdef ALBANY_MESHFIELDS +#include +#include +#include +#endif #include +using ExecutionSpace = PHX::Device::execution_space; +using MemorySpace = PHX::Device::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); + } + + #ifdef ALBANY_MESHFIELDS + 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"); + } + #endif + + 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) { @@ -385,58 +433,114 @@ 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, - "Error! Adaptation type '" << adapt_type << "' not supported.\n" - " - valid choices: None, Minimally-Oscillatory\n"); + 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"); 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) { + 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 + // 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()) { + 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, + "Error! Adaptation type '" << adapt_type << "' not supported.\n" + " - valid choices for 2D: None, SPR\n"); + + #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 = adapt_params.get("Adapt Ratio",0.1); + auto estimation = + MeshField::SPR::Estimation(*mesh, effectiveStrain, recoveredStrainField, adaptRatio); + + const auto [tgtLength, error] = MeshField::SPR::getSprSizeField(estimation, omf, coordFe); + const auto errorThreshold = adapt_params.get("Error Threshold",0.5); + 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); - return adapt_data; + if(verbose) printTriCount(*mesh, "beforeAdapt"); + adapt_data->type = AdaptationType::Topology; + } + return adapt_data; + #endif //ALBANY_MESHFIELDS + } else { //meshdim != 1 && meshdim != 2 + 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; + } } void OmegahDiscretization:: @@ -451,40 +555,71 @@ 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"); - - 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()); + 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"); + + 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, + // 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.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")) { + 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()); + } + 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); + if(verbose) std::cout << "adapting to scaled metric\n"; Omega_h::adapt(ohMesh.get(), opts); nelems = ohMesh->nglobal_ents(ohMesh->dim()); + if(verbose) 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 && 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; + + 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 = 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); + + if(verbose) printTriCount(*ohMesh, "afterAdapt"); } - std::string afterAdaptName = "after_adapt" + std::to_string(adaptCount) + ".vtk"; - Omega_h::vtk::write_parallel(afterAdaptName, ohMesh.get()); + 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/src/disc/omegah/Albany_OmegahGenericMesh.cpp b/src/disc/omegah/Albany_OmegahGenericMesh.cpp index e27c2fb247..89e2a11827 100644 --- a/src/disc/omegah/Albany_OmegahGenericMesh.cpp +++ b/src/disc/omegah/Albany_OmegahGenericMesh.cpp @@ -147,27 +147,24 @@ 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" - " - 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, - "[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; @@ -197,7 +194,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; diff --git a/src/landIce/CMakeLists.txt b/src/landIce/CMakeLists.txt index dccedb0b47..c84771a98b 100644 --- a/src/landIce/CMakeLists.txt +++ b/src/landIce/CMakeLists.txt @@ -218,9 +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_SUPPRESS_TRILINOS_WARNINGS) target_include_directories(landIce SYSTEM PUBLIC "${Trilinos_INCLUDE_DIRS};${Trilinos_TPL_INCLUDE_DIRS}") diff --git a/tests/corePDEs/AdvDiff/CMakeLists.txt b/tests/corePDEs/AdvDiff/CMakeLists.txt index 064d55662f..39ce1f716c 100644 --- a/tests/corePDEs/AdvDiff/CMakeLists.txt +++ b/tests/corePDEs/AdvDiff/CMakeLists.txt @@ -34,15 +34,23 @@ 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) - set_tests_properties(${testName}_1d_adapt_Omega_h PROPERTIES LABELS "Demo;Forward;Serial;Adapt") + 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 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") + 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() 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..e7f599aae2 --- /dev/null +++ b/tests/corePDEs/AdvDiff/input_2d_adapt_omegah.yaml @@ -0,0 +1,150 @@ +%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: + Number of Elements: [10,10] + Method: Omegah + Mesh Creation Method: Box2D + Mesh Adaptivity: + 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 # 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 + 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.661786260006e-02 + Absolute Tolerance: 2.0e-05 + Regression For Response 1: + Test Value: -8.182824963026e-05 + Absolute Tolerance: 3.0e-05 + Regression For Response 2: + Test Value: 4.462238299510e-03 + Absolute Tolerance: 2.0e-04 +...