From a75d2c3248ec2e542f52410e4db50f54dab4e443 Mon Sep 17 00:00:00 2001 From: "Henry R. Swantner" Date: Wed, 14 Nov 2018 15:17:20 -0700 Subject: [PATCH 1/5] New for GCC 7.3.0 --- ...lRequestLinuxGCC7.3.0TestingSettings.cmake | 22 ++++++++++++ .../std/sems/PullRequestGCC7.3.0TestingEnv.sh | 36 +++++++++++++++++++ 2 files changed, 58 insertions(+) create mode 100644 cmake/std/PullRequestLinuxGCC7.3.0TestingSettings.cmake create mode 100644 cmake/std/sems/PullRequestGCC7.3.0TestingEnv.sh diff --git a/cmake/std/PullRequestLinuxGCC7.3.0TestingSettings.cmake b/cmake/std/PullRequestLinuxGCC7.3.0TestingSettings.cmake new file mode 100644 index 000000000000..3879e4df5f54 --- /dev/null +++ b/cmake/std/PullRequestLinuxGCC7.3.0TestingSettings.cmake @@ -0,0 +1,22 @@ +# This file contains the options needed to both run the pull request testing +# for Trilinos for the Linux GCC 7.3.0 pull request testing builds, and to reproduce +# the errors reported by those builds. Prior to using this this file, the +# appropriate set of SEMS modules must be loaded and accessible through the +# SEMS NFS mount. (See the sems/PullRequestGCC*TestingEnv.sh files.) + +# Usage: cmake -C PullRequestLinuxGCC7.3.0TestingSettings.cmake + +# Misc options typically added by CI testing mode in TriBITS + +# Use the below option only when submitting to the dashboard +#set (CTEST_USE_LAUNCHERS ON CACHE BOOL "Set by default for PR testing") + +set (MPI_EXEC_PRE_NUMPROCS_FLAGS "--bind-to;none" CACHE STRING "Set by default for PR testing") +# NOTE: The above is a workaround for the problem of having threads on MPI +# ranks bind to the same cores (see #2422). + +# Disable just one Teko sub-unit test that fails with openmpi 1.10 (#2712) +set (Teko_DISABLE_LSCSTABALIZED_TPETRA_ALPAH_INV_D ON CACHE BOOL "Temporarily disabled in PR testing") + +include("${CMAKE_CURRENT_LIST_DIR}/PullRequestLinuxCommonTestingSettings.cmake") + diff --git a/cmake/std/sems/PullRequestGCC7.3.0TestingEnv.sh b/cmake/std/sems/PullRequestGCC7.3.0TestingEnv.sh new file mode 100644 index 000000000000..54844c260a3e --- /dev/null +++ b/cmake/std/sems/PullRequestGCC7.3.0TestingEnv.sh @@ -0,0 +1,36 @@ +# This script can be used to load the appropriate environment for the +# GCC 7.3 Pull Request testing build on a Linux machine that has access to +# the SEMS NFS mount. + +# usage: $ source PullRequestGCC7.3TestingEnv.sh + +source /projects/sems/modulefiles/utils/sems-modules-init.sh + +module load sems-gcc/7.3.0 +module load sems-openmpi/1.10.1 +module load sems-python/2.7.9 +module load sems-git/2.10.1 +module load sems-boost/1.63.0/base +module load sems-zlib/1.2.8/base +module load sems-hdf5/1.8.12/parallel +module load sems-netcdf/4.4.1/exo_parallel +module load sems-parmetis/4.0.3/parallel +module load sems-scotch/6.0.3/nopthread_64bit_parallel +module load sems-superlu/4.3/base + +# Load the SEMS CMake Module +# - One of the SEMS modules will load CMake 3.4.x also, +# so this will pull in the SEMS cmake 3.10.3 version +# for Trilinos compatibility. +module load sems-cmake/3.10.3 + +# Using CMake and Ninja modules from the ATDM project space. +# SEMS does not yet supply a recent enough version of CMake +# for the single configure/build/test capability. We are also +# using a custom version of Ninja (with Fortran support not +# available in main-line Ninja) to significantly speed up +# compile and link times. +module load atdm-env +module load atdm-cmake/3.11.1 +module load atdm-ninja_fortran/1.7.2 + From 9bb9e68c1927d83b7e8d83fe7e9cd45df73a479c Mon Sep 17 00:00:00 2001 From: "Henry R. Swantner" Date: Wed, 14 Nov 2018 16:04:46 -0700 Subject: [PATCH 2/5] Issue 3832: Added lines for GCC 7.3.0 --- cmake/std/PullRequestLinuxDriver-Test.sh | 9 +++++++++ cmake/std/sems/PullRequestGCC7.3.0TestingEnv.sh | 4 ++-- 2 files changed, 11 insertions(+), 2 deletions(-) diff --git a/cmake/std/PullRequestLinuxDriver-Test.sh b/cmake/std/PullRequestLinuxDriver-Test.sh index 68ce7c6a3f99..04fe6673a9a8 100755 --- a/cmake/std/PullRequestLinuxDriver-Test.sh +++ b/cmake/std/PullRequestLinuxDriver-Test.sh @@ -147,6 +147,13 @@ elif [ "Trilinos_pullrequest_gcc_4.9.3_SERIAL" == "${JOB_BASE_NAME:?}" ] ; then echo -e "There was an issue loading the gcc environment. The error code was: $ierror" exit $ierror fi +elif [ "Trilinos_pullrequest_gcc_7.3.0" == "${JOB_BASE_NAME:?}" ] ; then + source ${TRILINOS_DRIVER_SRC_DIR}/cmake/std/sems/PullRequestGCC7.3.0TestingEnv.sh + ierror=$? + if [[ $ierror != 0 ]]; then + echo -e "There was an issue loading the gcc environment. The error code was: $ierror" + exit $ierror + fi elif [ "Trilinos_pullrequest_intel_17.0.1" == "${JOB_BASE_NAME:?}" ] ; then source ${TRILINOS_DRIVER_SRC_DIR}/cmake/std/sems/PullRequestIntel17.0.1TestingEnv.sh ierror=$? @@ -230,6 +237,8 @@ else elif [ "Trilinos_pullrequest_gcc_4.9.3_SERIAL" == "${JOB_BASE_NAME:?}" ]; then # TODO: Update this to use a 4.9.3 SERIAL testing environment script. CONFIG_SCRIPT=PullRequestLinuxGCC4.9.3TestingSettingsSERIAL.cmake + elif [ "Trilinos_pullrequest_gcc_7.3.0" == "${JOB_BASE_NAME:?}" ]; then + CONFIG_SCRIPT=PullRequestLinuxGCC7.3.0TestingSettings.cmake fi fi diff --git a/cmake/std/sems/PullRequestGCC7.3.0TestingEnv.sh b/cmake/std/sems/PullRequestGCC7.3.0TestingEnv.sh index 54844c260a3e..fee3ed73f8a6 100644 --- a/cmake/std/sems/PullRequestGCC7.3.0TestingEnv.sh +++ b/cmake/std/sems/PullRequestGCC7.3.0TestingEnv.sh @@ -1,8 +1,8 @@ # This script can be used to load the appropriate environment for the -# GCC 7.3 Pull Request testing build on a Linux machine that has access to +# GCC 7.3.0 Pull Request testing build on a Linux machine that has access to # the SEMS NFS mount. -# usage: $ source PullRequestGCC7.3TestingEnv.sh +# usage: $ source PullRequestGCC7.3.0TestingEnv.sh source /projects/sems/modulefiles/utils/sems-modules-init.sh From 0376131b0c9e3e930918cbcaa8b4b521387684f8 Mon Sep 17 00:00:00 2001 From: Chris Siefert Date: Wed, 21 Nov 2018 13:57:42 -0700 Subject: [PATCH 3/5] Ctest: Gemina build fixes? --- .../geminga/TrilinosCTestDriverCore.geminga.gcc-cuda.cmake | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/cmake/ctest/drivers/geminga/TrilinosCTestDriverCore.geminga.gcc-cuda.cmake b/cmake/ctest/drivers/geminga/TrilinosCTestDriverCore.geminga.gcc-cuda.cmake index ed268c67ee80..2f40e331661a 100644 --- a/cmake/ctest/drivers/geminga/TrilinosCTestDriverCore.geminga.gcc-cuda.cmake +++ b/cmake/ctest/drivers/geminga/TrilinosCTestDriverCore.geminga.gcc-cuda.cmake @@ -102,6 +102,12 @@ MACRO(TRILINOS_SYSTEM_SPECIFIC_CTEST_DRIVER) "-DTPL_ENABLE_CUSPARSE:BOOL=ON" "-DTPL_ENABLE_HWLOC:BOOL=OFF" + # Host Blas is required (https://github.com/kokkos/kokkos-kernels/issues/347) for Kokkos-Kernels to build correctly + "-DTPL_ENABLE_BLAS:BOOL=ON" + "-DTPL_ENABLE_LAPACK:BOOL=ON" + "-DTPL_BLAS_LIBRARIES=/usr/lib64/libblas.so.3.2.1" + "-DTPL_LAPACK_LIBRARIES=/usr/lib64/liblapack.so.3.2.1" + ### PACKAGE CONFIGURATION ### "-DKokkos_ENABLE_Cuda:BOOL=ON" "-DKokkos_ENABLE_Cuda_UVM:BOOL=ON" From 5a53209942eff59dcf61e0356e56f06ee3a20227 Mon Sep 17 00:00:00 2001 From: Chris Siefert Date: Wed, 21 Nov 2018 14:01:43 -0700 Subject: [PATCH 4/5] Ctest: Try enabling Fortran on rocketman. --- .../rocketman/TrilinosCTestDriverCore.rocketman.gcc.cmake | 2 -- 1 file changed, 2 deletions(-) diff --git a/cmake/ctest/drivers/rocketman/TrilinosCTestDriverCore.rocketman.gcc.cmake b/cmake/ctest/drivers/rocketman/TrilinosCTestDriverCore.rocketman.gcc.cmake index fcb18faa181c..e2e187a900c9 100644 --- a/cmake/ctest/drivers/rocketman/TrilinosCTestDriverCore.rocketman.gcc.cmake +++ b/cmake/ctest/drivers/rocketman/TrilinosCTestDriverCore.rocketman.gcc.cmake @@ -88,8 +88,6 @@ MACRO(TRILINOS_SYSTEM_SPECIFIC_CTEST_DRIVER) "-DTrilinos_ENABLE_COMPLEX:BOOL=OFF" - "-DTrilinos_ENABLE_Fortran=OFF" - "-DHDF5_INCLUDE_DIRS:FILEPATH=$ENV{SEMS_HDF5_ROOT}/include" "-DHDF5_LIBRARY_DIRS:FILEPATH=$ENV{SEMS_HDF5_ROOT}/lib" From bf5b85c524a0fc5977e92564443bbffb4ef220a2 Mon Sep 17 00:00:00 2001 From: "Edward G. Phillips" Date: Tue, 20 Nov 2018 13:19:14 -0700 Subject: [PATCH 5/5] Panzer: Adapted CurlLaplacianExample into a mixed version that uses both HCurl and HDiv elements in 3D or HCurl and HVol elements in 2D. Exposed HVol basis type and defined its basis values. Modified L2Projection to allow HVol elements. Removed Mathematica file --- .../adapters-stk/example/CMakeLists.txt | 1 + .../MixedCurlLaplacianExample/CMakeLists.txt | 229 ++++++ .../Example_BCStrategy_Dirichlet_Constant.hpp | 80 ++ ...ple_BCStrategy_Dirichlet_Constant_impl.hpp | 144 ++++ .../Example_BCStrategy_Factory.hpp | 87 ++ .../Example_ClosureModel_Factory.hpp | 71 ++ ...e_ClosureModel_Factory_TemplateBuilder.hpp | 66 ++ .../Example_ClosureModel_Factory_impl.hpp | 322 ++++++++ .../Example_CurlLaplacianEquationSet.hpp | 94 +++ .../Example_CurlLaplacianEquationSet_impl.hpp | 321 ++++++++ .../Example_CurlSolution.hpp | 94 +++ .../Example_CurlSolution_impl.hpp | 108 +++ .../Example_EquationSetFactory.hpp | 93 +++ .../Example_SimpleSource.hpp | 94 +++ .../Example_SimpleSource_impl.hpp | 104 +++ .../convergence_rate.py | 87 ++ .../MixedCurlLaplacianExample/main.cpp | 768 ++++++++++++++++++ .../disc-fe/src/Panzer_BasisValues2.hpp | 5 + .../disc-fe/src/Panzer_BasisValues2_impl.hpp | 105 ++- .../src/Panzer_IntrepidBasisFactory.hpp | 11 + .../disc-fe/src/Panzer_L2Projection_impl.hpp | 4 +- .../panzer/disc-fe/src/Panzer_PureBasis.cpp | 5 + .../panzer/disc-fe/src/Panzer_PureBasis.hpp | 4 +- 23 files changed, 2891 insertions(+), 6 deletions(-) create mode 100644 packages/panzer/adapters-stk/example/MixedCurlLaplacianExample/CMakeLists.txt create mode 100644 packages/panzer/adapters-stk/example/MixedCurlLaplacianExample/Example_BCStrategy_Dirichlet_Constant.hpp create mode 100644 packages/panzer/adapters-stk/example/MixedCurlLaplacianExample/Example_BCStrategy_Dirichlet_Constant_impl.hpp create mode 100644 packages/panzer/adapters-stk/example/MixedCurlLaplacianExample/Example_BCStrategy_Factory.hpp create mode 100644 packages/panzer/adapters-stk/example/MixedCurlLaplacianExample/Example_ClosureModel_Factory.hpp create mode 100644 packages/panzer/adapters-stk/example/MixedCurlLaplacianExample/Example_ClosureModel_Factory_TemplateBuilder.hpp create mode 100644 packages/panzer/adapters-stk/example/MixedCurlLaplacianExample/Example_ClosureModel_Factory_impl.hpp create mode 100644 packages/panzer/adapters-stk/example/MixedCurlLaplacianExample/Example_CurlLaplacianEquationSet.hpp create mode 100644 packages/panzer/adapters-stk/example/MixedCurlLaplacianExample/Example_CurlLaplacianEquationSet_impl.hpp create mode 100644 packages/panzer/adapters-stk/example/MixedCurlLaplacianExample/Example_CurlSolution.hpp create mode 100644 packages/panzer/adapters-stk/example/MixedCurlLaplacianExample/Example_CurlSolution_impl.hpp create mode 100644 packages/panzer/adapters-stk/example/MixedCurlLaplacianExample/Example_EquationSetFactory.hpp create mode 100644 packages/panzer/adapters-stk/example/MixedCurlLaplacianExample/Example_SimpleSource.hpp create mode 100644 packages/panzer/adapters-stk/example/MixedCurlLaplacianExample/Example_SimpleSource_impl.hpp create mode 100644 packages/panzer/adapters-stk/example/MixedCurlLaplacianExample/convergence_rate.py create mode 100644 packages/panzer/adapters-stk/example/MixedCurlLaplacianExample/main.cpp diff --git a/packages/panzer/adapters-stk/example/CMakeLists.txt b/packages/panzer/adapters-stk/example/CMakeLists.txt index bd3ee7da746a..e39f1f84a9be 100644 --- a/packages/panzer/adapters-stk/example/CMakeLists.txt +++ b/packages/panzer/adapters-stk/example/CMakeLists.txt @@ -1,6 +1,7 @@ ADD_SUBDIRECTORY(square_mesh) ADD_SUBDIRECTORY(PoissonExample) ADD_SUBDIRECTORY(CurlLaplacianExample) +ADD_SUBDIRECTORY(MixedCurlLaplacianExample) ADD_SUBDIRECTORY(MixedPoissonExample) ADD_SUBDIRECTORY(PoissonInterfaceExample) ADD_SUBDIRECTORY(assembly_engine) diff --git a/packages/panzer/adapters-stk/example/MixedCurlLaplacianExample/CMakeLists.txt b/packages/panzer/adapters-stk/example/MixedCurlLaplacianExample/CMakeLists.txt new file mode 100644 index 000000000000..70c5be73db9e --- /dev/null +++ b/packages/panzer/adapters-stk/example/MixedCurlLaplacianExample/CMakeLists.txt @@ -0,0 +1,229 @@ + +INCLUDE_DIRECTORIES(${CMAKE_CURRENT_SOURCE_DIR}) + +SET(ASSEMBLY_EXAMPLE_SOURCES + main.cpp + ) + +TRIBITS_ADD_EXECUTABLE( + MixedCurlLaplacianExample + SOURCES ${ASSEMBLY_EXAMPLE_SOURCES} + ) + +TRIBITS_ADD_ADVANCED_TEST( + MixedCurlLaplacianExample + EXCLUDE_IF_NOT_TRUE ${${PARENT_PACKAGE_NAME}_ADD_EXPENSIVE_CUDA_TESTS} + TEST_0 EXEC MixedCurlLaplacianExample + ARGS --use-epetra --output-filename="base-curl-test-" + PASS_REGULAR_EXPRESSION "ALL PASSED: Epetra" + NUM_MPI_PROCS 4 + TEST_1 EXEC MixedCurlLaplacianExample + ARGS --use-tpetra --output-filename="base-curl-test-" + PASS_REGULAR_EXPRESSION "ALL PASSED: Tpetra" + NUM_MPI_PROCS 4 + COMM mpi + ) + +FOREACH( ORDER 1 2 ) + TRIBITS_ADD_ADVANCED_TEST( + MixedCurlLaplacianExample-ConvTest-Quad-Order-${ORDER} + EXCLUDE_IF_NOT_TRUE ${PARENT_PACKAGE_NAME}_ADD_EXPENSIVE_CUDA_TESTS + TEST_0 EXEC MixedCurlLaplacianExample + ARGS --use-tpetra --use-twod --cell="Quad" --x-elements=4 --y-elements=4 --z-elements=4 --basis-order=${ORDER} + PASS_REGULAR_EXPRESSION "ALL PASSED: Tpetra" + NUM_MPI_PROCS 4 + OUTPUT_FILE MPE-ConvTest-Quad-${ORDER}-04 + TEST_1 EXEC MixedCurlLaplacianExample + ARGS --use-tpetra --use-twod --cell="Quad" --x-elements=8 --y-elements=8 --z-elements=4 --basis-order=${ORDER} + PASS_REGULAR_EXPRESSION "ALL PASSED: Tpetra" + NUM_MPI_PROCS 4 + OUTPUT_FILE MPE-ConvTest-Quad-${ORDER}-08 + TEST_2 EXEC MixedCurlLaplacianExample + ARGS --use-tpetra --use-twod --cell="Quad" --x-elements=16 --y-elements=16 --z-elements=4 --basis-order=${ORDER} + PASS_REGULAR_EXPRESSION "ALL PASSED: Tpetra" + NUM_MPI_PROCS 4 + OUTPUT_FILE MPE-ConvTest-Quad-${ORDER}-16 + TEST_3 EXEC MixedCurlLaplacianExample + ARGS --use-tpetra --use-twod --cell="Quad" --x-elements=32 --y-elements=32 --z-elements=4 --basis-order=${ORDER} + PASS_REGULAR_EXPRESSION "ALL PASSED: Tpetra" + NUM_MPI_PROCS 4 + OUTPUT_FILE MPE-ConvTest-Quad-${ORDER}-32 + TEST_4 CMND python + ARGS ${CMAKE_CURRENT_SOURCE_DIR}/convergence_rate.py + ${ORDER} + MPE-ConvTest-Quad-${ORDER}- + 4 + 8 + 16 + 32 + PASS_REGULAR_EXPRESSION "Test Passed" + COMM mpi + ) +ENDFOREACH() + +FOREACH( ORDER 3 ) + TRIBITS_ADD_ADVANCED_TEST( + MixedCurlLaplacianExample-ConvTest-Quad-Order-${ORDER} + EXCLUDE_IF_NOT_TRUE ${PARENT_PACKAGE_NAME}_ADD_EXPENSIVE_CUDA_TESTS + TEST_0 EXEC MixedCurlLaplacianExample + ARGS --use-tpetra --use-twod --cell="Quad" --x-elements=4 --y-elements=4 --z-elements=4 --basis-order=${ORDER} + PASS_REGULAR_EXPRESSION "ALL PASSED: Tpetra" + NUM_MPI_PROCS 4 + OUTPUT_FILE MPE-ConvTest-Quad-${ORDER}-04 + TEST_1 EXEC MixedCurlLaplacianExample + ARGS --use-tpetra --use-twod --cell="Quad" --x-elements=8 --y-elements=8 --z-elements=4 --basis-order=${ORDER} + PASS_REGULAR_EXPRESSION "ALL PASSED: Tpetra" + NUM_MPI_PROCS 4 + OUTPUT_FILE MPE-ConvTest-Quad-${ORDER}-08 + TEST_2 EXEC MixedCurlLaplacianExample + ARGS --use-tpetra --use-twod --cell="Quad" --x-elements=16 --y-elements=16 --z-elements=4 --basis-order=${ORDER} + PASS_REGULAR_EXPRESSION "ALL PASSED: Tpetra" + NUM_MPI_PROCS 4 + OUTPUT_FILE MPE-ConvTest-Quad-${ORDER}-16 + TEST_4 CMND python + ARGS ${CMAKE_CURRENT_SOURCE_DIR}/convergence_rate.py + ${ORDER} + MPE-ConvTest-Quad-${ORDER}- + 4 + 8 + 16 + PASS_REGULAR_EXPRESSION "Test Passed" + COMM mpi + ) +ENDFOREACH() + +FOREACH( ORDER 1 ) + TRIBITS_ADD_ADVANCED_TEST( + MixedCurlLaplacianMultiblockExample-ConvTest-Quad-Order-${ORDER} + EXCLUDE_IF_NOT_TRUE ${PARENT_PACKAGE_NAME}_ADD_EXPENSIVE_CUDA_TESTS + TEST_0 EXEC MixedCurlLaplacianExample + ARGS --use-tpetra --use-twod --x-blocks=2 --cell="Quad" --x-elements=8 --y-elements=8 --z-elements=4 --basis-order=${ORDER} --output-filename="multiblock-" + PASS_REGULAR_EXPRESSION "ALL PASSED: Tpetra" + NUM_MPI_PROCS 4 + OUTPUT_FILE MPE-Multiblock-ConvTest-Quad-${ORDER}-08 + TEST_1 EXEC MixedCurlLaplacianExample + ARGS --use-tpetra --use-twod --x-blocks=2 --cell="Quad" --x-elements=16 --y-elements=16 --z-elements=4 --basis-order=${ORDER} --output-filename="multiblock-" + PASS_REGULAR_EXPRESSION "ALL PASSED: Tpetra" + NUM_MPI_PROCS 4 + OUTPUT_FILE MPE-Multiblock-ConvTest-Quad-${ORDER}-16 + TEST_2 EXEC MixedCurlLaplacianExample + ARGS --use-tpetra --use-twod --x-blocks=2 --cell="Quad" --x-elements=32 --y-elements=32 --z-elements=4 --basis-order=${ORDER} --output-filename="multiblock-" + PASS_REGULAR_EXPRESSION "ALL PASSED: Tpetra" + NUM_MPI_PROCS 4 + OUTPUT_FILE MPE-Multiblock-ConvTest-Quad-${ORDER}-32 + TEST_3 CMND python + ARGS ${CMAKE_CURRENT_SOURCE_DIR}/convergence_rate.py + ${ORDER} + MPE-Multiblock-ConvTest-Quad-${ORDER}- + 8 + 16 + 32 + PASS_REGULAR_EXPRESSION "Test Passed" + COMM mpi + ) +ENDFOREACH() + +FOREACH( ORDER 1 ) + TRIBITS_ADD_ADVANCED_TEST( + MixedCurlLaplacianExample-ConvTest-Tri-Order-${ORDER} + EXCLUDE_IF_NOT_TRUE ${PARENT_PACKAGE_NAME}_ADD_EXPENSIVE_CUDA_TESTS + TEST_0 EXEC MixedCurlLaplacianExample + ARGS --use-tpetra --use-twod --cell="Tri" --x-elements=8 --y-elements=8 --z-elements=4 --basis-order=${ORDER} + PASS_REGULAR_EXPRESSION "ALL PASSED: Tpetra" + NUM_MPI_PROCS 4 + OUTPUT_FILE MPE-ConvTest-Tri-${ORDER}-08 + TEST_1 EXEC MixedCurlLaplacianExample + ARGS --use-tpetra --use-twod --cell="Tri" --x-elements=16 --y-elements=16 --z-elements=4 --basis-order=${ORDER} + PASS_REGULAR_EXPRESSION "ALL PASSED: Tpetra" + NUM_MPI_PROCS 4 + OUTPUT_FILE MPE-ConvTest-Tri-${ORDER}-16 + TEST_2 EXEC MixedCurlLaplacianExample + ARGS --use-tpetra --use-twod --cell="Tri" --x-elements=32 --y-elements=32 --z-elements=4 --basis-order=${ORDER} + PASS_REGULAR_EXPRESSION "ALL PASSED: Tpetra" + NUM_MPI_PROCS 4 + OUTPUT_FILE MPE-ConvTest-Tri-${ORDER}-32 + TEST_3 EXEC MixedCurlLaplacianExample + ARGS --use-tpetra --use-twod --cell="Tri" --x-elements=64 --y-elements=64 --z-elements=4 --basis-order=${ORDER} + PASS_REGULAR_EXPRESSION "ALL PASSED: Tpetra" + NUM_MPI_PROCS 4 + OUTPUT_FILE MPE-ConvTest-Tri-${ORDER}-64 + TEST_4 CMND python + ARGS ${CMAKE_CURRENT_SOURCE_DIR}/convergence_rate.py + ${ORDER} + MPE-ConvTest-Tri-${ORDER}- + 8 + 16 + 32 + 64 + PASS_REGULAR_EXPRESSION "Test Passed" + COMM mpi + ) +ENDFOREACH() + +FOREACH( ORDER 2 ) + TRIBITS_ADD_ADVANCED_TEST( + MixedCurlLaplacianExample-ConvTest-Tri-Order-${ORDER} + EXCLUDE_IF_NOT_TRUE ${PARENT_PACKAGE_NAME}_ADD_EXPENSIVE_CUDA_TESTS + TEST_0 EXEC MixedCurlLaplacianExample + ARGS --use-tpetra --use-twod --cell="Tri" --x-elements=4 --y-elements=4 --z-elements=4 --basis-order=${ORDER} + PASS_REGULAR_EXPRESSION "ALL PASSED: Tpetra" + NUM_MPI_PROCS 4 + OUTPUT_FILE MPE-ConvTest-Tri-${ORDER}-04 + TEST_1 EXEC MixedCurlLaplacianExample + ARGS --use-tpetra --use-twod --cell="Tri" --x-elements=8 --y-elements=8 --z-elements=4 --basis-order=${ORDER} + PASS_REGULAR_EXPRESSION "ALL PASSED: Tpetra" + NUM_MPI_PROCS 4 + OUTPUT_FILE MPE-ConvTest-Tri-${ORDER}-08 + TEST_2 EXEC MixedCurlLaplacianExample + ARGS --use-tpetra --use-twod --cell="Tri" --x-elements=16 --y-elements=16 --z-elements=4 --basis-order=${ORDER} + PASS_REGULAR_EXPRESSION "ALL PASSED: Tpetra" + NUM_MPI_PROCS 4 + OUTPUT_FILE MPE-ConvTest-Tri-${ORDER}-16 + TEST_3 CMND python + ARGS ${CMAKE_CURRENT_SOURCE_DIR}/convergence_rate.py + ${ORDER} + MPE-ConvTest-Tri-${ORDER}- + 4 + 8 + 16 + PASS_REGULAR_EXPRESSION "Test Passed" + COMM mpi + ) +ENDFOREACH() + +# FOREACH( ORDER 1 2 3 4 ) +# TRIBITS_ADD_ADVANCED_TEST( +# CurlLaplacianExample-ConvTest-Hex-Order-${ORDER} +# EXCLUDE_IF_NOT_TRUE ${PARENT_PACKAGE_NAME}_ADD_EXPENSIVE_CUDA_TESTS +# TEST_0 EXEC CurlLaplacianExample +# ARGS --use-epetra --use-threed --x-elements=4 --y-elements=4 --z-elements=4 --basis-order=${ORDER} -- +# PASS_REGULAR_EXPRESSION "ALL PASSED: Epetra" +# NUM_MPI_PROCS 4 +# OUTPUT_FILE MPE-ConvTest-Hex-${ORDER}-04 +# TEST_1 EXEC CurlLaplacianExample +# ARGS --use-epetra --use-threed --x-elements=8 --y-elements=8 --z-elements=4 --basis-order=${ORDER} +# PASS_REGULAR_EXPRESSION "ALL PASSED: Epetra" +# NUM_MPI_PROCS 4 +# OUTPUT_FILE MPE-ConvTest-Hex-${ORDER}-08 +# TEST_2 EXEC CurlLaplacianExample +# ARGS --use-epetra --use-threed --x-elements=16 --y-elements=16 --z-elements=4 --basis-order=${ORDER} +# PASS_REGULAR_EXPRESSION "ALL PASSED: Epetra" +# NUM_MPI_PROCS 4 +# OUTPUT_FILE MPE-ConvTest-Hex-${ORDER}-16 +# TEST_3 EXEC CurlLaplacianExample +# ARGS --use-epetra --use-threed --x-elements=32 --y-elements=32 --z-elements=4 --basis-order=${ORDER} +# PASS_REGULAR_EXPRESSION "ALL PASSED: Epetra" +# NUM_MPI_PROCS 4 +# OUTPUT_FILE MPE-ConvTest-Hex-${ORDER}-32 +# TEST_4 CMND python +# ARGS ${CMAKE_CURRENT_SOURCE_DIR}/convergence_rate.py +# ${ORDER} +# MPE-ConvTest-Hex-${ORDER}- +# 4 +# 8 +# 16 +# 32 +# PASS_REGULAR_EXPRESSION "Test Passed" +# COMM mpi +# ) +# ENDFOREACH() diff --git a/packages/panzer/adapters-stk/example/MixedCurlLaplacianExample/Example_BCStrategy_Dirichlet_Constant.hpp b/packages/panzer/adapters-stk/example/MixedCurlLaplacianExample/Example_BCStrategy_Dirichlet_Constant.hpp new file mode 100644 index 000000000000..5e95296a7df4 --- /dev/null +++ b/packages/panzer/adapters-stk/example/MixedCurlLaplacianExample/Example_BCStrategy_Dirichlet_Constant.hpp @@ -0,0 +1,80 @@ +// @HEADER +// *********************************************************************** +// +// Panzer: A partial differential equation assembly +// engine for strongly coupled complex multiphysics systems +// Copyright (2011) Sandia Corporation +// +// Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation, +// the U.S. Government retains certain rights in this software. +// +// Redistribution and use in source and binary forms, with or without +// modification, are permitted provided that the following conditions are +// met: +// +// 1. Redistributions of source code must retain the above copyright +// notice, this list of conditions and the following disclaimer. +// +// 2. Redistributions in binary form must reproduce the above copyright +// notice, this list of conditions and the following disclaimer in the +// documentation and/or other materials provided with the distribution. +// +// 3. Neither the name of the Corporation nor the names of the +// contributors may be used to endorse or promote products derived from +// this software without specific prior written permission. +// +// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY +// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE +// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR +// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE +// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, +// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, +// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR +// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF +// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING +// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS +// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. +// +// Questions? Contact Roger P. Pawlowski (rppawlo@sandia.gov) and +// Eric C. Cyr (eccyr@sandia.gov) +// *********************************************************************** +// @HEADER + +#ifndef __Example_BC_Dirichlet_Constant_hpp__ +#define __Example_BC_Dirichlet_Constant_hpp__ + +#include +#include + +#include "Teuchos_RCP.hpp" +#include "Panzer_BCStrategy_Dirichlet_DefaultImpl.hpp" +#include "Panzer_Traits.hpp" +#include "Panzer_PureBasis.hpp" +#include "Phalanx_FieldManager.hpp" + +namespace Example { + + template + class BCStrategy_Dirichlet_Constant : public panzer::BCStrategy_Dirichlet_DefaultImpl { + public: + + BCStrategy_Dirichlet_Constant(const panzer::BC& bc, const Teuchos::RCP& global_data); + + void setup(const panzer::PhysicsBlock& side_pb, + const Teuchos::ParameterList& user_data); + + void buildAndRegisterEvaluators(PHX::FieldManager& fm, + const panzer::PhysicsBlock& pb, + const panzer::ClosureModelFactory_TemplateManager& factory, + const Teuchos::ParameterList& models, + const Teuchos::ParameterList& user_data) const; + + std::string residual_name; + Teuchos::RCP basis; + }; + +} + +#include "Example_BCStrategy_Dirichlet_Constant_impl.hpp" + +#endif diff --git a/packages/panzer/adapters-stk/example/MixedCurlLaplacianExample/Example_BCStrategy_Dirichlet_Constant_impl.hpp b/packages/panzer/adapters-stk/example/MixedCurlLaplacianExample/Example_BCStrategy_Dirichlet_Constant_impl.hpp new file mode 100644 index 000000000000..4fe16111034c --- /dev/null +++ b/packages/panzer/adapters-stk/example/MixedCurlLaplacianExample/Example_BCStrategy_Dirichlet_Constant_impl.hpp @@ -0,0 +1,144 @@ +// @HEADER +// *********************************************************************** +// +// Panzer: A partial differential equation assembly +// engine for strongly coupled complex multiphysics systems +// Copyright (2011) Sandia Corporation +// +// Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation, +// the U.S. Government retains certain rights in this software. +// +// Redistribution and use in source and binary forms, with or without +// modification, are permitted provided that the following conditions are +// met: +// +// 1. Redistributions of source code must retain the above copyright +// notice, this list of conditions and the following disclaimer. +// +// 2. Redistributions in binary form must reproduce the above copyright +// notice, this list of conditions and the following disclaimer in the +// documentation and/or other materials provided with the distribution. +// +// 3. Neither the name of the Corporation nor the names of the +// contributors may be used to endorse or promote products derived from +// this software without specific prior written permission. +// +// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY +// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE +// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR +// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE +// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, +// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, +// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR +// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF +// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING +// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS +// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. +// +// Questions? Contact Roger P. Pawlowski (rppawlo@sandia.gov) and +// Eric C. Cyr (eccyr@sandia.gov) +// *********************************************************************** +// @HEADER + +#include "Teuchos_ParameterList.hpp" +#include "Teuchos_RCP.hpp" +#include "Teuchos_Assert.hpp" +#include "Phalanx_DataLayout_MDALayout.hpp" +#include "Phalanx_FieldManager.hpp" +#include "Panzer_PhysicsBlock.hpp" + +#include "Panzer_PureBasis.hpp" + +// Evaluators +#include "Panzer_Constant.hpp" + +#include "Phalanx_MDField.hpp" +#include "Phalanx_DataLayout.hpp" +#include "Phalanx_DataLayout_MDALayout.hpp" + +// *********************************************************************** +template +Example::BCStrategy_Dirichlet_Constant:: +BCStrategy_Dirichlet_Constant(const panzer::BC& bc, const Teuchos::RCP& global_data) : + panzer::BCStrategy_Dirichlet_DefaultImpl(bc,global_data) +{ + TEUCHOS_ASSERT(this->m_bc.strategy() == "Constant"); +} + +// *********************************************************************** +template +void Example::BCStrategy_Dirichlet_Constant:: +setup(const panzer::PhysicsBlock& side_pb, + const Teuchos::ParameterList& /* user_data */) +{ + using Teuchos::RCP; + using std::vector; + using std::string; + using std::pair; + + // need the dof value to form the residual + this->required_dof_names.push_back(this->m_bc.equationSetName()); + + // unique residual name + this->residual_name = "Residual_" + this->m_bc.identifier(); + + // map residual to dof + this->residual_to_dof_names_map[residual_name] = this->m_bc.equationSetName(); + + // map residual to target field + this->residual_to_target_field_map[residual_name] = "Constant_" + this->m_bc.equationSetName(); + + // find the basis for this dof + const vector > >& dofs = side_pb.getProvidedDOFs(); + + for (vector > >::const_iterator dof_it = + dofs.begin(); dof_it != dofs.end(); ++dof_it) { + if (dof_it->first == this->m_bc.equationSetName()) + this->basis = dof_it->second; + } + + TEUCHOS_TEST_FOR_EXCEPTION(Teuchos::is_null(this->basis), std::runtime_error, + "Error the name \"" << this->m_bc.equationSetName() + << "\" is not a valid DOF for the boundary condition:\n" + << this->m_bc << "\n"); + +} + +// *********************************************************************** +template +void Example::BCStrategy_Dirichlet_Constant:: +buildAndRegisterEvaluators(PHX::FieldManager& fm, + const panzer::PhysicsBlock& /* pb */, + const panzer::ClosureModelFactory_TemplateManager& /* factory */, + const Teuchos::ParameterList& /* models */, + const Teuchos::ParameterList& /* user_data */) const +{ + using Teuchos::ParameterList; + using Teuchos::RCP; + using Teuchos::rcp; + + // provide a constant target value to map into residual + if(basis->isScalarBasis()) { + ParameterList p("BC Constant Dirichlet"); + p.set("Name", "Constant_" + this->m_bc.equationSetName()); + p.set("Data Layout", basis->functional); + p.set("Value", this->m_bc.params()->template get("Value")); + + RCP< PHX::Evaluator > op = + rcp(new panzer::Constant(p)); + + this->template registerEvaluator(fm, op); + } + else if(basis->isVectorBasis()) { + ParameterList p("BC Constant Dirichlet"); + p.set("Name", "Constant_" + this->m_bc.equationSetName()); + p.set("Data Layout", basis->functional_grad); + p.set("Value", this->m_bc.params()->template get("Value")); + + RCP< PHX::Evaluator > op = + rcp(new panzer::Constant(p)); + + this->template registerEvaluator(fm, op); + } + +} diff --git a/packages/panzer/adapters-stk/example/MixedCurlLaplacianExample/Example_BCStrategy_Factory.hpp b/packages/panzer/adapters-stk/example/MixedCurlLaplacianExample/Example_BCStrategy_Factory.hpp new file mode 100644 index 000000000000..53e4385ea04d --- /dev/null +++ b/packages/panzer/adapters-stk/example/MixedCurlLaplacianExample/Example_BCStrategy_Factory.hpp @@ -0,0 +1,87 @@ +// @HEADER +// *********************************************************************** +// +// Panzer: A partial differential equation assembly +// engine for strongly coupled complex multiphysics systems +// Copyright (2011) Sandia Corporation +// +// Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation, +// the U.S. Government retains certain rights in this software. +// +// Redistribution and use in source and binary forms, with or without +// modification, are permitted provided that the following conditions are +// met: +// +// 1. Redistributions of source code must retain the above copyright +// notice, this list of conditions and the following disclaimer. +// +// 2. Redistributions in binary form must reproduce the above copyright +// notice, this list of conditions and the following disclaimer in the +// documentation and/or other materials provided with the distribution. +// +// 3. Neither the name of the Corporation nor the names of the +// contributors may be used to endorse or promote products derived from +// this software without specific prior written permission. +// +// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY +// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE +// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR +// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE +// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, +// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, +// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR +// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF +// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING +// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS +// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. +// +// Questions? Contact Roger P. Pawlowski (rppawlo@sandia.gov) and +// Eric C. Cyr (eccyr@sandia.gov) +// *********************************************************************** +// @HEADER + +#ifndef __Example_BCStrategyFactory_hpp__ +#define __Example_BCStrategyFactory_hpp__ + +#include "Teuchos_RCP.hpp" +#include "Panzer_Traits.hpp" +#include "Panzer_BCStrategy_TemplateManager.hpp" +#include "Panzer_BCStrategy_Factory.hpp" +#include "Panzer_BCStrategy_Factory_Defines.hpp" + +// Add my bcstrategies here +#include "Example_BCStrategy_Dirichlet_Constant.hpp" + +namespace Example { + +PANZER_DECLARE_BCSTRATEGY_TEMPLATE_BUILDER(BCStrategy_Dirichlet_Constant, + BCStrategy_Dirichlet_Constant) + +struct BCStrategyFactory : public panzer::BCStrategyFactory { + + Teuchos::RCP > + buildBCStrategy(const panzer::BC& bc,const Teuchos::RCP& global_data) const + { + + Teuchos::RCP > bcs_tm = + Teuchos::rcp(new panzer::BCStrategy_TemplateManager); + + bool found = false; + + PANZER_BUILD_BCSTRATEGY_OBJECTS("Constant", + BCStrategy_Dirichlet_Constant) + + TEUCHOS_TEST_FOR_EXCEPTION(!found, std::logic_error, + "Error - the BC Strategy called \"" << bc.strategy() << + "\" is not a valid identifier in the BCStrategyFactory. Either add a " + "valid implementation to your factory or fix your input file. The " + "relevant boundary condition is:\n\n" << bc << std::endl); + + return bcs_tm; + } + +}; + +} + +#endif diff --git a/packages/panzer/adapters-stk/example/MixedCurlLaplacianExample/Example_ClosureModel_Factory.hpp b/packages/panzer/adapters-stk/example/MixedCurlLaplacianExample/Example_ClosureModel_Factory.hpp new file mode 100644 index 000000000000..7e6f9aa531d5 --- /dev/null +++ b/packages/panzer/adapters-stk/example/MixedCurlLaplacianExample/Example_ClosureModel_Factory.hpp @@ -0,0 +1,71 @@ +// @HEADER +// *********************************************************************** +// +// Panzer: A partial differential equation assembly +// engine for strongly coupled complex multiphysics systems +// Copyright (2011) Sandia Corporation +// +// Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation, +// the U.S. Government retains certain rights in this software. +// +// Redistribution and use in source and binary forms, with or without +// modification, are permitted provided that the following conditions are +// met: +// +// 1. Redistributions of source code must retain the above copyright +// notice, this list of conditions and the following disclaimer. +// +// 2. Redistributions in binary form must reproduce the above copyright +// notice, this list of conditions and the following disclaimer in the +// documentation and/or other materials provided with the distribution. +// +// 3. Neither the name of the Corporation nor the names of the +// contributors may be used to endorse or promote products derived from +// this software without specific prior written permission. +// +// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY +// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE +// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR +// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE +// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, +// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, +// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR +// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF +// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING +// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS +// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. +// +// Questions? Contact Roger P. Pawlowski (rppawlo@sandia.gov) and +// Eric C. Cyr (eccyr@sandia.gov) +// *********************************************************************** +// @HEADER + +#ifndef __Example_ClosureModelFactory_hpp__ +#define __Example_ClosureModelFactory_hpp__ + +#include "Panzer_ClosureModel_Factory.hpp" + +namespace Example { + +template +class ModelFactory : public panzer::ClosureModelFactory { + +public: + + Teuchos::RCP< std::vector< Teuchos::RCP > > > + buildClosureModels(const std::string& model_id, + const Teuchos::ParameterList& models, + const panzer::FieldLayoutLibrary& fl, + const Teuchos::RCP& ir, + const Teuchos::ParameterList& default_params, + const Teuchos::ParameterList& user_data, + const Teuchos::RCP& global_data, + PHX::FieldManager& fm) const; + +}; + +} + +#include "Example_ClosureModel_Factory_impl.hpp" + +#endif diff --git a/packages/panzer/adapters-stk/example/MixedCurlLaplacianExample/Example_ClosureModel_Factory_TemplateBuilder.hpp b/packages/panzer/adapters-stk/example/MixedCurlLaplacianExample/Example_ClosureModel_Factory_TemplateBuilder.hpp new file mode 100644 index 000000000000..69bdf35aceb7 --- /dev/null +++ b/packages/panzer/adapters-stk/example/MixedCurlLaplacianExample/Example_ClosureModel_Factory_TemplateBuilder.hpp @@ -0,0 +1,66 @@ +// @HEADER +// *********************************************************************** +// +// Panzer: A partial differential equation assembly +// engine for strongly coupled complex multiphysics systems +// Copyright (2011) Sandia Corporation +// +// Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation, +// the U.S. Government retains certain rights in this software. +// +// Redistribution and use in source and binary forms, with or without +// modification, are permitted provided that the following conditions are +// met: +// +// 1. Redistributions of source code must retain the above copyright +// notice, this list of conditions and the following disclaimer. +// +// 2. Redistributions in binary form must reproduce the above copyright +// notice, this list of conditions and the following disclaimer in the +// documentation and/or other materials provided with the distribution. +// +// 3. Neither the name of the Corporation nor the names of the +// contributors may be used to endorse or promote products derived from +// this software without specific prior written permission. +// +// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY +// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE +// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR +// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE +// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, +// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, +// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR +// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF +// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING +// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS +// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. +// +// Questions? Contact Roger P. Pawlowski (rppawlo@sandia.gov) and +// Eric C. Cyr (eccyr@sandia.gov) +// *********************************************************************** +// @HEADER + +#ifndef __Example_ClosureModel_Factory_TemplateBuilder_hpp__ +#define __Example_ClosureModel_Factory_TemplateBuilder_hpp__ + +#include +#include "Sacado_mpl_apply.hpp" +#include "Teuchos_RCP.hpp" +#include "Example_ClosureModel_Factory.hpp" + +namespace Example { + +class ClosureModelFactory_TemplateBuilder { +public: + + template + Teuchos::RCP build() const + { + return Teuchos::rcp( static_cast(new Example::ModelFactory) ); + } + +}; + +} + +#endif diff --git a/packages/panzer/adapters-stk/example/MixedCurlLaplacianExample/Example_ClosureModel_Factory_impl.hpp b/packages/panzer/adapters-stk/example/MixedCurlLaplacianExample/Example_ClosureModel_Factory_impl.hpp new file mode 100644 index 000000000000..af8cc4de4a24 --- /dev/null +++ b/packages/panzer/adapters-stk/example/MixedCurlLaplacianExample/Example_ClosureModel_Factory_impl.hpp @@ -0,0 +1,322 @@ +// @HEADER +// *********************************************************************** +// +// Panzer: A partial differential equation assembly +// engine for strongly coupled complex multiphysics systems +// Copyright (2011) Sandia Corporation +// +// Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation, +// the U.S. Government retains certain rights in this software. +// +// Redistribution and use in source and binary forms, with or without +// modification, are permitted provided that the following conditions are +// met: +// +// 1. Redistributions of source code must retain the above copyright +// notice, this list of conditions and the following disclaimer. +// +// 2. Redistributions in binary form must reproduce the above copyright +// notice, this list of conditions and the following disclaimer in the +// documentation and/or other materials provided with the distribution. +// +// 3. Neither the name of the Corporation nor the names of the +// contributors may be used to endorse or promote products derived from +// this software without specific prior written permission. +// +// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY +// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE +// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR +// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE +// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, +// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, +// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR +// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF +// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING +// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS +// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. +// +// Questions? Contact Roger P. Pawlowski (rppawlo@sandia.gov) and +// Eric C. Cyr (eccyr@sandia.gov) +// *********************************************************************** +// @HEADER + +#ifndef __Example_ClosureModelFactoryT_hpp__ +#define __Example_ClosureModelFactoryT_hpp__ + +#include +#include +#include + +#include "Panzer_IntegrationRule.hpp" +#include "Panzer_BasisIRLayout.hpp" +#include "Panzer_ScalarToVector.hpp" +#include "Panzer_String_Utilities.hpp" +#include "Panzer_Sum.hpp" +#include "Panzer_DotProduct.hpp" +#include "Panzer_Product.hpp" + +#include "Phalanx_FieldTag_Tag.hpp" +#include "Teuchos_ParameterEntry.hpp" +#include "Teuchos_TypeNameTraits.hpp" + +#include "Example_SimpleSource.hpp" +#include "Example_CurlSolution.hpp" + +// ******************************************************************** +// ******************************************************************** +template +Teuchos::RCP< std::vector< Teuchos::RCP > > > +Example::ModelFactory:: +buildClosureModels(const std::string& model_id, + const Teuchos::ParameterList& models, + const panzer::FieldLayoutLibrary& fl, + const Teuchos::RCP& ir, + const Teuchos::ParameterList& /* default_params */, + const Teuchos::ParameterList& /* user_data */, + const Teuchos::RCP& /* global_data */, + PHX::FieldManager& /* fm */) const +{ + using std::string; + using std::vector; + using Teuchos::RCP; + using Teuchos::rcp; + using Teuchos::ParameterList; + using PHX::Evaluator; + + RCP< vector< RCP > > > evaluators = + rcp(new vector< RCP > > ); + + if (!models.isSublist(model_id)) { + models.print(std::cout); + std::stringstream msg; + msg << "Falied to find requested model, \"" << model_id + << "\", for equation set:\n" << std::endl; + TEUCHOS_TEST_FOR_EXCEPTION(!models.isSublist(model_id), std::logic_error, msg.str()); + } + + std::vector > bases; + fl.uniqueBases(bases); + + const ParameterList& my_models = models.sublist(model_id); + + for (ParameterList::ConstIterator model_it = my_models.begin(); + model_it != my_models.end(); ++model_it) { + + bool found = false; + + const std::string key = model_it->first; + ParameterList input; + const Teuchos::ParameterEntry& entry = model_it->second; + const ParameterList& plist = Teuchos::getValue(entry); + + if (plist.isType("Value")) { + // at IP + { + input.set("Name", key); + input.set("Value", plist.get("Value")); + input.set("Data Layout", ir->dl_scalar); + RCP< Evaluator > e = + rcp(new panzer::Constant(input)); + evaluators->push_back(e); + } + // at BASIS + for (std::vector >::const_iterator basis_itr = bases.begin(); + basis_itr != bases.end(); ++basis_itr) { + input.set("Name", key); + input.set("Value", plist.get("Value")); + Teuchos::RCP basis = basisIRLayout(*basis_itr,*ir); + input.set("Data Layout", basis->functional); + RCP< Evaluator > e = + rcp(new panzer::Constant(input)); + evaluators->push_back(e); + } + found = true; + } + + if (plist.isType("Scalar Names")) { + { // at IP + std::string value = plist.get("Scalar Names"); + Teuchos::RCP > scalarNames = Teuchos::rcp(new std::vector); + panzer::StringTokenizer(*scalarNames,value,",",true); + + input.set("Vector Name", key); + input.set("Scalar Names", scalarNames.getConst()); // evaluator expects a so make sure it is + input.set("Data Layout Scalar", ir->dl_scalar); + input.set("Data Layout Vector", ir->dl_vector); + RCP< Evaluator > e = + rcp(new panzer::ScalarToVector(input)); + evaluators->push_back(e); + } + found = true; + } + + if (plist.isType("Type")) { + std::string type = plist.get("Type"); + if(type=="SIMPLE SOURCE") { + RCP< Evaluator > e = + rcp(new Example::SimpleSource(key,*ir)); + evaluators->push_back(e); + + found = true; + } + else if(type=="EFIELD_EXACT") { + RCP< Evaluator > e = + rcp(new Example::CurlSolution(key,*ir)); + evaluators->push_back(e); + + found = true; + } + else if(type=="L2 ERROR_CALC") { + { + std::vector values(2); + values[0] = plist.get("Field A"); + values[1] = plist.get("Field B"); + + std::vector scalars(2); + scalars[0] = 1.0; + scalars[1] = -1.0; + + Teuchos::ParameterList p; + p.set("Sum Name",key+"_DIFF"); // Name of sum + p.set > >("Values Names",Teuchos::rcpFromRef(values)); + p.set > >("Scalars",Teuchos::rcpFromRef(scalars)); + p.set("Data Layout",ir->dl_vector); + + RCP< Evaluator > e = + rcp(new panzer::Sum(p)); + + evaluators->push_back(e); + } + + { + Teuchos::RCP pr = ir; + + Teuchos::ParameterList p; + p.set("Result Name",key); + p.set("Vector A Name",key+"_DIFF"); + p.set("Vector B Name",key+"_DIFF"); + p.set("Point Rule",pr); + + RCP< Evaluator > e = + rcp(new panzer::DotProduct(p)); + + evaluators->push_back(e); + } + + found = true; + } + else if(type=="HCurl ERROR_CALC") { + // Compute L2 contribution + { + std::vector values(2); + values[0] = plist.get("Field A"); + values[1] = plist.get("Field B"); + + std::vector scalars(2); + scalars[0] = 1.0; + scalars[1] = -1.0; + + Teuchos::ParameterList p; + p.set("Sum Name",key+"_HCurl_ValueDiff"); // Name of sum + p.set > >("Values Names",Teuchos::rcpFromRef(values)); + p.set > >("Scalars",Teuchos::rcpFromRef(scalars)); + p.set("Data Layout",ir->dl_vector); + + RCP< Evaluator > e = + rcp(new panzer::Sum(p)); + + evaluators->push_back(e); + } + + { + Teuchos::RCP pr = ir; + + Teuchos::ParameterList p; + p.set("Result Name",key+"_L2"); + p.set("Vector A Name",key+"_HCurl_ValueDiff"); + p.set("Vector B Name",key+"_HCurl_ValueDiff"); + p.set("Point Rule",pr); + + RCP< Evaluator > e = + rcp(new panzer::DotProduct(p)); + + evaluators->push_back(e); + } + + // Compute HCurl contribution + { + std::vector values(2); + values[0] = "CURL_" + plist.get("Field A"); + values[1] = "CURL_" + plist.get("Field B"); + + std::vector scalars(2); + scalars[0] = 1.0; + scalars[1] = -1.0; + + Teuchos::ParameterList p; + p.set("Sum Name",key+"_HCurl_CurlDiff"); // Name of sum + p.set > >("Values Names",Teuchos::rcpFromRef(values)); + p.set > >("Scalars",Teuchos::rcpFromRef(scalars)); + p.set("Data Layout",ir->dl_scalar); + + RCP< Evaluator > e = + rcp(new panzer::Sum(p)); + + evaluators->push_back(e); + } + + { + Teuchos::RCP pr = ir; + std::vector values(2); + values[0] = key+"_HCurl_CurlDiff"; + values[1] = key+"_HCurl_CurlDiff"; + + Teuchos::ParameterList p; + p.set("Product Name",key+"_HCurlSemi"); + p.set("Values Names",Teuchos::rcpFromRef(values)); + p.set("Data Layout",ir->dl_scalar); + + RCP< Evaluator > e = + rcp(new panzer::Product(p)); + + evaluators->push_back(e); + } + + { + std::vector values(2); + values[0] = key+"_L2"; + values[1] = key+"_HCurlSemi"; + + Teuchos::ParameterList p; + p.set("Sum Name",key); // Name of sum + p.set > >("Values Names",Teuchos::rcpFromRef(values)); + p.set("Data Layout",ir->dl_scalar); + + RCP< Evaluator > e = + rcp(new panzer::Sum(p)); + + evaluators->push_back(e); + } + + found = true; + } + + + // none found! + } + + + if (!found) { + std::stringstream msg; + msg << "ClosureModelFactory failed to build evaluator for key \"" << key + << "\"\nin model \"" << model_id + << "\". Please correct the type or add support to the \nfactory." < +#include + +#include "Teuchos_RCP.hpp" +#include "Panzer_EquationSet_DefaultImpl.hpp" +#include "Panzer_Traits.hpp" +#include "Phalanx_FieldManager.hpp" + +namespace Example { + +/** The equation set serves two roles. The first is to let the panzer library + * know which fields this equation set defines and their names. It registers + * the evaluators required for a particular equation set. The level of the + * granularity is largely up to a user. For instance this could be the momentum + * or continuity equation in Navier-Stokes, or it could simply be the Navier-Stokes + * equations. + * + * Generally, this inherits from the panzer::EquationSet_DefaultImpl which takes + * care of adding the gather (extract basis coefficients from solution vector) and + * scatter (using element matrices and vectors distribute and sum their values + * to a global linear system) evaluators. These use data members that must be set by + * the user. + */ +template +class CurlLaplacianEquationSet : public panzer::EquationSet_DefaultImpl { +public: + + /** In the constructor you set all the fields provided by this + * equation set. + */ + CurlLaplacianEquationSet(const Teuchos::RCP& params, + const int& default_integration_order, + const panzer::CellData& cell_data, + const Teuchos::RCP& global_data, + const bool build_transient_support); + + /** The specific evaluators are registered with the field manager argument. + */ + void buildAndRegisterEquationSetEvaluators(PHX::FieldManager& fm, + const panzer::FieldLibrary& field_library, + const Teuchos::ParameterList& user_data) const; + +}; + +} + +#include "Example_CurlLaplacianEquationSet_impl.hpp" + +#endif diff --git a/packages/panzer/adapters-stk/example/MixedCurlLaplacianExample/Example_CurlLaplacianEquationSet_impl.hpp b/packages/panzer/adapters-stk/example/MixedCurlLaplacianExample/Example_CurlLaplacianEquationSet_impl.hpp new file mode 100644 index 000000000000..a41083d2fe11 --- /dev/null +++ b/packages/panzer/adapters-stk/example/MixedCurlLaplacianExample/Example_CurlLaplacianEquationSet_impl.hpp @@ -0,0 +1,321 @@ +// @HEADER +// *********************************************************************** +// +// Panzer: A partial differential equation assembly +// engine for strongly coupled complex multiphysics systems +// Copyright (2011) Sandia Corporation +// +// Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation, +// the U.S. Government retains certain rights in this software. +// +// Redistribution and use in source and binary forms, with or without +// modification, are permitted provided that the following conditions are +// met: +// +// 1. Redistributions of source code must retain the above copyright +// notice, this list of conditions and the following disclaimer. +// +// 2. Redistributions in binary form must reproduce the above copyright +// notice, this list of conditions and the following disclaimer in the +// documentation and/or other materials provided with the distribution. +// +// 3. Neither the name of the Corporation nor the names of the +// contributors may be used to endorse or promote products derived from +// this software without specific prior written permission. +// +// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY +// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE +// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR +// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE +// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, +// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, +// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR +// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF +// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING +// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS +// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. +// +// Questions? Contact Roger P. Pawlowski (rppawlo@sandia.gov) and +// Eric C. Cyr (eccyr@sandia.gov) +// *********************************************************************** +// @HEADER + +#ifndef USER_APP_EQUATIONSET_ENERGY_T_HPP +#define USER_APP_EQUATIONSET_ENERGY_T_HPP + +#include "Teuchos_ParameterList.hpp" +#include "Teuchos_StandardParameterEntryValidators.hpp" +#include "Teuchos_RCP.hpp" +#include "Teuchos_Assert.hpp" +#include "Phalanx_DataLayout_MDALayout.hpp" +#include "Phalanx_FieldManager.hpp" + +#include "Panzer_IntegrationRule.hpp" +#include "Panzer_BasisIRLayout.hpp" + +// include evaluators here +#include "Panzer_Integrator_BasisTimesVector.hpp" +#include "Panzer_Integrator_BasisTimesScalar.hpp" +#include "Panzer_Integrator_CurlBasisDotVector.hpp" +#include "Panzer_ScalarToVector.hpp" +#include "Panzer_Sum.hpp" +#include "Panzer_Constant.hpp" + +// *********************************************************************** +template +Example::CurlLaplacianEquationSet:: +CurlLaplacianEquationSet(const Teuchos::RCP& params, + const int& default_integration_order, + const panzer::CellData& cell_data, + const Teuchos::RCP& global_data, + const bool build_transient_support) : + panzer::EquationSet_DefaultImpl(params, default_integration_order, cell_data, global_data, build_transient_support ) +{ + // ******************** + // Validate and parse parameter list + // ******************** + { + Teuchos::ParameterList valid_parameters; + this->setDefaultValidParameters(valid_parameters); + + valid_parameters.set("Model ID","","Closure model id associated with this equaiton set"); + valid_parameters.set("EField Basis Type","HCurl","Type of Basis to use"); + valid_parameters.set("BField Basis Type","HDiv","Type of Basis to use"); + valid_parameters.set("EField Basis Order",1,"Order of the basis"); + valid_parameters.set("BField Basis Order",1,"Order of the basis"); + valid_parameters.set("Integration Order",-1,"Order of the integration rule"); + + params->validateParametersAndSetDefaults(valid_parameters); + } + + std::string e_basis_type = params->get("EField Basis Type"); + std::string b_basis_type = params->get("BField Basis Type"); + int e_basis_order = params->get("EField Basis Order"); + int b_basis_order = params->get("BField Basis Order"); + int integration_order = params->get("Integration Order"); + std::string model_id = params->get("Model ID"); + + // ******************** + // Panzer uses strings to match fields. In this section we define the + // name of the fields provided by this equation set. This is a bit strange + // in that this is not the fields necessarily required by this equation set. + // For instance for the momentum equations in Navier-Stokes only the velocity + // fields are added, the pressure field is added by continuity. + // + // In this case "EFIELD" is the lone field. We also name the curl + // for this field. These names automatically generate evaluators for "EFIELD" + // and "CURL_EFIELD" gathering the basis coefficients of "EFIELD" and + // the values of the EFIELD and CURL_EFIELD fields at quadrature points. + // + // After all the equation set evaluators are added to a given field manager, the + // panzer code adds in appropriate scatter evaluators to distribute the + // entries into the residual and the Jacobian operator. These operators will be + // "required" by the field manager and will serve as roots of evaluation tree. + // The leaves of this tree will include the gather evaluators whose job it is to + // gather the solution from a vector. + // ******************** + + // ******************** + // Assemble DOF names and Residual names + // ******************** + + this->addDOF("EFIELD",e_basis_type,e_basis_order,integration_order); + this->addDOFCurl("EFIELD"); + if (this->buildTransientSupport()) + this->addDOFTimeDerivative("EFIELD"); + + this->addDOF("BFIELD",b_basis_type,b_basis_order,integration_order); + if (this->buildTransientSupport()) + this->addDOFTimeDerivative("BFIELD"); + + // ******************** + // Build Basis Functions and Integration Rules + // ******************** + + this->addClosureModel(model_id); + + this->setupDOFs(); +} + +// *********************************************************************** +template +void Example::CurlLaplacianEquationSet:: +buildAndRegisterEquationSetEvaluators(PHX::FieldManager& fm, + const panzer::FieldLibrary& /* fl */, + const Teuchos::ParameterList& /* user_data */) const +{ + using Teuchos::ParameterList; + using Teuchos::RCP; + using Teuchos::rcp; + + Teuchos::RCP e_ir = this->getIntRuleForDOF("EFIELD"); + Teuchos::RCP b_ir = this->getIntRuleForDOF("BFIELD"); + Teuchos::RCP e_basis = this->getBasisIRLayoutForDOF("EFIELD"); + Teuchos::RCP b_basis = this->getBasisIRLayoutForDOF("BFIELD"); + bool scalar_b = b_basis->dimension() < 3; + + // ******************** + // Energy Equation + // ******************** + + // Transient Operator: Assembles \int \dot{T} v + if (this->buildTransientSupport()) { + ParameterList p("E Transient Residual"); + p.set("Residual Name", "RESIDUAL_EFIELD_TRANSIENT_OP"); // we are defining the name of this operator + p.set("Value Name", "DXDT_EFIELD"); // this field is constructed by the panzer library +// p.set("Test Field Name", "EFIELD"); + p.set("Basis", e_basis); + p.set("IR", e_ir); + p.set("Multiplier", 1.0); + + RCP< PHX::Evaluator > op = + rcp(new panzer::Integrator_BasisTimesVector(p)); + + this->template registerEvaluator(fm, op); + } + + // Transient Operator: Assembles \int \dot{T} v + if (this->buildTransientSupport()) { + ParameterList p("B Transient Residual"); + p.set("Residual Name", "RESIDUAL_BFIELD_TRANSIENT_OP"); // we are defining the name of this operator + p.set("Value Name", "DXDT_BFIELD"); // this field is constructed by the panzer library +// p.set("Test Field Name", "EFIELD"); + p.set("Basis", b_basis); + p.set("IR", b_ir); + p.set("Multiplier", 1.0); + + RCP< PHX::Evaluator > op = + rcp(new panzer::Integrator_BasisTimesVector(p)); + + this->template registerEvaluator(fm, op); + } + + // Diffusion Operator: Assembles \int \nabla T \cdot \nabla v + { + using panzer::EvaluatorStyle; + using panzer::Integrator_CurlBasisDotVector; + using panzer::Traits; + using PHX::Evaluator; + using std::string; + string resName("RESIDUAL_EFIELD"), valName("BFIELD"); + double thermalConductivity(1.0), multiplier(-thermalConductivity); + RCP> op = rcp(new + Integrator_CurlBasisDotVector(EvaluatorStyle::CONTRIBUTES, + resName, valName, *e_basis, *e_ir, multiplier)); + this->template registerEvaluator(fm, op); + } + + // Mass operator + { + ParameterList p("Source Residual"); + p.set("Residual Name", "RESIDUAL_EFIELD_MASS_OP"); + p.set("Value Name", "EFIELD"); + p.set("Basis", e_basis); + p.set("IR", e_ir); + p.set("Multiplier", -1.0); + + RCP< PHX::Evaluator > op = + rcp(new panzer::Integrator_BasisTimesVector(p)); + + this->template registerEvaluator(fm, op); + } + + // Source Operator + { + ParameterList p("Source Residual"); + p.set("Residual Name", "RESIDUAL_EFIELD_SOURCE_OP"); + p.set("Value Name", "SOURCE_EFIELD"); // this field must be provided by the closure model factory + // and specified by the user + p.set("Basis", e_basis); + p.set("IR", e_ir); + p.set("Multiplier", 1.0); + + RCP< PHX::Evaluator > op = + rcp(new panzer::Integrator_BasisTimesVector(p)); + + this->template registerEvaluator(fm, op); + } + + // Use a sum operator to form the overall residual for the equation + { + std::vector sum_names; + + // these are the names of the residual values to sum together + sum_names.push_back("RESIDUAL_EFIELD_MASS_OP"); + sum_names.push_back("RESIDUAL_EFIELD_SOURCE_OP"); + if (this->buildTransientSupport()) + sum_names.push_back("RESIDUAL_EFIELD_TRANSIENT_OP"); + + this->buildAndRegisterResidualSummationEvaluator(fm,"EFIELD",sum_names); + } + + // Diffusion Operator: Assembles \int \nabla T \cdot \nabla v + { + using panzer::EvaluatorStyle; + using panzer::Integrator_BasisTimesVector; + using panzer::Integrator_BasisTimesScalar; + using panzer::Traits; + using PHX::Evaluator; + using std::string; + ParameterList p("B Curl Residual"); + p.set("Residual Name", "RESIDUAL_BFIELD_CURL_OP"); + p.set("Value Name", "CURL_EFIELD"); + p.set("Basis", b_basis); + p.set("IR", b_ir); + p.set("Multiplier", 1.0); + RCP> op; + if(scalar_b){ + op = rcp(new + Integrator_BasisTimesScalar(p)); + } else { + op = rcp(new + Integrator_BasisTimesVector(p)); + } + this->template registerEvaluator(fm, op); + } + + // Mass operator + { + using panzer::EvaluatorStyle; + using panzer::Integrator_BasisTimesVector; + using panzer::Integrator_BasisTimesScalar; + using panzer::Traits; + using PHX::Evaluator; + using std::string; + ParameterList p("Source Residual"); + p.set("Residual Name", "RESIDUAL_BFIELD_MASS_OP"); + p.set("Value Name", "BFIELD"); + p.set("Basis", b_basis); + p.set("IR", b_ir); + p.set("Multiplier", -1.0); + + RCP< PHX::Evaluator > op; + if(scalar_b){ + op = rcp(new + Integrator_BasisTimesScalar(p)); + } else { + op = rcp(new + Integrator_BasisTimesVector(p)); + } + + this->template registerEvaluator(fm, op); + } + + // Use a sum operator to form the overall residual for the equation + { + std::vector sum_names; + + // these are the names of the residual values to sum together + sum_names.push_back("RESIDUAL_BFIELD_MASS_OP"); + sum_names.push_back("RESIDUAL_BFIELD_CURL_OP"); + if (this->buildTransientSupport()) + sum_names.push_back("RESIDUAL_BFIELD_TRANSIENT_OP"); + + this->buildAndRegisterResidualSummationEvaluator(fm,"BFIELD",sum_names); + } + +} + +// *********************************************************************** + +#endif diff --git a/packages/panzer/adapters-stk/example/MixedCurlLaplacianExample/Example_CurlSolution.hpp b/packages/panzer/adapters-stk/example/MixedCurlLaplacianExample/Example_CurlSolution.hpp new file mode 100644 index 000000000000..6174d7e111a4 --- /dev/null +++ b/packages/panzer/adapters-stk/example/MixedCurlLaplacianExample/Example_CurlSolution.hpp @@ -0,0 +1,94 @@ +// @HEADER +// *********************************************************************** +// +// Panzer: A partial differential equation assembly +// engine for strongly coupled complex multiphysics systems +// Copyright (2011) Sandia Corporation +// +// Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation, +// the U.S. Government retains certain rights in this software. +// +// Redistribution and use in source and binary forms, with or without +// modification, are permitted provided that the following conditions are +// met: +// +// 1. Redistributions of source code must retain the above copyright +// notice, this list of conditions and the following disclaimer. +// +// 2. Redistributions in binary form must reproduce the above copyright +// notice, this list of conditions and the following disclaimer in the +// documentation and/or other materials provided with the distribution. +// +// 3. Neither the name of the Corporation nor the names of the +// contributors may be used to endorse or promote products derived from +// this software without specific prior written permission. +// +// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY +// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE +// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR +// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE +// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, +// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, +// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR +// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF +// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING +// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS +// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. +// +// Questions? Contact Roger P. Pawlowski (rppawlo@sandia.gov) and +// Eric C. Cyr (eccyr@sandia.gov) +// *********************************************************************** +// @HEADER + +#ifndef __Example_CurlSolution_hpp__ +#define __Example_CurlSolution_hpp__ + +#include "PanzerAdaptersSTK_config.hpp" + +#include "Phalanx_Evaluator_WithBaseImpl.hpp" +#include "Phalanx_Evaluator_Derived.hpp" +#include "Phalanx_FieldManager.hpp" + +#include "Panzer_Dimension.hpp" +#include "Panzer_FieldLibrary.hpp" + +#include + +#include "Panzer_Evaluator_WithBaseImpl.hpp" + +namespace Example { + + using panzer::Cell; + using panzer::Point; + using panzer::Dim; + +/** The analytic solution to the mixed poisson equation for the sine source. + */ +template +class CurlSolution : public panzer::EvaluatorWithBaseImpl, + public PHX::EvaluatorDerived { + +public: + CurlSolution(const std::string & name, + const panzer::IntegrationRule & ir); + + void postRegistrationSetup(typename Traits::SetupData d, + PHX::FieldManager& fm); + + void evaluateFields(typename Traits::EvalData d); + + +private: + typedef typename EvalT::ScalarT ScalarT; + + // Simulation solution + PHX::MDField solution; + PHX::MDField solution_curl; + int ir_degree, ir_index; +}; + +} + +#include "Example_CurlSolution_impl.hpp" + +#endif diff --git a/packages/panzer/adapters-stk/example/MixedCurlLaplacianExample/Example_CurlSolution_impl.hpp b/packages/panzer/adapters-stk/example/MixedCurlLaplacianExample/Example_CurlSolution_impl.hpp new file mode 100644 index 000000000000..6bbea230436d --- /dev/null +++ b/packages/panzer/adapters-stk/example/MixedCurlLaplacianExample/Example_CurlSolution_impl.hpp @@ -0,0 +1,108 @@ +// @HEADER +// *********************************************************************** +// +// Panzer: A partial differential equation assembly +// engine for strongly coupled complex multiphysics systems +// Copyright (2011) Sandia Corporation +// +// Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation, +// the U.S. Government retains certain rights in this software. +// +// Redistribution and use in solution and binary forms, with or without +// modification, are permitted provided that the following conditions are +// met: +// +// 1. Redistributions of solution code must retain the above copyright +// notice, this list of conditions and the following disclaimer. +// +// 2. Redistributions in binary form must reproduce the above copyright +// notice, this list of conditions and the following disclaimer in the +// documentation and/or other materials provided with the distribution. +// +// 3. Neither the name of the Corporation nor the names of the +// contributors may be used to endorse or promote products derived from +// this software without specific prior written permission. +// +// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY +// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE +// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR +// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE +// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, +// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, +// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR +// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF +// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING +// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS +// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. +// +// Questions? Contact Roger P. Pawlowski (rppawlo@sandia.gov) and +// Eric C. Cyr (eccyr@sandia.gov) +// *********************************************************************** +// @HEADER + +#ifndef __Example_CurlSolution_impl_hpp__ +#define __Example_CurlSolution_impl_hpp__ + +#include + +#include "Panzer_BasisIRLayout.hpp" +#include "Panzer_Workset.hpp" +#include "Panzer_Workset_Utilities.hpp" + +namespace Example { + +//********************************************************************** +template +CurlSolution::CurlSolution(const std::string & name, + const panzer::IntegrationRule & ir) +{ + using Teuchos::RCP; + + Teuchos::RCP dl_vector = ir.dl_vector; + Teuchos::RCP dl_scalar = ir.dl_scalar; + ir_degree = ir.cubature_degree; + + solution = PHX::MDField(name, dl_vector); + solution_curl = PHX::MDField("CURL_"+name, dl_scalar); + + // only thing that works + TEUCHOS_ASSERT(ir.spatial_dimension==2); + + this->addEvaluatedField(solution); + this->addEvaluatedField(solution_curl); + + std::string n = "Curl Solution"; + this->setName(n); +} + +//********************************************************************** +template +void CurlSolution::postRegistrationSetup(typename Traits::SetupData sd, + PHX::FieldManager& /* fm */) +{ + ir_index = panzer::getIntegrationRuleIndex(ir_degree,(*sd.worksets_)[0], this->wda); +} + +//********************************************************************** +template +void CurlSolution::evaluateFields(typename Traits::EvalData workset) +{ + using panzer::index_t; + for (index_t cell = 0; cell < workset.num_cells; ++cell) { + for (int point = 0; point < solution.extent_int(1); ++point) { + + const double & x = this->wda(workset).int_rules[ir_index]->ip_coordinates(cell,point,0); + const double & y = this->wda(workset).int_rules[ir_index]->ip_coordinates(cell,point,1); + + solution(cell,point,0) = -(y-1.0)*y + cos(2.0*M_PI*x)*sin(2.0*M_PI*y); + solution(cell,point,1) = -(x-1.0)*x + sin(2.0*M_PI*x)*cos(2.0*M_PI*y); + + solution_curl(cell,point) = -2.0*x+2.0*y; + } + } +} + +//********************************************************************** +} + +#endif diff --git a/packages/panzer/adapters-stk/example/MixedCurlLaplacianExample/Example_EquationSetFactory.hpp b/packages/panzer/adapters-stk/example/MixedCurlLaplacianExample/Example_EquationSetFactory.hpp new file mode 100644 index 000000000000..1290ec65d9b1 --- /dev/null +++ b/packages/panzer/adapters-stk/example/MixedCurlLaplacianExample/Example_EquationSetFactory.hpp @@ -0,0 +1,93 @@ +// @HEADER +// *********************************************************************** +// +// Panzer: A partial differential equation assembly +// engine for strongly coupled complex multiphysics systems +// Copyright (2011) Sandia Corporation +// +// Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation, +// the U.S. Government retains certain rights in this software. +// +// Redistribution and use in source and binary forms, with or without +// modification, are permitted provided that the following conditions are +// met: +// +// 1. Redistributions of source code must retain the above copyright +// notice, this list of conditions and the following disclaimer. +// +// 2. Redistributions in binary form must reproduce the above copyright +// notice, this list of conditions and the following disclaimer in the +// documentation and/or other materials provided with the distribution. +// +// 3. Neither the name of the Corporation nor the names of the +// contributors may be used to endorse or promote products derived from +// this software without specific prior written permission. +// +// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY +// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE +// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR +// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE +// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, +// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, +// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR +// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF +// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING +// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS +// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. +// +// Questions? Contact Roger P. Pawlowski (rppawlo@sandia.gov) and +// Eric C. Cyr (eccyr@sandia.gov) +// *********************************************************************** +// @HEADER + +#ifndef __CurlLaplacianExample_EquationSetFactory_hpp__ +#define __CurlLaplacianExample_EquationSetFactory_hpp__ + +#include "Panzer_EquationSet_Factory.hpp" +#include "Panzer_EquationSet_Factory_Defines.hpp" +#include "Panzer_CellData.hpp" + +// Add my equation sets here +#include "Example_CurlLaplacianEquationSet.hpp" + +namespace Example { + +// A macro that defines a class to make construction of the equation sets easier +// - The equation set is constructed over a list of automatic differention types +PANZER_DECLARE_EQSET_TEMPLATE_BUILDER(CurlLaplacianEquationSet, CurlLaplacianEquationSet) + +// A user written factory that creates each equation set. The key member here +// is buildEquationSet +class EquationSetFactory : public panzer::EquationSetFactory { +public: + + Teuchos::RCP > + buildEquationSet(const Teuchos::RCP& params, + const int& default_integration_order, + const panzer::CellData& cell_data, + const Teuchos::RCP& global_data, + const bool build_transient_support) const + { + Teuchos::RCP > eq_set= + Teuchos::rcp(new panzer::EquationSet_TemplateManager); + + bool found = false; // this is used by PANZER_BUILD_EQSET_OBJECTS + + // macro checks if(ies.name=="CurlLaplacian") then an EquationSet_Energy object is constructed + PANZER_BUILD_EQSET_OBJECTS("CurlLaplacian", CurlLaplacianEquationSet) + + // make sure your equation set has been found + if(!found) { + std::string msg = "Error - the \"Equation Set\" called \"" + params->get("Type") + + "\" is not a valid equation set identifier. Please supply the correct factory.\n"; + TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error, msg); + } + + return eq_set; + } + +}; + +} + +#endif diff --git a/packages/panzer/adapters-stk/example/MixedCurlLaplacianExample/Example_SimpleSource.hpp b/packages/panzer/adapters-stk/example/MixedCurlLaplacianExample/Example_SimpleSource.hpp new file mode 100644 index 000000000000..d753728bcd5c --- /dev/null +++ b/packages/panzer/adapters-stk/example/MixedCurlLaplacianExample/Example_SimpleSource.hpp @@ -0,0 +1,94 @@ +// @HEADER +// *********************************************************************** +// +// Panzer: A partial differential equation assembly +// engine for strongly coupled complex multiphysics systems +// Copyright (2011) Sandia Corporation +// +// Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation, +// the U.S. Government retains certain rights in this software. +// +// Redistribution and use in source and binary forms, with or without +// modification, are permitted provided that the following conditions are +// met: +// +// 1. Redistributions of source code must retain the above copyright +// notice, this list of conditions and the following disclaimer. +// +// 2. Redistributions in binary form must reproduce the above copyright +// notice, this list of conditions and the following disclaimer in the +// documentation and/or other materials provided with the distribution. +// +// 3. Neither the name of the Corporation nor the names of the +// contributors may be used to endorse or promote products derived from +// this software without specific prior written permission. +// +// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY +// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE +// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR +// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE +// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, +// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, +// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR +// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF +// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING +// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS +// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. +// +// Questions? Contact Roger P. Pawlowski (rppawlo@sandia.gov) and +// Eric C. Cyr (eccyr@sandia.gov) +// *********************************************************************** +// @HEADER + +#ifndef EXAMPLE_SIMPLE_SOURCE_DECL_HPP +#define EXAMPLE_SIMPLE_SOURCE_DECL_HPP + +#include "PanzerAdaptersSTK_config.hpp" + +#include "Phalanx_config.hpp" +#include "Phalanx_Evaluator_WithBaseImpl.hpp" +#include "Phalanx_Evaluator_Derived.hpp" +#include "Phalanx_FieldManager.hpp" + +#include "Panzer_Dimension.hpp" +#include "Panzer_FieldLibrary.hpp" + +#include + +#include "Panzer_Evaluator_WithBaseImpl.hpp" + +namespace Example { + + using panzer::Cell; + using panzer::Point; + using panzer::Dim; + +/** A source for the curl Laplacian that results in the solution + */ +template +class SimpleSource : public panzer::EvaluatorWithBaseImpl, + public PHX::EvaluatorDerived { + +public: + SimpleSource(const std::string & name, + const panzer::IntegrationRule & ir); + + void postRegistrationSetup(typename Traits::SetupData d, + PHX::FieldManager& fm); + + void evaluateFields(typename Traits::EvalData d); + + +private: + typedef typename EvalT::ScalarT ScalarT; + + // Simulation source + PHX::MDField source; + int ir_degree, ir_index; +}; + +} + +#include "Example_SimpleSource_impl.hpp" + +#endif diff --git a/packages/panzer/adapters-stk/example/MixedCurlLaplacianExample/Example_SimpleSource_impl.hpp b/packages/panzer/adapters-stk/example/MixedCurlLaplacianExample/Example_SimpleSource_impl.hpp new file mode 100644 index 000000000000..d9766a64a2f4 --- /dev/null +++ b/packages/panzer/adapters-stk/example/MixedCurlLaplacianExample/Example_SimpleSource_impl.hpp @@ -0,0 +1,104 @@ +// @HEADER +// *********************************************************************** +// +// Panzer: A partial differential equation assembly +// engine for strongly coupled complex multiphysics systems +// Copyright (2011) Sandia Corporation +// +// Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation, +// the U.S. Government retains certain rights in this software. +// +// Redistribution and use in source and binary forms, with or without +// modification, are permitted provided that the following conditions are +// met: +// +// 1. Redistributions of source code must retain the above copyright +// notice, this list of conditions and the following disclaimer. +// +// 2. Redistributions in binary form must reproduce the above copyright +// notice, this list of conditions and the following disclaimer in the +// documentation and/or other materials provided with the distribution. +// +// 3. Neither the name of the Corporation nor the names of the +// contributors may be used to endorse or promote products derived from +// this software without specific prior written permission. +// +// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY +// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE +// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR +// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE +// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, +// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, +// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR +// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF +// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING +// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS +// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. +// +// Questions? Contact Roger P. Pawlowski (rppawlo@sandia.gov) and +// Eric C. Cyr (eccyr@sandia.gov) +// *********************************************************************** +// @HEADER + +#ifndef EXAMPLE_SIMPLE_SOURCE_IMPL_HPP +#define EXAMPLE_SIMPLE_SOURCE_IMPL_HPP + +#include + +#include "Panzer_BasisIRLayout.hpp" +#include "Panzer_Workset.hpp" +#include "Panzer_Workset_Utilities.hpp" + +namespace Example { + +//********************************************************************** +template +SimpleSource::SimpleSource(const std::string & name, + const panzer::IntegrationRule & ir) +{ + using Teuchos::RCP; + + Teuchos::RCP data_layout = ir.dl_vector; + ir_degree = ir.cubature_degree; + + source = PHX::MDField(name, data_layout); + + this->addEvaluatedField(source); + + std::string n = "Simple Source"; + this->setName(n); +} + +//********************************************************************** +template +void SimpleSource::postRegistrationSetup(typename Traits::SetupData sd, + PHX::FieldManager& /* fm */) +{ + ir_index = panzer::getIntegrationRuleIndex(ir_degree,(*sd.worksets_)[0], this->wda); +} + +//********************************************************************** +template +void SimpleSource::evaluateFields(typename Traits::EvalData workset) +{ + using panzer::index_t; + for (index_t cell = 0; cell < workset.num_cells; ++cell) { + for (int point = 0; point < source.extent_int(1); ++point) { + + const double & x = this->wda(workset).int_rules[ir_index]->ip_coordinates(cell,point,0); + const double & y = this->wda(workset).int_rules[ir_index]->ip_coordinates(cell,point,1); + + source(cell,point,0) = 2.0+y-y*y + cos(2.0*M_PI*x)*sin(2.0*M_PI*y); + source(cell,point,1) = 2.0+x-x*x + sin(2.0*M_PI*x)*cos(2.0*M_PI*y); + + // if three d + if(source.extent(2)==3) + source(cell,point,2) = 0.0; + } + } +} + +//********************************************************************** +} + +#endif diff --git a/packages/panzer/adapters-stk/example/MixedCurlLaplacianExample/convergence_rate.py b/packages/panzer/adapters-stk/example/MixedCurlLaplacianExample/convergence_rate.py new file mode 100644 index 000000000000..b72bf1141c0f --- /dev/null +++ b/packages/panzer/adapters-stk/example/MixedCurlLaplacianExample/convergence_rate.py @@ -0,0 +1,87 @@ +import math +import sys +import re + +try: + if not len(sys.argv)>4: + print('Expected at least three solution files with the line \"Error = \%f\" in them') + raise Exception("Failure: Command line - [basis-order] [filename-prefix] [coarsest mesh res] [second mesh res] [third mesh res] ...") + + basis_order = float(sys.argv[1]) + file_prefix = sys.argv[2] + + mesh_res = [] + for i in range(3,len(sys.argv)): + mesh_res += [int(sys.argv[i])] + + max_res = max(mesh_res) + + #orderPat = re.compile('.*Basis Order = (.*)$') + l2ErrorPat = re.compile('.*L2 Error = (.*)$') + h1ErrorPat = re.compile('.*HCurl Error = (.*)$') + + #basis_order = 0 + l2_error_value = 0 + hcurl_error_value = 0 + prev_l2_error_value = 0 + prev_hcurl_error_value = 0 + + l2_error_values = [] + hcurl_error_values = [] + rate_values = [] + + for cnt, res in enumerate(mesh_res): + if max_res>=100: + filename = '%s%03d' % (file_prefix,res) + elif max_res>=10: + filename = '%s%02d' % (file_prefix,res) + else: + filename = '%s%d' % (file_prefix,res) + + print("opening \"" + filename +"\"") + f = open(filename); + + # look for the "Error = ..." line + for line in f: + #orderMatch = orderPat.match(line) + #if orderMatch!=None: + # basis_order = float(orderMatch.group(1)) + + errorMatch = l2ErrorPat.match(line) + if errorMatch!=None: + l2_error_value = float(errorMatch.group(1)) + l2_error_values += [l2_error_value] + + errorMatch = h1ErrorPat.match(line) + if errorMatch!=None: + hcurl_error_value = float(errorMatch.group(1)) + hcurl_error_values += [hcurl_error_value] + + if cnt>0: + prev_error = prev_l2_error_value + prev_hcurl_error_value + error = l2_error_value + hcurl_error_value + rate = math.log(prev_error/error)/math.log(2) + rate_values += [rate] + diff = rate - basis_order + + if rate > 0.90*basis_order: + print('%s: Convergence order of %f is within 10 percent of %f: PASSED' % (filename, rate, basis_order)) + else: + print('%s: Convergence order of %f is not within 10 percent of %f: FAILED' % (filename, rate, basis_order)) + raise Exception('Exception') + + prev_l2_error_value = l2_error_value + prev_hcurl_error_value = hcurl_error_value + + # end for res + + print('L2 errors : ', l2_error_values) + print('Curl errors : ', hcurl_error_values) + print('HCurl errors : ', [l2_error_values[i] + hcurl_error_values[i] for i in range(len(l2_error_values))]) + print('Conv rate : ', rate_values) + +except Exception as e: + print("Test Failed: ") + print(e) +else: + print("Test Passed") diff --git a/packages/panzer/adapters-stk/example/MixedCurlLaplacianExample/main.cpp b/packages/panzer/adapters-stk/example/MixedCurlLaplacianExample/main.cpp new file mode 100644 index 000000000000..0ef9d1f97e27 --- /dev/null +++ b/packages/panzer/adapters-stk/example/MixedCurlLaplacianExample/main.cpp @@ -0,0 +1,768 @@ +// @HEADER +// *********************************************************************** +// +// Panzer: A partial differential equation assembly +// engine for strongly coupled complex multiphysics systems +// Copyright (2011) Sandia Corporation +// +// Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation, +// the U.S. Government retains certain rights in this software. +// +// Redistribution and use in source and binary forms, with or without +// modification, are permitted provided that the following conditions are +// met: +// +// 1. Redistributions of source code must retain the above copyright +// notice, this list of conditions and the following disclaimer. +// +// 2. Redistributions in binary form must reproduce the above copyright +// notice, this list of conditions and the following disclaimer in the +// documentation and/or other materials provided with the distribution. +// +// 3. Neither the name of the Corporation nor the names of the +// contributors may be used to endorse or promote products derived from +// this software without specific prior written permission. +// +// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY +// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE +// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR +// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE +// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, +// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, +// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR +// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF +// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING +// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS +// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. +// +// Questions? Contact Roger P. Pawlowski (rppawlo@sandia.gov) and +// Eric C. Cyr (eccyr@sandia.gov) +// *********************************************************************** +// @HEADER + +#include +#include +#include +#include +#include +#include + +#include "Teuchos_DefaultComm.hpp" +#include "Teuchos_GlobalMPISession.hpp" +#include "Teuchos_CommandLineProcessor.hpp" + +#include "PanzerAdaptersSTK_config.hpp" +#include "Panzer_GlobalData.hpp" +#include "Panzer_Workset_Builder.hpp" +#include "Panzer_WorksetContainer.hpp" +#include "Panzer_AssemblyEngine.hpp" +#include "Panzer_AssemblyEngine_TemplateManager.hpp" +#include "Panzer_AssemblyEngine_TemplateBuilder.hpp" +#include "Panzer_LinearObjFactory.hpp" +#include "Panzer_BlockedEpetraLinearObjFactory.hpp" +#include "Panzer_TpetraLinearObjFactory.hpp" +#include "Panzer_DOFManagerFactory.hpp" +#include "Panzer_FieldManagerBuilder.hpp" +#include "Panzer_PureBasis.hpp" +#include "Panzer_GlobalData.hpp" +#include "Panzer_ResponseLibrary.hpp" +#include "Panzer_ResponseEvaluatorFactory_Functional.hpp" +#include "Panzer_Response_Functional.hpp" +#include "Panzer_CheckBCConsistency.hpp" + +#include "PanzerAdaptersSTK_config.hpp" +#include "Panzer_STK_WorksetFactory.hpp" +#include "Panzer_STKConnManager.hpp" +#include "Panzer_STK_Version.hpp" +#include "Panzer_STK_Interface.hpp" +#include "Panzer_STK_SquareQuadMeshFactory.hpp" +#include "Panzer_STK_SquareTriMeshFactory.hpp" +#include "Panzer_STK_CubeHexMeshFactory.hpp" +#include "Panzer_STK_SetupUtilities.hpp" +#include "Panzer_STK_Utilities.hpp" +#include "Panzer_STK_ResponseEvaluatorFactory_SolutionWriter.hpp" + +#include "Epetra_MpiComm.h" + +#include "EpetraExt_RowMatrixOut.h" +#include "EpetraExt_VectorOut.h" + +#include "BelosPseudoBlockGmresSolMgr.hpp" +#include "BelosTpetraAdapter.hpp" + +#include "Example_BCStrategy_Factory.hpp" +#include "Example_ClosureModel_Factory_TemplateBuilder.hpp" +#include "Example_EquationSetFactory.hpp" + +#include "AztecOO.h" + +#include + +using Teuchos::RCP; +using Teuchos::rcp; + + +// +// This is the Mathematica code used to generate this example. +// It also generates a plot of the vector field so it is clear +// what the solution is doing. +// +// Needs["VectorAnalysis`"] +// +// phi0[x_,y_]=(1-x)*(1-y) +// phi1[x_,y_]=x*(1-y) +// phi2[x_,y_]=x*y +// phi3[x_,y_]=y*(1-x) +// +// psi0[x_,y_]={1-y,0,0} +// psi1[x_,y_]={0,x,0} +// psi2[x_,y_]={y,0,0} +// psi3[x_,y_]={0,1-x,0} +// +// u[x_,y_]=phi2[x,y]*psi0[x,y]+phi3[x,y]*psi1[x,y]+phi0[x,y]*psi2[x,y]+phi1[x,y]*psi3[x,y] +// f[x_,y_]=u[x,y]+Curl[Curl[u[x,y],Cartesian[x,y,z]],Cartesian[x,y,z]] +// +// TwoDVec[g_]={g[[1]],g[[2]]} +// +// DotProduct[u[0.5,0],{1,0,0}] +// DotProduct[u[1,0.5],{0,1,0}] +// DotProduct[u[0.5,1],{1,0,0}] +// DotProduct[u[0,0.5],{0,1,0}] +// +// Out[118]= 0. +// Out[119]= 0. +// Out[120]= 0. +// Out[121]= 0. +// +// VectorPlot[TwoDVec[u[x,y]],{x,0,1},{y,0,1}] +// Simplify[u[x,y]] +// Simplify[f[x,y]] +// +// Out[144]= {-(-1+y) y,-(-1+x) x,0} +// Out[145]= {2+y-y^2,2+x-x^2,0} +// + +void testInitialization(const Teuchos::RCP& ipb, + std::vector& bcs, + const std::vector& eBlockNames, + int basis_order, + const bool threeD); + +void solveEpetraSystem(panzer::LinearObjContainer & container); +void solveTpetraSystem(panzer::LinearObjContainer & container); + +// calls MPI_Init and MPI_Finalize +int main(int argc,char * argv[]) +{ + using std::endl; + using Teuchos::RCP; + using Teuchos::rcp_dynamic_cast; + using panzer::StrPureBasisPair; + using panzer::StrPureBasisComp; + + Kokkos::initialize(argc,argv); + + Teuchos::GlobalMPISession mpiSession(&argc,&argv); + RCP Comm = Teuchos::rcp(new Epetra_MpiComm(MPI_COMM_WORLD)); + Teuchos::RCP > comm = Teuchos::rcp(new Teuchos::MpiComm(MPI_COMM_WORLD)); + Teuchos::FancyOStream out(Teuchos::rcpFromRef(std::cout)); + out.setOutputToRootOnly(0); + out.setShowProcRank(true); + + const auto stackedTimer = Teuchos::rcp(new Teuchos::StackedTimer("Panzer MixedPoisson Test")); + Teuchos::TimeMonitor::setStackedTimer(stackedTimer); + stackedTimer->start("Curl Laplacian"); + + // Build command line processor + //////////////////////////////////////////////////// + + bool useTpetra = true; + bool threeD = false; + int x_blocks = 1; + int x_elements=20,y_elements=20,z_elements=20; + std::string celltype = "Quad"; // or "Tri" (2d), Hex or Tet (3d) + double x_size=1.,y_size=1.,z_size=1.; + int basis_order = 1; + std::string output_filename="output_"; + + Teuchos::CommandLineProcessor clp; + clp.throwExceptions(false); + clp.setDocString("This example solves curl laplacian problem with Hex and Tet inline mesh with high order.\n"); + + clp.setOption("cell",&celltype); + clp.setOption("use-tpetra","use-epetra",&useTpetra); + clp.setOption("use-threed","use-twod",&threeD); + clp.setOption("x-blocks",&x_blocks); + clp.setOption("x-elements",&x_elements); + clp.setOption("y-elements",&y_elements); + clp.setOption("z-elements",&z_elements); + clp.setOption("x-size",&x_size); + clp.setOption("y-size",&y_size); + clp.setOption("z-size",&z_size); + clp.setOption("basis-order",&basis_order); + clp.setOption("output-filename",&output_filename); + + // parse commandline argument + Teuchos::CommandLineProcessor::EParseCommandLineReturn r_parse= clp.parse( argc, argv ); + if (r_parse == Teuchos::CommandLineProcessor::PARSE_HELP_PRINTED) return 0; + if (r_parse != Teuchos::CommandLineProcessor::PARSE_SUCCESSFUL ) return -1; + + // variable declarations + //////////////////////////////////////////////////// + + // factory definitions + Teuchos::RCP eqset_factory = + Teuchos::rcp(new Example::EquationSetFactory); // where poison equation is defined + Example::BCStrategyFactory bc_factory; // where boundary conditions are defined + + // construction of uncommitted (no elements) mesh + //////////////////////////////////////////////////////// + + RCP mesh_factory; + if(threeD) { + mesh_factory = rcp(new panzer_stk::CubeHexMeshFactory); + + // set mesh factory parameters + RCP pl = rcp(new Teuchos::ParameterList); + pl->set("X Blocks",x_blocks); + pl->set("Y Blocks",1); + pl->set("Z Blocks",1); + pl->set("X Elements",x_elements/x_blocks); + pl->set("Y Elements",y_elements); + pl->set("Z Elements",z_elements); + pl->set("Xf",x_size); + pl->set("Yf",y_size); + pl->set("Zf",z_size); + mesh_factory->setParameterList(pl); + } + else { + if (celltype == "Quad") mesh_factory = Teuchos::rcp(new panzer_stk::SquareQuadMeshFactory); + else if (celltype == "Tri") mesh_factory = Teuchos::rcp(new panzer_stk::SquareTriMeshFactory); + else + throw std::runtime_error("not supported celltype argument: try Quad or Tri"); + + // set mesh factory parameters + RCP pl = rcp(new Teuchos::ParameterList); + pl->set("X Blocks",x_blocks); + pl->set("Y Blocks",1); + pl->set("X Elements",x_elements/x_blocks); + pl->set("Y Elements",y_elements); + pl->set("Xf",x_size); + pl->set("Yf",y_size); + mesh_factory->setParameterList(pl); + } + + RCP mesh = mesh_factory->buildUncommitedMesh(MPI_COMM_WORLD); + + // other declarations + const std::size_t workset_size = 8; + + // construct input physics and physics block + //////////////////////////////////////////////////////// + + Teuchos::RCP ipb = Teuchos::parameterList("Physics Blocks"); + std::vector bcs; + std::vector > physicsBlocks; + { + bool build_transient_support = false; + + std::vector eBlockNames; + mesh->getElementBlockNames(eBlockNames); + testInitialization(ipb, bcs, eBlockNames, basis_order, threeD); + + const panzer::CellData volume_cell_data(workset_size, mesh->getCellTopology(eBlockNames[0])); + + // GobalData sets ostream and parameter interface to physics + Teuchos::RCP gd = panzer::createGlobalData(); + + // Can be overridden by the equation set + int default_integration_order = 4; + + // the physics block knows how to build and register evaluator with the field manager + for (const auto& block : eBlockNames) { + RCP pb + = rcp(new panzer::PhysicsBlock(ipb, block, + default_integration_order, + volume_cell_data, + eqset_factory, + gd, + build_transient_support)); + + // we can have more than one physics block, one per element block + physicsBlocks.push_back(pb); + } + } + + // finish building mesh, set required field variables and mesh bulk data + //////////////////////////////////////////////////////////////////////// + + for (const auto pb : physicsBlocks) { + const std::vector & blockFields = pb->getProvidedDOFs(); + + // insert all fields into a set + std::set fieldNames; + fieldNames.insert(blockFields.begin(),blockFields.end()); + + // build string for modifiying vectors + std::vector dimenStr(3); + dimenStr[0] = "X"; dimenStr[1] = "Y"; dimenStr[2] = "Z"; + + // add basis to DOF manager: block specific + std::set::const_iterator fieldItr; + for (fieldItr=fieldNames.begin();fieldItr!=fieldNames.end();++fieldItr) { + Teuchos::RCP basis = fieldItr->second; + if(basis->getElementSpace()==panzer::PureBasis::HGRAD) + mesh->addSolutionField(fieldItr->first,pb->elementBlockID()); + else if(basis->getElementSpace()==panzer::PureBasis::HCURL) { + for(int i=0;idimension();i++) + mesh->addCellField(fieldItr->first+dimenStr[i],pb->elementBlockID()); + } + else if(basis->getElementSpace()==panzer::PureBasis::HDIV) { + for(int i=0;idimension();i++) + mesh->addCellField(fieldItr->first+dimenStr[i],pb->elementBlockID()); + } + else if(basis->getElementSpace()==panzer::PureBasis::HVOL) { + for(int i=0;idimension();i++) + mesh->addCellField(fieldItr->first+dimenStr[i],pb->elementBlockID()); + } + } + } + mesh_factory->completeMeshConstruction(*mesh,MPI_COMM_WORLD); + + + // check that the bcs exist in the mesh + ///////////////////////////////////////////////////////////// + { + std::vector eBlockNames; + mesh->getElementBlockNames(eBlockNames); + std::vector sidesetNames; + mesh->getSidesetNames(sidesetNames); + panzer::checkBCConsistency(eBlockNames,sidesetNames,bcs); + } + + // build DOF Manager and linear object factory + ///////////////////////////////////////////////////////////// + + RCP dofManager; + Teuchos::RCP > linObjFactory; + + // build the connection manager + if(!useTpetra) { + const Teuchos::RCP > conn_manager = Teuchos::rcp(new panzer_stk::STKConnManager(mesh)); + + panzer::DOFManagerFactory globalIndexerFactory; + RCP > dofManager_int + = globalIndexerFactory.buildUniqueGlobalIndexer(Teuchos::opaqueWrapper(MPI_COMM_WORLD),physicsBlocks,conn_manager); + dofManager = dofManager_int; + + // construct some linear algebra object, build object to pass to evaluators + linObjFactory = Teuchos::rcp(new panzer::BlockedEpetraLinearObjFactory(comm.getConst(),dofManager_int)); + } + else { + const Teuchos::RCP > conn_manager + = Teuchos::rcp(new panzer_stk::STKConnManager(mesh)); + + panzer::DOFManagerFactory globalIndexerFactory; + RCP > dofManager_long + = globalIndexerFactory.buildUniqueGlobalIndexer(Teuchos::opaqueWrapper(MPI_COMM_WORLD),physicsBlocks,conn_manager); + dofManager = dofManager_long; + + // construct some linear algebra object, build object to pass to evaluators + linObjFactory = Teuchos::rcp(new panzer::TpetraLinearObjFactory(comm,dofManager_long)); + } + + // build worksets + //////////////////////////////////////////////////////// + + Teuchos::RCP wkstFactory + = Teuchos::rcp(new panzer_stk::WorksetFactory(mesh)); // build STK workset factory + Teuchos::RCP wkstContainer // attach it to a workset container (uses lazy evaluation) + = Teuchos::rcp(new panzer::WorksetContainer); + wkstContainer->setFactory(wkstFactory); + for(size_t i=0;isetNeeds(physicsBlocks[i]->elementBlockID(),physicsBlocks[i]->getWorksetNeeds()); + wkstContainer->setWorksetSize(workset_size); + wkstContainer->setGlobalIndexer(dofManager); + + // Setup STK response library for writing out the solution fields + //////////////////////////////////////////////////////////////////////// + Teuchos::RCP > stkIOResponseLibrary + = Teuchos::rcp(new panzer::ResponseLibrary(wkstContainer,dofManager,linObjFactory)); + + { + // get a vector of all the element blocks + std::vector eBlocks; + { + // get all element blocks and add them to the list + std::vector eBlockNames; + mesh->getElementBlockNames(eBlockNames); + for(std::size_t i=0;iaddResponse("Main Field Output",eBlocks,builder); + } + + // Setup response library for checking the error in this manufactered solution + //////////////////////////////////////////////////////////////////////// + + Teuchos::RCP > errorResponseLibrary + = Teuchos::rcp(new panzer::ResponseLibrary(wkstContainer,dofManager,linObjFactory)); + + { + std::vector eBlocks; + mesh->getElementBlockNames(eBlocks); + + panzer::FunctionalResponse_Builder builder; + + builder.comm = MPI_COMM_WORLD; + builder.cubatureDegree = 10; + builder.requiresCellIntegral = true; + builder.quadPointField = "EFIELD_ERROR"; + + errorResponseLibrary->addResponse("L2 Error",eBlocks,builder); + + builder.comm = MPI_COMM_WORLD; + builder.cubatureDegree = 10; + builder.requiresCellIntegral = true; + builder.quadPointField = "EFIELD_HCURL_ERROR"; + + errorResponseLibrary->addResponse("HCurl Error",eBlocks,builder); + } + + // setup closure model + ///////////////////////////////////////////////////////////// + + // Add in the application specific closure model factory + panzer::ClosureModelFactory_TemplateManager cm_factory; + Example::ClosureModelFactory_TemplateBuilder cm_builder; + cm_factory.buildObjects(cm_builder); + + Teuchos::ParameterList closure_models("Closure Models"); + { + closure_models.sublist("solid").sublist("SOURCE_EFIELD").set("Type","SIMPLE SOURCE"); // a constant source + // SOURCE_EFIELD field is required by the CurlLaplacianEquationSet + + // required for error calculation + closure_models.sublist("solid").sublist("EFIELD_ERROR").set("Type","L2 ERROR_CALC"); + closure_models.sublist("solid").sublist("EFIELD_ERROR").set("Field A","EFIELD"); + closure_models.sublist("solid").sublist("EFIELD_ERROR").set("Field B","EFIELD_EXACT"); + + closure_models.sublist("solid").sublist("EFIELD_HCURL_ERROR").set("Type","HCurl ERROR_CALC"); + closure_models.sublist("solid").sublist("EFIELD_HCURL_ERROR").set("Field A","EFIELD"); + closure_models.sublist("solid").sublist("EFIELD_HCURL_ERROR").set("Field B","EFIELD_EXACT"); + + closure_models.sublist("solid").sublist("EFIELD_EXACT").set("Type","EFIELD_EXACT"); + // closure_models.sublist("solid").sublist("EFIELD_EXACT").set("Type","SIMPLE SOURCE"); + } + + Teuchos::ParameterList user_data("User Data"); // user data can be empty here + + // setup field manager builder + ///////////////////////////////////////////////////////////// + + Teuchos::RCP fmb = + Teuchos::rcp(new panzer::FieldManagerBuilder); + fmb->setWorksetContainer(wkstContainer); + fmb->setupVolumeFieldManagers(physicsBlocks,cm_factory,closure_models,*linObjFactory,user_data); + fmb->setupBCFieldManagers(bcs,physicsBlocks,*eqset_factory,cm_factory,bc_factory,closure_models, + *linObjFactory,user_data); + + fmb->writeVolumeGraphvizDependencyFiles("volume",physicsBlocks); + fmb->writeBCGraphvizDependencyFiles("saucy"); + + // setup assembly engine + ///////////////////////////////////////////////////////////// + + // build assembly engine: The key piece that brings together everything and + // drives and controls the assembly process. Just add + // matrices and vectors + panzer::AssemblyEngine_TemplateManager ae_tm; + panzer::AssemblyEngine_TemplateBuilder builder(fmb,linObjFactory); + ae_tm.buildObjects(builder); + + // Finalize construcition of STK writer response library + ///////////////////////////////////////////////////////////// + { + user_data.set("Workset Size",workset_size); + stkIOResponseLibrary->buildResponseEvaluators(physicsBlocks, + cm_factory, + closure_models, + user_data,true,"io"); + + user_data.set("Workset Size",workset_size); + errorResponseLibrary->buildResponseEvaluators(physicsBlocks, + cm_factory, + closure_models, + user_data,true,"error"); + } + + // assemble linear system + ///////////////////////////////////////////////////////////// + + // build linear algebra objects: Ghost is for parallel assembly, it contains + // local element contributions summed, the global IDs + // are not unique. The non-ghosted or "global" + // container will contain the sum over all processors + // of the ghosted objects. The global indices are unique. + RCP ghostCont = linObjFactory->buildGhostedLinearObjContainer(); + RCP container = linObjFactory->buildLinearObjContainer(); + linObjFactory->initializeGhostedContainer(panzer::LinearObjContainer::X | + panzer::LinearObjContainer::F | + panzer::LinearObjContainer::Mat,*ghostCont); + linObjFactory->initializeContainer(panzer::LinearObjContainer::X | + panzer::LinearObjContainer::F | + panzer::LinearObjContainer::Mat,*container); + ghostCont->initialize(); + container->initialize(); + + // Actually evaluate + ///////////////////////////////////////////////////////////// + + panzer::AssemblyEngineInArgs input(ghostCont,container); + input.alpha = 0; + input.beta = 1; + + // evaluate physics: This does both the Jacobian and residual at once + ae_tm.getAsObject()->evaluate(input); + + // solve linear system + ///////////////////////////////////////////////////////////// + + if(useTpetra) + solveTpetraSystem(*container); + else + solveEpetraSystem(*container); + + // output data (optional) + ///////////////////////////////////////////////////////////// + + // write out solution + if(true) { + // fill STK mesh objects + Teuchos::RCP resp = stkIOResponseLibrary->getResponse("Main Field Output"); + panzer::AssemblyEngineInArgs respInput(ghostCont,container); + respInput.alpha = 0; + respInput.beta = 1; + + stkIOResponseLibrary->addResponsesToInArgs(respInput); + stkIOResponseLibrary->evaluate(respInput); + + // write to exodus + // --------------- + // Due to multiple instances of this test being run at the same + // time (one for each order), we need to differentiate output to + // prevent race conditions on output file. Multiple runs for the + // same order are ok as they are staged one after another in the + // ADD_ADVANCED_TEST cmake macro. + std::ostringstream filename; + filename << output_filename << basis_order << ".exo"; + mesh->writeToExodus(filename.str()); + } + + // compute error norm + ///////////////////////////////////////////////////////////// + + if(true) { + Teuchos::FancyOStream lout(Teuchos::rcpFromRef(std::cout)); + lout.setOutputToRootOnly(0); + + panzer::AssemblyEngineInArgs respInput(ghostCont,container); + respInput.alpha = 0; + respInput.beta = 1; + + Teuchos::RCP l2_resp = errorResponseLibrary->getResponse("L2 Error"); + Teuchos::RCP > l2_resp_func = + Teuchos::rcp_dynamic_cast >(l2_resp); + Teuchos::RCP > l2_respVec = Thyra::createMember(l2_resp_func->getVectorSpace()); + l2_resp_func->setVector(l2_respVec); + + Teuchos::RCP h1_resp = errorResponseLibrary->getResponse("HCurl Error"); + Teuchos::RCP > h1_resp_func = + Teuchos::rcp_dynamic_cast >(h1_resp); + Teuchos::RCP > h1_respVec = Thyra::createMember(h1_resp_func->getVectorSpace()); + h1_resp_func->setVector(h1_respVec); + + errorResponseLibrary->addResponsesToInArgs(respInput); + errorResponseLibrary->evaluate(respInput); + + lout << "L2 Error = " << sqrt(l2_resp_func->value) << std::endl; + lout << "HCurl Error = " << sqrt(h1_resp_func->value) << std::endl; + } + + stackedTimer->stop("Curl Laplacian"); + Teuchos::StackedTimer::OutputOptions options; + options.output_fraction = true; + options.output_minmax = true; + options.output_histogram = true; + options.num_histogram = 5; + stackedTimer->report(std::cout, Teuchos::DefaultComm::getComm(), options); + + // all done! + ///////////////////////////////////////////////////////////// + + if(useTpetra) + out << "ALL PASSED: Tpetra" << std::endl; + else + out << "ALL PASSED: Epetra" << std::endl; + + return 0; +} + +void solveEpetraSystem(panzer::LinearObjContainer & container) +{ + // convert generic linear object container to epetra container + panzer::EpetraLinearObjContainer & ep_container + = Teuchos::dyn_cast(container); + + ep_container.get_x()->PutScalar(0.0); + + // Setup the linear solve: notice A is used directly + Epetra_LinearProblem problem(&*ep_container.get_A(),&*ep_container.get_x(),&*ep_container.get_f()); + + // build the solver + AztecOO solver(problem); + solver.SetAztecOption(AZ_solver,AZ_gmres); // we don't push out dirichlet conditions + solver.SetAztecOption(AZ_precond,AZ_none); + solver.SetAztecOption(AZ_kspace,3000); + solver.SetAztecOption(AZ_output,1); + solver.SetAztecOption(AZ_precond,AZ_Jacobi); + + // solve the linear system + solver.Iterate(50000,1e-9); + + // we have now solved for the residual correction from + // zero in the context of a Newton solve. + // J*e = -r = -(f - J*0) where f = J*u + // Therefore we have J*e=-J*u which implies e = -u + // thus we will scale the solution vector + ep_container.get_x()->Scale(-1.0); +} + +void solveTpetraSystem(panzer::LinearObjContainer & container) +{ + typedef panzer::TpetraLinearObjContainer LOC; + + LOC & tp_container = Teuchos::dyn_cast(container); + + // do stuff + // Wrap the linear problem to solve in a Belos::LinearProblem + // object. The "X" argument of the LinearProblem constructor is + // only copied shallowly and will be overwritten by the solve, so we + // make a deep copy here. That way we can compare the result + // against the original X_guess. + typedef Tpetra::MultiVector MV; + typedef Tpetra::Operator OP; + typedef Belos::LinearProblem ProblemType; + Teuchos::RCP problem(new ProblemType(tp_container.get_A(), tp_container.get_x(), tp_container.get_f())); + TEUCHOS_ASSERT(problem->setProblem()); + + typedef Belos::PseudoBlockGmresSolMgr SolverType; + + Teuchos::ParameterList belosList; + belosList.set( "Num Blocks", 3000 ); // Maximum number of blocks in Krylov factorization + belosList.set( "Block Size", 1 ); // Blocksize to be used by iterative solver + belosList.set( "Maximum Iterations", 50000 ); // Maximum number of iterations allowed + belosList.set( "Maximum Restarts", 20 ); // Maximum number of restarts allowed + belosList.set( "Convergence Tolerance", 1e-9 ); // Relative convergence tolerance requested + belosList.set( "Verbosity", Belos::Errors + Belos::Warnings + Belos::TimingDetails + Belos::StatusTestDetails ); + belosList.set( "Output Frequency", 1 ); + belosList.set( "Output Style", 1 ); + + SolverType solver(problem, Teuchos::rcpFromRef(belosList)); + + Belos::ReturnType result = solver.solve(); + if (result == Belos::Converged) + std::cout << "Result: Converged." << std::endl; + else { + TEUCHOS_ASSERT(false); // FAILURE! + } + + // scale by -1 + tp_container.get_x()->scale(-1.0); + + tp_container.get_A()->resumeFill(); // where does this go? +} + +void testInitialization(const Teuchos::RCP& ipb, + std::vector& bcs, + const std::vector& eBlockNames, + int basis_order, + const bool threeD) +{ + { + std::string b_basis_type = threeD ? "HDiv" : "HVol"; + int b_basis_order = threeD ? basis_order : basis_order-1; + + Teuchos::ParameterList& p = ipb->sublist("CurlLapacian Physics"); + p.set("Type","CurlLaplacian"); + p.set("Model ID","solid"); + p.set("EField Basis Type","HCurl"); + p.set("BField Basis Type",b_basis_type); + p.set("EField Basis Order",basis_order); + p.set("BField Basis Order",b_basis_order); + p.set("Integration Order",10); + } + + std::size_t bc_id = 0; + for (const auto& block : eBlockNames) { + panzer::BCType bctype = panzer::BCT_Dirichlet; + std::string sideset_id = "left"; + std::string element_block_id = block; + std::string dof_name = "EFIELD"; + std::string strategy = "Constant"; + double value = 0.0; + Teuchos::ParameterList p; + p.set("Value",value); + panzer::BC bc(bc_id++, bctype, sideset_id, element_block_id, dof_name, + strategy, p); + bcs.push_back(bc); + } + + // multiblock splits on x only + for (const auto& block : eBlockNames) { + panzer::BCType bctype = panzer::BCT_Dirichlet; + std::string sideset_id = "top"; + std::string element_block_id = block; + std::string dof_name = "EFIELD"; + std::string strategy = "Constant"; + double value = 0.0; + Teuchos::ParameterList p; + p.set("Value",value); + panzer::BC bc(bc_id++, bctype, sideset_id, element_block_id, dof_name, + strategy, p); + bcs.push_back(bc); + } + + for (const auto& block : eBlockNames) { + panzer::BCType bctype = panzer::BCT_Dirichlet; + std::string sideset_id = "right"; + std::string element_block_id = block; + std::string dof_name = "EFIELD"; + std::string strategy = "Constant"; + double value = 0.0; + Teuchos::ParameterList p; + p.set("Value",value); + panzer::BC bc(bc_id++, bctype, sideset_id, element_block_id, dof_name, + strategy, p); + bcs.push_back(bc); + } + + // multiblock splits on x only + for (const auto& block : eBlockNames) { + panzer::BCType bctype = panzer::BCT_Dirichlet; + std::string sideset_id = "bottom"; + std::string element_block_id = block; + std::string dof_name = "EFIELD"; + std::string strategy = "Constant"; + double value = 0.0; + Teuchos::ParameterList p; + p.set("Value",value); + panzer::BC bc(bc_id++, bctype, sideset_id, element_block_id, dof_name, + strategy, p); + bcs.push_back(bc); + } + +} diff --git a/packages/panzer/disc-fe/src/Panzer_BasisValues2.hpp b/packages/panzer/disc-fe/src/Panzer_BasisValues2.hpp index 6664dce5367a..c2b57a247bdc 100644 --- a/packages/panzer/disc-fe/src/Panzer_BasisValues2.hpp +++ b/packages/panzer/disc-fe/src/Panzer_BasisValues2.hpp @@ -175,6 +175,11 @@ namespace panzer { const PHX::MDField & jac_inv, const PHX::MDField & weighted_measure); + void evaluateValues_HVol(const PHX::MDField & cub_points, + const PHX::MDField & jac_det, + const PHX::MDField & jac_inv, + const PHX::MDField & weighted_measure); + void evaluateValues_HGrad(const PHX::MDField & cub_points, const PHX::MDField & jac_inv, const PHX::MDField & weighted_measure); diff --git a/packages/panzer/disc-fe/src/Panzer_BasisValues2_impl.hpp b/packages/panzer/disc-fe/src/Panzer_BasisValues2_impl.hpp index f07ddfe6cfd5..a615dd831429 100644 --- a/packages/panzer/disc-fe/src/Panzer_BasisValues2_impl.hpp +++ b/packages/panzer/disc-fe/src/Panzer_BasisValues2_impl.hpp @@ -127,6 +127,20 @@ evaluateValues(const PHX::MDField & basis_vector.get_view()); } } + else if(elmtspace==PureBasis::HVOL) + { + Intrepid2::FunctionSpaceTools:: + HVOLtransformVALUE(basis_scalar.get_view(), + jac_det.get_view(), + basis_ref_scalar.get_view()); + + if(build_weighted) { + Intrepid2::FunctionSpaceTools:: + multiplyMeasure(weighted_basis_scalar.get_view(), + weighted_measure.get_view(), + basis_scalar.get_view()); + } + } else { TEUCHOS_ASSERT(false); } if(elmtspace==PureBasis::HGRAD && compute_derivatives) { @@ -260,6 +274,8 @@ evaluateValues(const PHX::MDField & if(elmtspace == PureBasis::CONST){ evaluateValues_Const(cub_points,jac_inv,weighted_measure); + } else if(elmtspace == PureBasis::HVOL){ + evaluateValues_HVol(cub_points,jac_det,jac_inv,weighted_measure); } else if(elmtspace == PureBasis::HGRAD){ evaluateValues_HGrad(cub_points,jac_inv,weighted_measure); } else if(elmtspace == PureBasis::HCURL){ @@ -353,6 +369,66 @@ evaluateValues_Const(const PHX::MDField +void panzer::BasisValues2:: +evaluateValues_HVol(const PHX::MDField & cub_points, + const PHX::MDField & jac_det, + const PHX::MDField & jac_inv, + const PHX::MDField & weighted_measure) +{ + + TEUCHOS_ASSERT(getElementSpace() == PureBasis::HVOL); + + typedef Intrepid2::FunctionSpaceTools fst; + MDFieldArrayFactory af("",ddims_,true); + + const panzer::PureBasis & basis = *(basis_layout->getBasis()); + + const int num_points = basis_layout->numPoints(); + const int num_basis = basis.cardinality(); + const int num_dim = basis_layout->dimension(); + const int num_cells = basis_layout->numCells(); + + auto cell_basis_scalar = af.buildStaticArray("cell_basis_scalar",1,num_basis,num_points); + auto cell_cub_points = af.buildStaticArray("cell_cub_points",num_points,num_dim); + auto cell_jac_det = af.buildStaticArray("cell_jac_det",1,num_points); + + auto cell_basis_ref_scalar = af.buildStaticArray("cell_basis_ref_scalar",num_basis,num_points); + + for(int cell=0;cellgetValues(cell_basis_ref_scalar.get_view(),cell_cub_points.get_view(),Intrepid2::OPERATOR_VALUE); + + // ============================================= + // Transform reference values to physical values + + fst::HVOLtransformVALUE(cell_basis_scalar.get_view(),cell_jac_det.get_view(),cell_basis_ref_scalar.get_view()); + for(int b=0;b @@ -701,6 +777,31 @@ evaluateValuesCV(const PHX::MDField for (int ip = 0; ip < num_ip; ++ip) basis_scalar(icell,b,ip) = dyn_basis_ref_scalar(b,ip); + } + if(elmtspace==PureBasis::HVOL) { + ArrayDynamic dyn_basis_ref_scalar = af.buildArray("dyn_basis_ref_scalar",num_card,num_ip); + + intrepid_basis->getValues(dyn_basis_ref_scalar.get_view(), + dyn_cub_points.get_view(), + Intrepid2::OPERATOR_VALUE); + + int one_cell= 1; + ArrayDynamic dyn_basis_scalar = af.buildArray("dyn_basis_vector",one_cell,num_card,num_ip); + ArrayDynamic dyn_jac = af.buildArray("dyn_jac",one_cell,num_ip,num_dim,num_dim); + ArrayDynamic dyn_jac_det = af.buildArray("dyn_jac_det",one_cell,num_ip); + + int cellInd = 0; + for (int ip = 0; ip < num_ip; ++ip) + dyn_jac_det(cellInd,ip) = jac_det(icell,ip); + + Intrepid2::FunctionSpaceTools::HVOLtransformVALUE(dyn_basis_scalar.get_view(), + dyn_jac_det.get_view(), + dyn_basis_ref_scalar.get_view()); + + for (int b = 0; b < num_card; ++b) + for (int ip = 0; ip < num_ip; ++ip) + basis_scalar(icell,b,ip) = dyn_basis_scalar(0,b,ip); + } if(elmtspace==PureBasis::HGRAD) { ArrayDynamic dyn_basis_ref_scalar = af.buildArray("dyn_basis_ref_scalar",num_card,num_ip); @@ -901,7 +1002,7 @@ evaluateReferenceValues(const PHX::MDField & cub_points,bool comp dyn_cub_points(ip,d) = cub_points(ip,d); PureBasis::EElementSpace elmtspace = getElementSpace(); - if(elmtspace==PureBasis::HGRAD || elmtspace==PureBasis::CONST) { + if(elmtspace==PureBasis::HGRAD || elmtspace==PureBasis::CONST || elmtspace==PureBasis::HVOL) { ArrayDynamic dyn_basis_ref_scalar = af.buildArray("dyn_basis_ref_scalar",num_card,num_quad); intrepid_basis->getValues(dyn_basis_ref_scalar.get_view(), @@ -1513,7 +1614,7 @@ setupArrays(const Teuchos::RCP& layout, weighted_div_basis = af.buildStaticArray("weighted_div_basis",numcells,card,num_quad); } } - else if(elmtspace==panzer::PureBasis::CONST) { + else if(elmtspace==panzer::PureBasis::CONST || elmtspace==panzer::PureBasis::HVOL) { // CONST is a nodal field // build values diff --git a/packages/panzer/disc-fe/src/Panzer_IntrepidBasisFactory.hpp b/packages/panzer/disc-fe/src/Panzer_IntrepidBasisFactory.hpp index ae78635273fa..6f0a3f48549c 100644 --- a/packages/panzer/disc-fe/src/Panzer_IntrepidBasisFactory.hpp +++ b/packages/panzer/disc-fe/src/Panzer_IntrepidBasisFactory.hpp @@ -48,6 +48,8 @@ #include #include "Intrepid2_HVOL_C0_FEM.hpp" +#include "Intrepid2_HVOL_TRI_Cn_FEM.hpp" +#include "Intrepid2_HVOL_QUAD_Cn_FEM.hpp" #include "Teuchos_RCP.hpp" #include "Intrepid2_Basis.hpp" @@ -137,6 +139,15 @@ namespace panzer { if ( (basis_type == "Const") && (basis_order == 0) ) basis = Teuchos::rcp( new Intrepid2::Basis_HVOL_C0_FEM(cell_topology) ); + else if ( (basis_type == "HVol") && (basis_order == 0) ) + basis = Teuchos::rcp( new Intrepid2::Basis_HVOL_C0_FEM(cell_topology) ); + + else if ( (basis_type == "HVol") && (cell_type == "Quadrilateral") && (basis_order > 0) ) + basis = Teuchos::rcp( new Intrepid2::Basis_HVOL_QUAD_Cn_FEM(basis_order, point_type) ); + + else if ( (basis_type == "HVol") && (cell_type == "Triangle") && (basis_order > 0) ) + basis = Teuchos::rcp( new Intrepid2::Basis_HVOL_TRI_Cn_FEM(basis_order, point_type) ); + else if ( (basis_type == "HGrad") && (cell_type == "Hexahedron") && (basis_order == 1) ) basis = Teuchos::rcp( new Intrepid2::Basis_HGRAD_HEX_C1_FEM ); diff --git a/packages/panzer/disc-fe/src/Panzer_L2Projection_impl.hpp b/packages/panzer/disc-fe/src/Panzer_L2Projection_impl.hpp index 250889ceb4f7..a313e057ae27 100644 --- a/packages/panzer/disc-fe/src/Panzer_L2Projection_impl.hpp +++ b/packages/panzer/disc-fe/src/Panzer_L2Projection_impl.hpp @@ -94,7 +94,7 @@ namespace panzer { auto M = ghostedMatrix->getLocalMatrix(); const int fieldIndex = targetGlobalIndexer_->getFieldNum(targetBasisDescriptor_.getType()); - const bool is_scalar = targetBasisDescriptor_.getType()=="HGrad" || targetBasisDescriptor_.getType()=="Const"; + const bool is_scalar = targetBasisDescriptor_.getType()=="HGrad" || targetBasisDescriptor_.getType()=="Const" || targetBasisDescriptor_.getType()=="HVol"; // Loop over element blocks and fill mass matrix if(is_scalar){ @@ -348,7 +348,7 @@ namespace panzer { Kokkos::View sourceUnweightedScalarBasis; Kokkos::View sourceUnweightedVectorBasis; bool useRankThreeBasis = false; // default to gradient or vector basis - if ( (sourceBasisDescriptor.getType() == "HGrad") || (sourceBasisDescriptor.getType() == "Const") ) { + if ( (sourceBasisDescriptor.getType() == "HGrad") || (sourceBasisDescriptor.getType() == "Const") || (sourceBasisDescriptor.getType() == "HVol") ) { if (directionIndex == -1) { // Project dof value sourceUnweightedScalarBasis = sourceBasisValues.basis_scalar.get_static_view(); useRankThreeBasis = true; diff --git a/packages/panzer/disc-fe/src/Panzer_PureBasis.cpp b/packages/panzer/disc-fe/src/Panzer_PureBasis.cpp index bf67823cbd9f..0aee31cae0bc 100644 --- a/packages/panzer/disc-fe/src/Panzer_PureBasis.cpp +++ b/packages/panzer/disc-fe/src/Panzer_PureBasis.cpp @@ -122,6 +122,8 @@ void panzer::PureBasis::initialize(const std::string & in_basis_type,const int i element_space_ = HDIV; else if(basis_type_=="Const") element_space_ = CONST; + else if(basis_type_=="HVol") + element_space_ = HVOL; else { TEUCHOS_TEST_FOR_EXCEPTION(true,std::invalid_argument, "PureBasis::initializeIntrospection - Invalid basis name \"" << basis_type_ << "\""); } @@ -130,6 +132,9 @@ void panzer::PureBasis::initialize(const std::string & in_basis_type,const int i case CONST: basis_rank_ = 0; break; + case HVOL: + basis_rank_ = 0; + break; case HGRAD: basis_rank_ = 0; break; diff --git a/packages/panzer/disc-fe/src/Panzer_PureBasis.hpp b/packages/panzer/disc-fe/src/Panzer_PureBasis.hpp index 1dbc3f1faebf..a7cd7229981e 100644 --- a/packages/panzer/disc-fe/src/Panzer_PureBasis.hpp +++ b/packages/panzer/disc-fe/src/Panzer_PureBasis.hpp @@ -61,7 +61,7 @@ namespace panzer { class PureBasis { public: - typedef enum { HGRAD=0, HCURL=1, HDIV=2, CONST=3 } EElementSpace; + typedef enum { HGRAD=0, HCURL=1, HDIV=2, HVOL=3, CONST=4 } EElementSpace; /** Build a basis given a type, order and CellData object \param[in] basis_type String name that describes the type of basis @@ -139,7 +139,7 @@ namespace panzer { { return getElementSpace()==HCURL || getElementSpace()==HDIV; } bool isScalarBasis() const - { return getElementSpace()==HGRAD || getElementSpace()==CONST; } + { return getElementSpace()==HGRAD || getElementSpace()==CONST || getElementSpace()==HVOL; } int getBasisRank() const { return basis_rank_; }