diff --git a/.codedocs b/.codedocs index 8cb73be71..d8d13cd4c 100644 --- a/.codedocs +++ b/.codedocs @@ -1,5 +1,4 @@ # CodeDocs.xyz Configuration File -EXCLUDE_PATTERNS = */deprecated/* -DOXYFILE = doc/doxygen.cfg.in +DOXYFILE = doc/doxygen_codedocs.cfg diff --git a/.travis.yml b/.travis.yml index b4788012b..9058a4787 100644 --- a/.travis.yml +++ b/.travis.yml @@ -1 +1,36 @@ +os: linux +dist: bionic language: c++ + +install: +- chmod u+x ci/travis-ci.sh + +stages: +- doc +#- build and test +#- deploy + +jobs: + include: +# - stage: build and test +# env: TITLE="Build CCORE (Linux)." +# script: ./ci/travis-ci.sh BUILD_CCORE +# - state: build and test +# env: TITLE="Run C/C++ static code analyser." +# script: ./ci/travis-ci.sh ANALYSE_CCORE +# - stage: build and test +# env: TITLE="Run unit-tests for CCORE." +# script: ./ci/travis-ci.sh UT_CCORE +# - stage: build and test +# env: TITLE="Run memory leak analyser (valgrind) for CCORE." +# script: ./ci/travis-ci.sh VALGRIND_CCORE +# - stage: build and test +# env: TITLE="PyClustering unit and integration testing." +# script: ./ci/travis-ci.sh TEST_PYCLUSTERING +# - stage: build and test +# os: osx +# env: TITLE="Build CCORE (MacOS) and intergration testing." +# script: ./ci/travis-ci.sh BUILD_TEST_CCORE_MACOS + - stage: doc + env: TITLE="Build documentation." + script: ./ci/travis-ci.sh DOCUMENTATION diff --git a/CMakeLists.txt b/CMakeLists.txt index c05b558aa..908b083bf 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -15,45 +15,48 @@ set(MPIMAX 4 CACHE STRING "Default number of processors used in ctest mpirun -np set(EXECUTABLE_OUTPUT_PATH ${PROJECT_BINARY_DIR}/bin) set(CLEAN_UP_FILES ./bin/* ./CMakeCache.txt) -# Find deal.ii library -find_package(deal.II 9.0.1 QUIET - HINTS ${DEAL_II_DIR} ../ ../../ $ENV{DEAL_II_DIR} - ) -IF(NOT ${deal.II_FOUND}) - message(FATAL_ERROR "\n" - "*** Could not locate a (sufficiently recent) version of deal.II. ***\n\n" - "You may want to either pass a flag -DDEAL_II_DIR=/path/to/deal.II to cmake\n" - "or set an environment variable \"DEAL_II_DIR\" that contains this path." - ) -ENDIF() -DEAL_II_INITIALIZE_CACHED_VARIABLES() - -# Make sure deal.II has the proper external dependencies -IF(NOT ${DEAL_II_WITH_MPI}) - message(FATAL_ERROR "\n" "*** deal.II needs to be configured with -DDEAL_II_WITH_MPI=ON ***\n\n") -ENDIF() -IF(NOT ${DEAL_II_WITH_TRILINOS}) - message(FATAL_ERROR "\n" "*** deal.II needs to be configured with -DDEAL_II_WITH_TRILINOS=ON ***\n\n") -ENDIF() -IF(NOT ${DEAL_II_WITH_P4EST}) - message(FATAL_ERROR "\n" "*** deal.II needs to be configured with -DDEAL_II_WITH_P4EST=ON ***\n\n") -ENDIF() -IF(NOT ${DEAL_II_WITH_METIS}) - message(FATAL_ERROR "\n" "*** deal.II needs to be configured with -DDEAL_II_WITH_METIS=ON ***\n\n") -ENDIF() - -# Use ld.gold for faster linking -option(USE_LD_GOLD "Use GNU gold linker" OFF) -if(USE_LD_GOLD AND "${CMAKE_C_COMPILER_ID}" STREQUAL "GNU") - execute_process(COMMAND ${CMAKE_C_COMPILER} -fuse-ld=gold -Wl,--version OUTPUT_VARIABLE stdout ERROR_QUIET) - if("${stdout}" MATCHES "GNU gold") - set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -fuse-ld=gold") - set(CMAKE_C_FLAGS "${CMAKE_C_FLAGS} -fuse-ld=gold") - else() - message(WARNING "GNU gold linker isn't available, using the default system linker.") - endif() -endif() - +option(DOC_ONLY "Only generate code for documentation sake. Not compilable." OFF) +if(NOT DOC_ONLY) + # Find deal.ii library + find_package(deal.II 9.0.1 QUIET + HINTS ${DEAL_II_DIR} ../ ../../ $ENV{DEAL_II_DIR} + ) + IF(NOT ${deal.II_FOUND}) + message(FATAL_ERROR "\n" + "*** Could not locate a (sufficiently recent) version of deal.II. ***\n\n" + "You may want to either pass a flag -DDEAL_II_DIR=/path/to/deal.II to cmake\n" + "or set an environment variable \"DEAL_II_DIR\" that contains this path." + ) + ENDIF() + DEAL_II_INITIALIZE_CACHED_VARIABLES() + + # Make sure deal.II has the proper external dependencies + IF(NOT ${DEAL_II_WITH_MPI}) + message(FATAL_ERROR "\n" "*** deal.II needs to be configured with -DDEAL_II_WITH_MPI=ON ***\n\n") + ENDIF() + IF(NOT ${DEAL_II_WITH_TRILINOS}) + message(FATAL_ERROR "\n" "*** deal.II needs to be configured with -DDEAL_II_WITH_TRILINOS=ON ***\n\n") + ENDIF() + IF(NOT ${DEAL_II_WITH_P4EST}) + message(FATAL_ERROR "\n" "*** deal.II needs to be configured with -DDEAL_II_WITH_P4EST=ON ***\n\n") + ENDIF() + IF(NOT ${DEAL_II_WITH_METIS}) + message(FATAL_ERROR "\n" "*** deal.II needs to be configured with -DDEAL_II_WITH_METIS=ON ***\n\n") + ENDIF() + + # Use ld.gold for faster linking + option(USE_LD_GOLD "Use GNU gold linker" OFF) + if(USE_LD_GOLD AND "${CMAKE_C_COMPILER_ID}" STREQUAL "GNU") + execute_process(COMMAND ${CMAKE_C_COMPILER} -fuse-ld=gold -Wl,--version OUTPUT_VARIABLE stdout ERROR_QUIET) + if("${stdout}" MATCHES "GNU gold") + set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -fuse-ld=gold") + set(CMAKE_C_FLAGS "${CMAKE_C_FLAGS} -fuse-ld=gold") + else() + message(WARNING "GNU gold linker isn't available, using the default system linker.") + endif() + endif() +endif() # DOC_ONLY + add_custom_target(1D COMMAND ${CMAKE_COMMAND} -E echo "Makes all 1D executables, including tests. Allows ctest -R 1D") add_custom_target(2D COMMAND ${CMAKE_COMMAND} -E echo "Makes all 2D executables, including tests. Allows ctest -R 2D") add_custom_target(3D COMMAND ${CMAKE_COMMAND} -E echo "Makes all 3D executables, including tests. Allows ctest -R 3D") diff --git a/INSTALL.md b/INSTALL.md index 428ad264b..879fdf7cf 100644 --- a/INSTALL.md +++ b/INSTALL.md @@ -1,3 +1,64 @@ +## deal.II + +This is the main library being used by this code. Since we are using advanced features, we only support the `master` of deal.II, which means it will have to be installed from source by cloning their [repository](https://github.com/dealii/dealii). + +Most of the packages are readily available through `apt install`. p4est might need to be installed from source since the current p4est version available on `apt` is lower than what is required by deal.II (p4est 2.0). A small set of instructions is available [here](https://www.dealii.org/current/external-libs/p4est.html). + +There is an [example script](doc/install_dealii.sh) for what has been used to install deal.II. You may need to provide the path such as `-DTRILINOS_DIR` if CMake does not find the package on its own. + +The deal.II library has been setup with the following options: + +~~~~ + deal.II configuration: + CMAKE_BUILD_TYPE: DebugRelease + BUILD_SHARED_LIBS: ON + CMAKE_INSTALL_PREFIX: /home/ddong/Codes/dealii/install + CMAKE_SOURCE_DIR: /home/ddong/Codes/dealii + (version 9.0.1, shortrev 097bf59e49) + CMAKE_BINARY_DIR: /home/ddong/Codes/dealii/build + CMAKE_CXX_COMPILER: GNU 7.3.0 on platform Linux x86_64 + /usr/bin/mpicxx + + Configured Features (DEAL_II_ALLOW_BUNDLED = ON, DEAL_II_ALLOW_AUTODETECTION = ON): + ( DEAL_II_WITH_64BIT_INDICES = OFF ) + DEAL_II_WITH_ADOLC set up with external dependencies + DEAL_II_WITH_ARPACK set up with external dependencies + ( DEAL_II_WITH_ASSIMP = OFF ) + DEAL_II_WITH_BOOST set up with external dependencies + ( DEAL_II_WITH_CUDA = OFF ) + DEAL_II_WITH_CXX14 = ON + DEAL_II_WITH_CXX17 = ON + DEAL_II_WITH_GMSH set up with external dependencies + DEAL_II_WITH_GSL set up with external dependencies + ( DEAL_II_WITH_HDF5 = OFF ) + DEAL_II_WITH_LAPACK set up with external dependencies + DEAL_II_WITH_METIS set up with external dependencies + DEAL_II_WITH_MPI set up with external dependencies + DEAL_II_WITH_MUPARSER set up with bundled packages + ( DEAL_II_WITH_NANOFLANN = OFF ) + ( DEAL_II_WITH_NETCDF = OFF ) + DEAL_II_WITH_OPENCASCADE set up with external dependencies + DEAL_II_WITH_P4EST set up with external dependencies + DEAL_II_WITH_PETSC set up with external dependencies + DEAL_II_WITH_SCALAPACK set up with external dependencies + DEAL_II_WITH_SLEPC set up with external dependencies + DEAL_II_WITH_SUNDIALS set up with external dependencies + ( DEAL_II_WITH_THREADS = OFF ) + DEAL_II_WITH_TRILINOS set up with external dependencies + DEAL_II_WITH_UMFPACK set up with external dependencies + DEAL_II_WITH_ZLIB set up with external dependencies + + Component configuration: + ( DEAL_II_COMPONENT_DOCUMENTATION = OFF ) + DEAL_II_COMPONENT_EXAMPLES + ( DEAL_II_COMPONENT_PACKAGE = OFF ) + ( DEAL_II_COMPONENT_PYTHON_BINDINGS = OFF ) + + Detailed information (compiler flags, feature configuration) can be found in detailed.log + + Run $ make info to print a help message with a list of top level targets +~~~~ + ## Installation of PHiLiP on Beluga cluster This section is aimed McGill's group who use Compute Canada's Beluga cluster. @@ -72,65 +133,6 @@ since it takes a very long time. If you have any questions, feel free to contact me. -## deal.II - -This is the main library being used by this code. Most of the packages are readily available through apt. p4est might need to be installed from source since the apt version is lower than what is required by deal.II. - -There is an [example script](install_dealii.sh) for what has been used to install deal.II. You may need to provide the path such as `-DTRILINOS_DIR` if CMake does not find the package on its own. - -The deal.II library has been setup with the following options: - -~~~~ - deal.II configuration: - CMAKE_BUILD_TYPE: DebugRelease - BUILD_SHARED_LIBS: ON - CMAKE_INSTALL_PREFIX: /home/ddong/Codes/dealii/install - CMAKE_SOURCE_DIR: /home/ddong/Codes/dealii - (version 9.0.1, shortrev 097bf59e49) - CMAKE_BINARY_DIR: /home/ddong/Codes/dealii/build - CMAKE_CXX_COMPILER: GNU 7.3.0 on platform Linux x86_64 - /usr/bin/mpicxx - - Configured Features (DEAL_II_ALLOW_BUNDLED = ON, DEAL_II_ALLOW_AUTODETECTION = ON): - ( DEAL_II_WITH_64BIT_INDICES = OFF ) - DEAL_II_WITH_ADOLC set up with external dependencies - DEAL_II_WITH_ARPACK set up with external dependencies - ( DEAL_II_WITH_ASSIMP = OFF ) - DEAL_II_WITH_BOOST set up with external dependencies - ( DEAL_II_WITH_CUDA = OFF ) - DEAL_II_WITH_CXX14 = ON - DEAL_II_WITH_CXX17 = ON - DEAL_II_WITH_GMSH set up with external dependencies - DEAL_II_WITH_GSL set up with external dependencies - ( DEAL_II_WITH_HDF5 = OFF ) - DEAL_II_WITH_LAPACK set up with external dependencies - DEAL_II_WITH_METIS set up with external dependencies - DEAL_II_WITH_MPI set up with external dependencies - DEAL_II_WITH_MUPARSER set up with bundled packages - ( DEAL_II_WITH_NANOFLANN = OFF ) - ( DEAL_II_WITH_NETCDF = OFF ) - DEAL_II_WITH_OPENCASCADE set up with external dependencies - DEAL_II_WITH_P4EST set up with external dependencies - DEAL_II_WITH_PETSC set up with external dependencies - DEAL_II_WITH_SCALAPACK set up with external dependencies - DEAL_II_WITH_SLEPC set up with external dependencies - DEAL_II_WITH_SUNDIALS set up with external dependencies - ( DEAL_II_WITH_THREADS = OFF ) - DEAL_II_WITH_TRILINOS set up with external dependencies - DEAL_II_WITH_UMFPACK set up with external dependencies - DEAL_II_WITH_ZLIB set up with external dependencies - - Component configuration: - ( DEAL_II_COMPONENT_DOCUMENTATION = OFF ) - DEAL_II_COMPONENT_EXAMPLES - ( DEAL_II_COMPONENT_PACKAGE = OFF ) - ( DEAL_II_COMPONENT_PYTHON_BINDINGS = OFF ) - - Detailed information (compiler flags, feature configuration) can be found in detailed.log - - Run $ make info to print a help message with a list of top level targets -~~~~ - To compile deal.II on the Beluga cluster, load the following modules. ~~~~ diff --git a/README.md b/README.md index 6c5895dc0..e4dc0e69f 100644 --- a/README.md +++ b/README.md @@ -16,6 +16,16 @@ The math supporting this code can be viewed in this **very rough draft in progre - Supported elements: LINEs, QUADs, HEXs since it uses deal.II - Supported refinements: h (size) or p (order). +## Documentation + +The code itself is documented using Doxygen, and the latest documentation is hosted on [codedocs.xyz](https://codedocs.xyz/dougshidong/PHiLiP/). + +Since deal.II is heavily used, their [documentation](https://www.dealii.org/developer/doxygen/deal.II/index.html) is probably the most useful. + +Another great ressource is the [deal.II Google Groups](https://groups.google.com/forum/#!forum/dealii), where developers are actively answering questions. + +Finally, I am also always available to answer questions regarding the code by e-mail at doug.shi-dong@mail.mcgill.ca + ## Building/Running the Code The code has been succesfully built in the following environments: @@ -55,10 +65,12 @@ The html documentation can be accessed by pointing a browser at `ROOT/doc/html/i ## Testing +A list of currently known failing tests is kept in the [GitHub issues](https://github.com/dougshidong/PHiLiP/issues?q=is%3Aissue+is%3Aopen+label%3Atestfail) with `testfail` tags. + Testing can be performed using CMake's `ctest` functionality. After successfully compiling the project, all tests can be run by executing: ```sh -$ ROOT$ make test (which is equivalent to ROOT$ make test) +$ ROOT$ ctest (which is equivalent to ROOT$ make test) ``` An alternative make target is provided to run tests with --output-on-failure: @@ -87,7 +99,10 @@ For a serial run, you may simply use gdb as intended ROOT$ gdb --args commmand_to_launch_test GDB$ run (Executes the program. Can re-launch the program if you forgot to put breakpoints.) ``` -For example `--args /home/ddong/Codes/PHiLiP_temp/PHiLiP/build_debug/bin/PHiLiP_2D "-i" "/home/ddong/Codes/PHiLiP_temp/PHiLiP/build_de bug/tests/adv ection_implicit/2d_advection_implicit_strong.prm`. +For example +``` +gdb --args /home/ddong/Codes/PHiLiP_temp/PHiLiP/build_debug/bin/PHiLiP_2D "-i" "/home/ddong/Codes/PHiLiP_temp/PHiLiP/build_debug/tests/advection_implicit/2d_advection_implicit_strong.prm +``` Additional useful commands are: @@ -99,6 +114,14 @@ GDB$ next (Execute the next line of code in the function. Will NOT go into the f GDB$ quit ``` +### Memory + +Memory leaks can be detected using Valgrind's tool `memcheck`. The application must be compiled in `Debug` mode. For example + +``` +valgrind --leak-check=full --track-origins=yes /home/ddong/Codes/PHiLiP/build_debug/bin/2D_HighOrder_MappingFEField +``` + ### Parallel debugging If the error only occurs when using parallelism, you can use the following example command @@ -107,6 +130,45 @@ mpirun -np 2 xterm -hold -e gdb -ex 'break MPI_Abort' -ex run --args /home/ddong ``` This launches 2 xterm processes, each of which will launch gdb processes that will run the code and will have a breakpoint when MPI_Abort is encountered. +## Performance + +Problems tend to show up in the 3D version if an algorithm has been implemented inefficiently. It is therefore highly recommended that a 3D test accompanies the implemented features. + +### Computational + +Computational bottlenecks can be inspected using Valgrind's tool `callgrind`. It is used as such + +``` +valgrind --tool=callgrind /home/ddong/Codes/PHiLiP/build_release/bin/2D_RBF_mesh_movement +``` + +This will result in a `callgrind.out.#####`. A visualizer such as `kcachegrind` (available through `apt`) can then be used to sort through the results. For example: + +```kcachegrind callgrind.out.24250``` + +### Memory + +Apart from memory leaks, it is possible that some required allocations demand too much memory. Valgrind also has a tool for this called `massif`. For example + +```valgrind --tool=massif /home/ddong/Codes/PHiLiP/build_debug/bin/3D_RBF_mesh_movement``` + +will generate a `massif.out.#####` file that can be visualized using `massif-visualizer` (available through `apt`) as + +```massif-visualizer massif.out.18580``` + + +## Contributing checklist + +1. A unit test, integration test, or regression test accompanies the feature. Tests longer than a few seconds should be tagged as with the suffix MEDIUM, and tests a minute or longer should be tagged with LONG. + * A unit test is often most appropriate, and is aimed at testing a single component of the code. See the test on [Euler's primitive to conservative conversion](https://github.com/dougshidong/PHiLiP/blob/master/tests/unit_tests/euler_unit_test/euler_convert_primitive_conservative.cpp) + * An integration test runs the entire main program by taking an input file and calling PHiLiP_1/2/3D. It should be derived from the [`TestBase` class](https://github.com/dougshidong/PHiLiP/blob/master/src/testing/tests.h), and have a control file located in the [integration test directory](https://github.com/dougshidong/PHiLiP/tree/master/tests/integration_tests_control_files). Since integrations tests uses multiple components, they usually take longer. Furthermore, the cause of failure is sometimes less obvious. A good suggestion is to use an existing test control file, and only change 1 parameter to help pinpoint issues when it fails. + * A regression test stores previously computed data to validate future results. Note that this type of test is rarely appropriate since valid changes in the code can fail this type of test. If implemented, a script/code should be made available such that newly computed results can replace the old results. See [file1](https://github.com/dougshidong/PHiLiP/blob/master/tests/unit_tests/regression/jacobian_matrix_regression.cpp) and [file2](https://github.com/dougshidong/PHiLiP/blob/master/tests/unit_tests/regression/matrix_data/copy_matrices.sh) +2. The feature has been documented. + * Function and member variable documentation should be presented in the associated header file. `make doc` should generate a html file in the `/path_to_build/doc/html/index.html` that can be opened used your browser of choice. + * Comments in the code as appropriate. +3. The `master` branch of `https://github.com/dougshidong/PHiLiP` has been merged into your fork and merge conflicts have been resolved. +4. The entire `ctest` suite has been run in `Release` mode and the short/medium length tests have been run in `Debug` mode (using `ctest -E LONG`). Please save the log for the next point. +5. Submit a pull request with a log of the tests. # License diff --git a/ci/travis-ci.sh b/ci/travis-ci.sh new file mode 100755 index 000000000..851185218 --- /dev/null +++ b/ci/travis-ci.sh @@ -0,0 +1,245 @@ +# +# @authors Andrei Novikov (pyclustering@yandex.ru) +# @date 2014-2019 +# @copyright GNU Public License +# +# GNU_PUBLIC_LICENSE +# pyclustering is free software: you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# pyclustering is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with this program. If not, see . +# + + +# CCORE_X64_BINARY_FOLDER=pyclustering/core/x64/linux +# CCORE_X64_BINARY_PATH=$CCORE_X64_BINARY_FOLDER/ccore.so +# +# CCORE_X86_BINARY_FOLDER=pyclustering/core/x86/linux +# CCORE_X86_BINARY_PATH=$CCORE_X86_BINARY_FOLDER/ccore.so +# +DOXYGEN_FILTER=( "warning: Unexpected new line character" ) + + +print_error() { + echo "[PHiLiP CI] ERROR: $1" +} + + +print_info() { + echo "[PHiLiP CI] INFO: $1" +} + + +check_failure() { + if [ $? -ne 0 ] ; then + if [ -z $1 ] ; then + print_error $1 + else + print_error "Failure exit code is detected." + fi + exit 1 + fi +} + + +check_error_log_file() { + problems_amount=$(cat $1 | wc -l) + printf "Total amount of errors and warnings: '%d'\n" "$problems_amount" + + if [ $problems_amount -ne 0 ] ; then + print_info "List of warnings and errors:" + cat $1 + + print_error $2 + exit 1 + fi +} + + +build_ccore() { + cd $TRAVIS_BUILD_DIR/ccore/ + + [ -f stderr.log ] && rm stderr.log + [ -f stdout.log ] && rm stdout.log + + if [ "$1" == "x64" ]; then + make ccore_x64 > >(tee -a stdout.log) 2> >(tee -a stderr.log >&2) + check_error_log_file stderr.log "Building CCORE (x64): FAILURE." + elif [ "$1" == "x86" ]; then + make ccore_x86 > >(tee -a stdout.log) 2> >(tee -a stderr.log >&2) + check_error_log_file stderr.log "Building CCORE (x86): FAILURE." + else + print_error "Unknown CCORE platform is specified." + exit 1 + fi + + cd - +} + + +run_build_ccore_job() { + print_info "CCORE (C++ code building):" + print_info "- Build CCORE library for x64 platform." + print_info "- Build CCORE library for x86 platform." + + #install requirement for the job + print_info "Install requirement for CCORE building." + + sudo apt-get install -qq g++-5 + sudo apt-get install -qq g++-5-multilib + sudo update-alternatives --install /usr/bin/g++ g++ /usr/bin/g++-5 50 + + # show info + g++ --version + gcc --version + + # build ccore library + build_ccore x64 + build_ccore x86 + + print_info "Upload ccore x64 binary." + upload_binary x64 linux + + print_info "Upload ccore x86 binary." + upload_binary x86 linux +} + +run_valgrind_ccore_job() { + print_info "VALGRIND CCORE (C++ code valgrind shock checking):" + print_info "- Run unit-tests of pyclustering." + print_info "- Shock memory leakage detection by valgrind." + + # install requirements for the job + sudo apt-get install -qq g++-5 + sudo apt-get install -qq g++-5-multilib + sudo apt-get install -qq valgrind + sudo update-alternatives --install /usr/bin/g++ g++ /usr/bin/g++-5 50 + + # build and run unit-test project under valgrind to check memory leakage + cd ccore/ + + make valgrind_shock + check_failure "CCORE shock memory leakage status: FAILURE." +} + +run_integration_test_job() { + print_info "INTEGRATION TESTING ('ccore' <-> 'pyclustering' for platform '$1')." + print_info "- Build CCORE library." + print_info "- Run integration tests of pyclustering." + + PLATFORM_TARGET=$1 + + # install requirements for the job + install_miniconda $PLATFORM_TARGET + + sudo apt-get install -qq g++-5 gcc-5 + sudo apt-get install -qq g++-5-multilib gcc-5-multilib + sudo update-alternatives --install /usr/bin/g++ g++ /usr/bin/g++-5 50 + sudo update-alternatives --install /usr/bin/gcov gcov /usr/bin/gcov-5 50 + + # build ccore library + build_ccore $PLATFORM_TARGET + + # run integration tests + python pyclustering/tests/tests_runner.py --integration +} + + +run_doxygen_job() { + print_info "DOXYGEN (documentation generation)." + print_info "- Generate documentation and check for warnings." + + sudo apt update + print_info "Install cmake" + sudo apt install cmake + + print_info "Install doxygen" + sudo apt install doxygen + + print_info "Install requirements for doxygen." + sudo apt install graphviz + sudo apt install texlive + + rm -rf build_doc; + mkdir build_doc; cd build_doc; + cmake ../ -DDOC_ONLY=ON + + print_info "Prepare log files." + report_file=doxygen_problems.log + report_file_filtered=doxygen_problems_filtered.log + + rm -f $report_file + rm -f $report_file_filtered + + print_info "Generate documentation." + doxygen --version +#make doc 2>&1 | tee $report_file > >(while read line; do echo -e "\e[01;31m$line\e[0m" >&2; done) + make doc > $report_file 2>&1 + awk '/error|warning/ {print}' $report_file > $report_file_filtered + cat $report_file_filtered + + check_error_log_file $report_file_filtered "Building doxygen documentation: FAILURE." + print_info "Finished doxygen documentation test. SUCCESS." +} + + +set -e +set -x + +if [[ $TRAVIS_COMMIT_MESSAGE == *"[no-build]"* ]]; then + print_info "Option '[no-build]' is detected, sources will not be built, checked, verified and published." + exit 0 +fi + +if [[ $TRAVIS_COMMIT_MESSAGE == *"[build-only-osx]"* ]]; then + if [[ $1 == BUILD_TEST_CCORE_MACOS ]]; then + print_info "Option '[build-only-osx]' is detected, mac os build will be started." + else + print_info "Option '[build-only-osx]' is detected, sources will not be built, checked, verified and published." + exit 0 + fi +fi + +case $1 in + BUILD_CCORE) + run_build_ccore_job ;; + + ANALYSE_CCORE) + run_analyse_ccore_job ;; + + UT_CCORE) + run_ut_ccore_job ;; + + VALGRIND_CCORE) + run_valgrind_ccore_job ;; + + TEST_PYCLUSTERING) + run_test_pyclustering_job ;; + + IT_CCORE_X86) + run_integration_test_job x86 ;; + + IT_CCORE_X64) + run_integration_test_job x64 ;; + + BUILD_TEST_CCORE_MACOS) + run_build_test_ccore_macos_job ;; + + DOCUMENTATION) + run_doxygen_job ;; + + DEPLOY) + run_deploy_job ;; + + *) + print_error "Unknown target is specified: '$1'" + exit 1 ;; +esac diff --git a/cmakedocs.cmake b/cmakedocs.cmake new file mode 100644 index 000000000..95aed7cd4 --- /dev/null +++ b/cmakedocs.cmake @@ -0,0 +1,130 @@ +project(PHiLiP) + +CMAKE_MINIMUM_REQUIRED(VERSION 3.3.0) +set(CMAKE_CXX_COMPILER mpicxx) +set(CMAKE_CXX_STANDARD_REQUIRED ON) +set(CMAKE_CXX_STANDARD 17) + +set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -Wall -Werror") + +set(MPIMAX 4 CACHE STRING "Default number of processors used in ctest mpirun -np MPIMAX. Not the same as ctest -jX") + + +#set(CMAKE_CXX_FLAGS_DEBUG "${CMAKE_CXX_FLAGS_DEBUG} -Og -g") + +set(EXECUTABLE_OUTPUT_PATH ${PROJECT_BINARY_DIR}/bin) +set(CLEAN_UP_FILES ./bin/* ./CMakeCache.txt) + +# Use ld.gold for faster linking +option(USE_LD_GOLD "Use GNU gold linker" OFF) +if(USE_LD_GOLD AND "${CMAKE_C_COMPILER_ID}" STREQUAL "GNU") + execute_process(COMMAND ${CMAKE_C_COMPILER} -fuse-ld=gold -Wl,--version OUTPUT_VARIABLE stdout ERROR_QUIET) + if("${stdout}" MATCHES "GNU gold") + set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -fuse-ld=gold") + set(CMAKE_C_FLAGS "${CMAKE_C_FLAGS} -fuse-ld=gold") + else() + message(WARNING "GNU gold linker isn't available, using the default system linker.") + endif() +endif() + +add_custom_target(1D COMMAND ${CMAKE_COMMAND} -E echo "Makes all 1D executables, including tests. Allows ctest -R 1D") +add_custom_target(2D COMMAND ${CMAKE_COMMAND} -E echo "Makes all 2D executables, including tests. Allows ctest -R 2D") +add_custom_target(3D COMMAND ${CMAKE_COMMAND} -E echo "Makes all 3D executables, including tests. Allows ctest -R 3D") +SET(_dimension_targets "# +# make 1D - to build all 1D executables, including tests, allow 'ctest -R 1D' +# make 2D - to build all 2D executables, including tests, allow 'ctest -R 2D' +# make 3D - to build all 3D executables, including tests, allow 'ctest -R 3D' ") +SET(_philip_targets "# +# make PHiLiP_1D - to build main program (wihtout tests) in 1D +# make PHiLiP_2D - to build main program (wihtout tests) in 2D +# make PHiLiP_3D - to build main program (wihtout tests) in 3D ") +add_custom_target(unit_tests) + +# Source code +include_directories(src) +add_subdirectory(src) + +# Test directory +enable_testing() +add_subdirectory(tests) + +# Documentation +add_subdirectory(doc) + +########## +# Add a few commands to make life easier + +#include(ProcessorCount) +#ProcessorCount(NPROC) + +# Define custom targets to easily switch the build type: +add_custom_target(debug + COMMAND ${CMAKE_COMMAND} -DCMAKE_BUILD_TYPE=Debug ${CMAKE_SOURCE_DIR} + COMMAND ${CMAKE_COMMAND} --build ${CMAKE_BINARY_DIR} --target all + COMMENT "Switch CMAKE_BUILD_TYPE to Debug" + ) + +add_custom_target(release + COMMAND ${CMAKE_COMMAND} -DCMAKE_BUILD_TYPE=Release ${CMAKE_SOURCE_DIR} + COMMAND ${CMAKE_COMMAND} --build ${CMAKE_BINARY_DIR} --target all + COMMENT "Switch CMAKE_BUILD_TYPE to Release" + ) +IF(${DEAL_II_BUILD_TYPE} MATCHES "DebugRelease") +SET(_switch_targets "# +# make debug - to switch the build type to 'Debug' +# make release - to switch the build type to 'Release' " +) +ENDIF() + +# Print out some usage information to file: +FILE(WRITE ${CMAKE_BINARY_DIR}${CMAKE_FILES_DIRECTORY}/print_usage.cmake +"MESSAGE( +\"### +# +# Project ${TARGET} set up with ${DEAL_II_PACKAGE_NAME}-${DEAL_II_PACKAGE_VERSION} found at +# ${DEAL_II_PATH} +# +# CMAKE_BUILD_TYPE: ${CMAKE_BUILD_TYPE} +# +# You can now run +# make - to compile and link all the program and tests. +${_philip_targets} +${_dimension_targets} +${_switch_targets} +") + FILE(APPEND ${CMAKE_BINARY_DIR}${CMAKE_FILES_DIRECTORY}/print_usage.cmake +"# +# make clean - to remove the generated executable as well as +# all intermediate compilation files +# make info - to view this message again +# +# Note that you may append -j N, where N represents the number of processors +# to compile with. If you are compiling on 8GB of memory, be careful not to +# use too many processors, as you will run out of RAM. If you have (m) GB of +# available memory, you can use around (m-1) processor without having to use the swap. +# +###\")" +) + +add_custom_target(info + COMMAND ${CMAKE_COMMAND} -P ${CMAKE_BINARY_DIR}${CMAKE_FILES_DIRECTORY}/print_usage.cmake +) + +# Print this message once: +#IF(NOT USAGE_PRINTED) +# INCLUDE(${CMAKE_BINARY_DIR}${CMAKE_FILES_DIRECTORY}/print_usage.cmake) +# SET(USAGE_PRINTED TRUE CACHE INTERNAL "") +#ELSE() +# MESSAGE(STATUS "Run make info to print a detailed help message") +#ENDIF() +INCLUDE(${CMAKE_BINARY_DIR}${CMAKE_FILES_DIRECTORY}/print_usage.cmake) +MESSAGE(STATUS "Run make info to print a detailed help message") + +MESSAGE(STATUS " +### +# +# Using MPIMAX = ${MPIMAX}, which sets the default values used by ctest for the MPI runs. +# Can use cmake ../ -DMPIMAX=XX to change this default value. +# +### +") diff --git a/doc/doxygen.cfg.in b/doc/doxygen.cfg.in index b59277755..d805e64da 100644 --- a/doc/doxygen.cfg.in +++ b/doc/doxygen.cfg.in @@ -4,9 +4,9 @@ # Project related configuration options #--------------------------------------------------------------------------- DOXYFILE_ENCODING = UTF-8 -PROJECT_NAME = "@CMAKE_PROJECT_NAME@" +PROJECT_NAME = "[P]arallel [Hi]gh-order [Li]brary for [P]DEs" PROJECT_NUMBER = "@PHiLiP_VERSION@" -PROJECT_BRIEF = "(P)arallel (Hi)gh-order (Li)brary for PDEs" +PROJECT_BRIEF = "Parallel High-Order Library for PDEs through hp-adaptive Discontinuous Galerkin methods" PROJECT_LOGO = OUTPUT_DIRECTORY = CREATE_SUBDIRS = NO @@ -40,7 +40,7 @@ BUILTIN_STL_SUPPORT = NO CPP_CLI_SUPPORT = NO SIP_SUPPORT = NO IDL_PROPERTY_SUPPORT = YES -DISTRIBUTE_GROUP_DOC = NO +DISTRIBUTE_GROUP_DOC = YES GROUP_NESTED_COMPOUNDS = NO SUBGROUPING = YES INLINE_GROUPED_CLASSES = NO @@ -50,7 +50,7 @@ LOOKUP_CACHE_SIZE = 0 #--------------------------------------------------------------------------- # Build related configuration options #--------------------------------------------------------------------------- -EXTRACT_ALL = YES +EXTRACT_ALL = NO EXTRACT_PRIVATE = YES EXTRACT_PACKAGE = NO EXTRACT_STATIC = YES @@ -94,7 +94,7 @@ QUIET = NO WARNINGS = YES WARN_IF_UNDOCUMENTED = YES WARN_IF_DOC_ERROR = YES -WARN_NO_PARAMDOC = YES +WARN_NO_PARAMDOC = NO WARN_AS_ERROR = NO WARN_FORMAT = "$file:$line: $text" WARN_LOGFILE = @@ -109,13 +109,14 @@ FILE_PATTERNS = *.cpp \ *.h.in \ *.dox \ # *.md \ # doxygen was given errors for relative paths in markdown files. - *.mk +# *.mk RECURSIVE = YES EXCLUDE = @PROJECT_SOURCE_DIR@/README.md EXCLUDE_SYMLINKS = NO EXCLUDE_PATTERNS = */build/* \ */build_*/* \ */src/deprecated/* \ + */src/dummy/* \ */deprecated/* \ */external/* EXCLUDE_SYMBOLS = diff --git a/doc/doxygen_codedocs.cfg b/doc/doxygen_codedocs.cfg new file mode 100644 index 000000000..fbe523698 --- /dev/null +++ b/doc/doxygen_codedocs.cfg @@ -0,0 +1,330 @@ +# Doxyfile 1.8.11 + +#--------------------------------------------------------------------------- +# Project related configuration options +#--------------------------------------------------------------------------- +DOXYFILE_ENCODING = UTF-8 +PROJECT_NAME = "[P]arallel [Hi]gh-order [Li]brary for [P]DEs" +PROJECT_NUMBER = "Latest" +PROJECT_BRIEF = "Parallel High-Order Library for PDEs through hp-adaptive Discontinuous Galerkin methods" +PROJECT_LOGO = +OUTPUT_DIRECTORY = +CREATE_SUBDIRS = NO +ALLOW_UNICODE_NAMES = NO +OUTPUT_LANGUAGE = English +BRIEF_MEMBER_DESC = YES +REPEAT_BRIEF = YES +ABBREVIATE_BRIEF = +ALWAYS_DETAILED_SEC = NO +INLINE_INHERITED_MEMB = NO +FULL_PATH_NAMES = YES +STRIP_FROM_PATH = ./ +STRIP_FROM_INC_PATH = +SHORT_NAMES = NO +JAVADOC_AUTOBRIEF = NO +QT_AUTOBRIEF = NO +MULTILINE_CPP_IS_BRIEF = NO +INHERIT_DOCS = YES +SEPARATE_MEMBER_PAGES = NO +TAB_SIZE = 4 +ALIASES = +TCL_SUBST = +OPTIMIZE_OUTPUT_FOR_C = NO +OPTIMIZE_OUTPUT_JAVA = NO +OPTIMIZE_FOR_FORTRAN = NO +OPTIMIZE_OUTPUT_VHDL = NO +EXTENSION_MAPPING = +MARKDOWN_SUPPORT = YES +AUTOLINK_SUPPORT = YES +BUILTIN_STL_SUPPORT = NO +CPP_CLI_SUPPORT = NO +SIP_SUPPORT = NO +IDL_PROPERTY_SUPPORT = YES +DISTRIBUTE_GROUP_DOC = YES +GROUP_NESTED_COMPOUNDS = NO +SUBGROUPING = YES +INLINE_GROUPED_CLASSES = NO +INLINE_SIMPLE_STRUCTS = NO +TYPEDEF_HIDES_STRUCT = NO +LOOKUP_CACHE_SIZE = 0 +#--------------------------------------------------------------------------- +# Build related configuration options +#--------------------------------------------------------------------------- +EXTRACT_ALL = NO +EXTRACT_PRIVATE = YES +EXTRACT_PACKAGE = NO +EXTRACT_STATIC = YES +EXTRACT_LOCAL_CLASSES = YES +EXTRACT_LOCAL_METHODS = NO +EXTRACT_ANON_NSPACES = NO +HIDE_UNDOC_MEMBERS = NO +HIDE_UNDOC_CLASSES = NO +HIDE_FRIEND_COMPOUNDS = NO +HIDE_IN_BODY_DOCS = NO +INTERNAL_DOCS = NO +CASE_SENSE_NAMES = YES +HIDE_SCOPE_NAMES = NO +HIDE_COMPOUND_REFERENCE= NO +SHOW_INCLUDE_FILES = YES +SHOW_GROUPED_MEMB_INC = NO +FORCE_LOCAL_INCLUDES = NO +INLINE_INFO = YES +SORT_MEMBER_DOCS = YES +SORT_BRIEF_DOCS = NO +SORT_MEMBERS_CTORS_1ST = NO +SORT_GROUP_NAMES = NO +SORT_BY_SCOPE_NAME = NO +STRICT_PROTO_MATCHING = NO +GENERATE_TODOLIST = YES +GENERATE_TESTLIST = YES +GENERATE_BUGLIST = YES +GENERATE_DEPRECATEDLIST= YES +ENABLED_SECTIONS = +MAX_INITIALIZER_LINES = 30 +SHOW_USED_FILES = YES +SHOW_FILES = YES +SHOW_NAMESPACES = YES +FILE_VERSION_FILTER = +LAYOUT_FILE = ./doc/DoxygenLayout.xml +CITE_BIB_FILES = ./doc/code.bib +#--------------------------------------------------------------------------- +# Configuration options related to warning and progress messages +#--------------------------------------------------------------------------- +QUIET = NO +WARNINGS = YES +WARN_IF_UNDOCUMENTED = YES +WARN_IF_DOC_ERROR = YES +WARN_NO_PARAMDOC = NO +WARN_AS_ERROR = NO +WARN_FORMAT = "$file:$line: $text" +WARN_LOGFILE = +#--------------------------------------------------------------------------- +# Configuration options related to the input files +#--------------------------------------------------------------------------- +INPUT = ./ +INPUT_ENCODING = UTF-8 +FILE_PATTERNS = *.cpp \ + *.h \ + *.hpp \ + *.h.in \ + *.dox \ +# *.md \ # doxygen was given errors for relative paths in markdown files. +# *.mk +RECURSIVE = YES +EXCLUDE = ./README.md +EXCLUDE_SYMLINKS = NO +EXCLUDE_PATTERNS = */build/* \ + */build_*/* \ + */src/deprecated/* \ + */src/dummy/* \ + */deprecated/* \ + */external/* +EXCLUDE_SYMBOLS = +EXAMPLE_PATH = +EXAMPLE_PATTERNS = +EXAMPLE_RECURSIVE = NO +IMAGE_PATH = +INPUT_FILTER = +FILTER_PATTERNS = +FILTER_SOURCE_FILES = NO +FILTER_SOURCE_PATTERNS = +USE_MDFILE_AS_MAINPAGE = +#--------------------------------------------------------------------------- +# Configuration options related to source browsing +#--------------------------------------------------------------------------- +SOURCE_BROWSER = YES +INLINE_SOURCES = NO +STRIP_CODE_COMMENTS = YES +REFERENCED_BY_RELATION = NO +REFERENCES_RELATION = NO +REFERENCES_LINK_SOURCE = YES +SOURCE_TOOLTIPS = YES +USE_HTAGS = NO +VERBATIM_HEADERS = YES +CLANG_ASSISTED_PARSING = NO +CLANG_OPTIONS = +#--------------------------------------------------------------------------- +# Configuration options related to the alphabetical class index +#--------------------------------------------------------------------------- +ALPHABETICAL_INDEX = YES +COLS_IN_ALPHA_INDEX = 5 +IGNORE_PREFIX = +#--------------------------------------------------------------------------- +# Configuration options related to the HTML output +#--------------------------------------------------------------------------- +GENERATE_HTML = YES +HTML_OUTPUT = html +HTML_FILE_EXTENSION = .html +HTML_HEADER = +HTML_FOOTER = +HTML_STYLESHEET = +HTML_EXTRA_STYLESHEET = +HTML_EXTRA_FILES = +HTML_COLORSTYLE_HUE = 240 +HTML_COLORSTYLE_SAT = 255 +HTML_COLORSTYLE_GAMMA = 40 +HTML_TIMESTAMP = YES +HTML_DYNAMIC_SECTIONS = NO +HTML_INDEX_NUM_ENTRIES = 100 +GENERATE_DOCSET = NO +DOCSET_FEEDNAME = "Doxygen generated docs" +DOCSET_BUNDLE_ID = org.doxygen.Project +DOCSET_PUBLISHER_ID = org.doxygen.Publisher +DOCSET_PUBLISHER_NAME = Publisher +GENERATE_HTMLHELP = NO +CHM_FILE = +HHC_LOCATION = +GENERATE_CHI = NO +CHM_INDEX_ENCODING = +BINARY_TOC = NO +TOC_EXPAND = NO +GENERATE_QHP = NO +QCH_FILE = +QHP_NAMESPACE = org.doxygen.Project +QHP_VIRTUAL_FOLDER = doc +QHP_CUST_FILTER_NAME = +QHP_CUST_FILTER_ATTRS = +QHP_SECT_FILTER_ATTRS = +QHG_LOCATION = +GENERATE_ECLIPSEHELP = NO +ECLIPSE_DOC_ID = org.doxygen.Project +DISABLE_INDEX = NO +GENERATE_TREEVIEW = YES +ENUM_VALUES_PER_LINE = 4 +TREEVIEW_WIDTH = 250 +EXT_LINKS_IN_WINDOW = NO +FORMULA_FONTSIZE = 10 +FORMULA_TRANSPARENT = YES +USE_MATHJAX = YES +MATHJAX_FORMAT = HTML-CSS +MATHJAX_RELPATH = http://cdn.mathjax.org/mathjax/latest +MATHJAX_EXTENSIONS = +MATHJAX_CODEFILE = +SEARCHENGINE = YES +SERVER_BASED_SEARCH = NO +EXTERNAL_SEARCH = NO +SEARCHENGINE_URL = +SEARCHDATA_FILE = searchdata.xml +EXTERNAL_SEARCH_ID = +EXTRA_SEARCH_MAPPINGS = +#--------------------------------------------------------------------------- +# Configuration options related to the LaTeX output +#--------------------------------------------------------------------------- +GENERATE_LATEX = YES +LATEX_OUTPUT = latex +LATEX_CMD_NAME = latex +MAKEINDEX_CMD_NAME = makeindex +COMPACT_LATEX = NO +PAPER_TYPE = a4 +EXTRA_PACKAGES = +LATEX_HEADER = +LATEX_FOOTER = +LATEX_EXTRA_STYLESHEET = +LATEX_EXTRA_FILES = +PDF_HYPERLINKS = YES +USE_PDFLATEX = YES +LATEX_BATCHMODE = NO +LATEX_HIDE_INDICES = NO +LATEX_SOURCE_CODE = NO +LATEX_BIB_STYLE = plain +LATEX_TIMESTAMP = YES +#--------------------------------------------------------------------------- +# Configuration options related to the RTF output +#--------------------------------------------------------------------------- +GENERATE_RTF = NO +RTF_OUTPUT = rtf +COMPACT_RTF = NO +RTF_HYPERLINKS = NO +RTF_STYLESHEET_FILE = +RTF_EXTENSIONS_FILE = +RTF_SOURCE_CODE = NO +#--------------------------------------------------------------------------- +# Configuration options related to the man page output +#--------------------------------------------------------------------------- +GENERATE_MAN = NO +MAN_OUTPUT = man +MAN_EXTENSION = .3 +MAN_SUBDIR = +MAN_LINKS = NO +#--------------------------------------------------------------------------- +# Configuration options related to the XML output +#--------------------------------------------------------------------------- +GENERATE_XML = NO +XML_OUTPUT = xml +XML_PROGRAMLISTING = YES +#--------------------------------------------------------------------------- +# Configuration options related to the DOCBOOK output +#--------------------------------------------------------------------------- +GENERATE_DOCBOOK = NO +DOCBOOK_OUTPUT = docbook +DOCBOOK_PROGRAMLISTING = NO +#--------------------------------------------------------------------------- +# Configuration options for the AutoGen Definitions output +#--------------------------------------------------------------------------- +GENERATE_AUTOGEN_DEF = NO +#--------------------------------------------------------------------------- +# Configuration options related to the Perl module output +#--------------------------------------------------------------------------- +GENERATE_PERLMOD = NO +PERLMOD_LATEX = NO +PERLMOD_PRETTY = YES +PERLMOD_MAKEVAR_PREFIX = +#--------------------------------------------------------------------------- +# Configuration options related to the preprocessor +#--------------------------------------------------------------------------- +ENABLE_PREPROCESSING = YES +MACRO_EXPANSION = NO +EXPAND_ONLY_PREDEF = NO +SEARCH_INCLUDES = YES +INCLUDE_PATH = +INCLUDE_FILE_PATTERNS = +PREDEFINED = TYPE_RC:=2 TYPE_COMPLEX:=2 +EXPAND_AS_DEFINED = +SKIP_FUNCTION_MACROS = YES +#--------------------------------------------------------------------------- +# Configuration options related to external references +#--------------------------------------------------------------------------- +TAGFILES = +GENERATE_TAGFILE = +ALLEXTERNALS = NO +EXTERNAL_GROUPS = YES +EXTERNAL_PAGES = YES +PERL_PATH = /usr/bin/perl +#--------------------------------------------------------------------------- +# Configuration options related to the dot tool +#--------------------------------------------------------------------------- +CLASS_DIAGRAMS = YES +MSCGEN_PATH = +DIA_PATH = +HIDE_UNDOC_RELATIONS = NO +HAVE_DOT = @DOXYGEN_DOT_FOUND@ +DOT_NUM_THREADS = 0 +DOT_FONTNAME = Helvetica +DOT_FONTSIZE = 10 +DOT_FONTPATH = +CLASS_GRAPH = YES +COLLABORATION_GRAPH = YES +GROUP_GRAPHS = YES +UML_LOOK = YES +UML_LIMIT_NUM_FIELDS = 100 +TEMPLATE_RELATIONS = YES +INCLUDE_GRAPH = YES +INCLUDED_BY_GRAPH = YES +CALL_GRAPH = NO +CALLER_GRAPH = NO +GRAPHICAL_HIERARCHY = YES +DIRECTORY_GRAPH = YES +DOT_IMAGE_FORMAT = svg +INTERACTIVE_SVG = NO +DOT_PATH = @DOXYGEN_DOT_EXECUTABLE@ +DOTFILE_DIRS = +MSCFILE_DIRS = +DIAFILE_DIRS = +PLANTUML_JAR_PATH = +PLANTUML_INCLUDE_PATH = +DOT_GRAPH_MAX_NODES = 50 +MAX_DOT_GRAPH_DEPTH = 0 +DOT_TRANSPARENT = NO +DOT_MULTI_TARGETS = NO +GENERATE_LEGEND = YES +DOT_CLEANUP = YES diff --git a/doc/getting_started.dox b/doc/getting_started.dox deleted file mode 100644 index acdea7549..000000000 --- a/doc/getting_started.dox +++ /dev/null @@ -1,7 +0,0 @@ -/** \file - -\getting_started - - -*/ - diff --git a/doc/mainpage.dox b/doc/mainpage.dox index 6b0ca6a99..673a98927 100644 --- a/doc/mainpage.dox +++ b/doc/mainpage.dox @@ -2,12 +2,7 @@ \mainpage -These pages serve as the main programmer's reference manual for PHiLiP and were automatically generated from the -source using Doxygen. - - -A good place to start is the Getting Started page, which provides a getting started -file which contains detailed instructions for building your understanding of the implementation. +These pages serve as the main programmer's reference manual for PHiLiP and were automatically generated from the source using Doxygen. The project home contains additional information relating to the code functionality and status. diff --git a/passing_tests.log b/passing_tests.log new file mode 100644 index 000000000..e76d54fff --- /dev/null +++ b/passing_tests.log @@ -0,0 +1,232 @@ +2019/10/22 fb80cd060f51fdf69034fe4f9b2ee7a5652eb5bb + +Test project /home/ddong/Codes/PHiLiP/build_release + Start 1: 1D_numerical_flux_conservation + 1/53 Test #1: 1D_numerical_flux_conservation ......................................... Passed 0.57 sec + Start 2: 2D_numerical_flux_conservation + 2/53 Test #2: 2D_numerical_flux_conservation ......................................... Passed 0.44 sec + Start 3: 3D_numerical_flux_conservation + 3/53 Test #3: 3D_numerical_flux_conservation ......................................... Passed 0.43 sec + Start 4: 1D_jacobian_matrix_regression + 4/53 Test #4: 1D_jacobian_matrix_regression .......................................... Passed 0.50 sec + Start 5: 2D_jacobian_matrix_regression + 5/53 Test #5: 2D_jacobian_matrix_regression .......................................... Passed 0.54 sec + Start 6: 3D_jacobian_matrix_regression + 6/53 Test #6: 3D_jacobian_matrix_regression .......................................... Passed 2.68 sec + Start 7: 1D_euler_convert_primitive_conservative + 7/53 Test #7: 1D_euler_convert_primitive_conservative ................................ Passed 0.22 sec + Start 8: 2D_euler_convert_primitive_conservative + 8/53 Test #8: 2D_euler_convert_primitive_conservative ................................ Passed 0.22 sec + Start 9: 3D_euler_convert_primitive_conservative + 9/53 Test #9: 3D_euler_convert_primitive_conservative ................................ Passed 0.65 sec + Start 10: 1D_euler_manufactured_solution_source +10/53 Test #10: 1D_euler_manufactured_solution_source .................................. Passed 0.21 sec + Start 11: 2D_euler_manufactured_solution_source +11/53 Test #11: 2D_euler_manufactured_solution_source .................................. Passed 0.23 sec + Start 12: 3D_euler_manufactured_solution_source +12/53 Test #12: 3D_euler_manufactured_solution_source .................................. Passed 0.72 sec + Start 13: 1D_euler_convective_jacobian +13/53 Test #13: 1D_euler_convective_jacobian ........................................... Passed 0.23 sec + Start 14: 2D_euler_convective_jacobian +14/53 Test #14: 2D_euler_convective_jacobian ........................................... Passed 0.25 sec + Start 15: 3D_euler_convective_jacobian +15/53 Test #15: 3D_euler_convective_jacobian ........................................... Passed 1.41 sec + Start 16: 2D_HighOrder_MappingFEField +16/53 Test #16: 2D_HighOrder_MappingFEField ............................................ Passed 0.65 sec + Start 17: 3D_HighOrder_MappingFEField +17/53 Test #17: 3D_HighOrder_MappingFEField ............................................ Passed 4.54 sec + Start 18: 2D_RBF_mesh_movement +18/53 Test #18: 2D_RBF_mesh_movement ................................................... Passed 0.63 sec + Start 19: 3D_RBF_mesh_movement +19/53 Test #19: 3D_RBF_mesh_movement ................................................... Passed 2.06 sec + Start 20: 2D_make_cells_valid +20/53 Test #20: 2D_make_cells_valid .................................................... Passed 0.63 sec + Start 21: 1D_ADVECTION_IMPLICIT_MANUFACTURED_SOLUTION +21/53 Test #21: 1D_ADVECTION_IMPLICIT_MANUFACTURED_SOLUTION ............................ Passed 0.53 sec + Start 22: 2D_ADVECTION_IMPLICIT_MANUFACTURED_SOLUTION +22/53 Test #22: 2D_ADVECTION_IMPLICIT_MANUFACTURED_SOLUTION ............................ Passed 3.72 sec + Start 23: MPI_2D_ADVECTION_EXPLICIT_MANUFACTURED_SOLUTION_LONG +23/53 Test #23: MPI_2D_ADVECTION_EXPLICIT_MANUFACTURED_SOLUTION_LONG ................... Passed 136.48 sec + Start 24: MPI_2D_ADVECTION_IMPLICIT_MANUFACTURED_SOLUTION +24/53 Test #24: MPI_2D_ADVECTION_IMPLICIT_MANUFACTURED_SOLUTION ........................ Passed 2.71 sec + Start 25: MPI_3D_ADVECTION_IMPLICIT_MANUFACTURED_SOLUTION +25/53 Test #25: MPI_3D_ADVECTION_IMPLICIT_MANUFACTURED_SOLUTION ........................ Passed 11.19 sec + Start 26: 1D_ADVECTION_IMPLICIT_MANUFACTURED_SOLUTION_STRONG +26/53 Test #26: 1D_ADVECTION_IMPLICIT_MANUFACTURED_SOLUTION_STRONG ..................... Passed 0.58 sec + Start 27: 2D_ADVECTION_IMPLICIT_MANUFACTURED_SOLUTION_STRONG +27/53 Test #27: 2D_ADVECTION_IMPLICIT_MANUFACTURED_SOLUTION_STRONG ..................... Passed 6.27 sec + Start 28: MPI_3D_ADVECTION_IMPLICIT_MANUFACTURED_SOLUTION_STRONG_MEDIUM +28/53 Test #28: MPI_3D_ADVECTION_IMPLICIT_MANUFACTURED_SOLUTION_STRONG_MEDIUM .......... Passed 35.05 sec + Start 29: 1D_ADVECTION_IMPLICIT_COLLOCATED_MANUFACTURED_SOLUTION +29/53 Test #29: 1D_ADVECTION_IMPLICIT_COLLOCATED_MANUFACTURED_SOLUTION ................. Passed 0.62 sec + Start 30: 1D_DIFFUSION_IMPLICIT_MANUFACTURED_SOLUTION +30/53 Test #30: 1D_DIFFUSION_IMPLICIT_MANUFACTURED_SOLUTION ............................ Passed 0.59 sec + Start 31: 2D_DIFFUSION_IMPLICIT_MANUFACTURED_SOLUTION +31/53 Test #31: 2D_DIFFUSION_IMPLICIT_MANUFACTURED_SOLUTION ............................ Passed 1.64 sec + Start 32: MPI_3D_DIFFUSION_IMPLICIT_MANUFACTURED_SOLUTION_MEDIUM +32/53 Test #32: MPI_3D_DIFFUSION_IMPLICIT_MANUFACTURED_SOLUTION_MEDIUM ................. Passed 42.64 sec + Start 33: 1D_CONVECTION_DIFFUSION_IMPLICIT_MANUFACTURED_SOLUTION +33/53 Test #33: 1D_CONVECTION_DIFFUSION_IMPLICIT_MANUFACTURED_SOLUTION ................. Passed 0.67 sec + Start 34: 2D_CONVECTION_DIFFUSION_IMPLICIT_MANUFACTURED_SOLUTION +34/53 Test #34: 2D_CONVECTION_DIFFUSION_IMPLICIT_MANUFACTURED_SOLUTION ................. Passed 1.73 sec + Start 35: 2D_CONVECTION_DIFFUSION_IMPLICIT_COLLOCATED_MANUFACTURED_SOLUTION +35/53 Test #35: 2D_CONVECTION_DIFFUSION_IMPLICIT_COLLOCATED_MANUFACTURED_SOLUTION ...... Passed 1.80 sec + Start 36: MPI_3D_CONVECTION_DIFFUSION_IMPLICIT_MANUFACTURED_SOLUTION_MEDIUM +36/53 Test #36: MPI_3D_CONVECTION_DIFFUSION_IMPLICIT_MANUFACTURED_SOLUTION_MEDIUM ...... Passed 40.54 sec + Start 37: 1D_ADVECTION_VECTOR_VALUED_IMPLICIT_MANUFACTURED_SOLUTION +37/53 Test #37: 1D_ADVECTION_VECTOR_VALUED_IMPLICIT_MANUFACTURED_SOLUTION .............. Passed 0.68 sec + Start 38: 2D_ADVECTION_VECTOR_VALUED_IMPLICIT_MANUFACTURED_SOLUTION +38/53 Test #38: 2D_ADVECTION_VECTOR_VALUED_IMPLICIT_MANUFACTURED_SOLUTION .............. Passed 9.39 sec + Start 39: MPI_3D_ADVECTION_VECTOR_VALUED_IMPLICIT_MANUFACTURED_SOLUTION_MEDIUM +39/53 Test #39: MPI_3D_ADVECTION_VECTOR_VALUED_IMPLICIT_MANUFACTURED_SOLUTION_MEDIUM ... Passed 16.11 sec + Start 40: 1D_BURGERS_INVISCID_IMPLICIT_MANUFACTURED_SOLUTION +40/53 Test #40: 1D_BURGERS_INVISCID_IMPLICIT_MANUFACTURED_SOLUTION ..................... Passed 0.65 sec + Start 41: 1D_BURGERS_INVISCID_IMPLICIT_COLLOCATED_MANUFACTURED_SOLUTION +41/53 Test #41: 1D_BURGERS_INVISCID_IMPLICIT_COLLOCATED_MANUFACTURED_SOLUTION .......... Passed 0.61 sec + Start 42: 1D_burgers_stability_LONG +42/53 Test #42: 1D_burgers_stability_LONG .............................................. Passed 198.05 sec + Start 43: MPI_3D_EULER_SPLIT_TAYLOR_GREEN +43/53 Test #43: MPI_3D_EULER_SPLIT_TAYLOR_GREEN ........................................***Failed 0.98 sec + Start 44: MPI_2D_ADVECTION_EXPLICIT_PERIODIC_LONG +44/53 Test #44: MPI_2D_ADVECTION_EXPLICIT_PERIODIC_LONG ................................ Passed 24.06 sec + Start 45: 1D_EULER_LAXFRIEDRICHS_MANUFACTURED_SOLUTION +45/53 Test #45: 1D_EULER_LAXFRIEDRICHS_MANUFACTURED_SOLUTION ........................... Passed 0.94 sec + Start 46: 1D_EULER_LAXFRIEDRICHS_COLLOCATED_MANUFACTURED_SOLUTION +46/53 Test #46: 1D_EULER_LAXFRIEDRICHS_COLLOCATED_MANUFACTURED_SOLUTION ................ Passed 0.91 sec + Start 47: 2D_EULER_LAXFRIEDRICHS_MANUFACTURED_SOLUTION_LONG +47/53 Test #47: 2D_EULER_LAXFRIEDRICHS_MANUFACTURED_SOLUTION_LONG ...................... Passed 27.95 sec + Start 48: 2D_EULER_LAXFRIEDRICHS_COLLOCATED_MANUFACTURED_SOLUTION +48/53 Test #48: 2D_EULER_LAXFRIEDRICHS_COLLOCATED_MANUFACTURED_SOLUTION ................***Failed 35.32 sec + Start 49: MPI_2D_EULER_LAXFRIEDRICHS_MANUFACTURED_SOLUTION_MEDIUM +49/53 Test #49: MPI_2D_EULER_LAXFRIEDRICHS_MANUFACTURED_SOLUTION_MEDIUM ................ Passed 11.77 sec + Start 50: 1D_EULER_ROE_MANUFACTURED_SOLUTION_LONG +50/53 Test #50: 1D_EULER_ROE_MANUFACTURED_SOLUTION_LONG ................................ Passed 0.91 sec + Start 51: MPI_2D_EULER_ROE_MANUFACTURED_SOLUTION_MEDIUM +51/53 Test #51: MPI_2D_EULER_ROE_MANUFACTURED_SOLUTION_MEDIUM .......................... Passed 10.06 sec + Start 52: MPI_2D_EULER_INTEGRATION_CYLINDER_LONG +52/53 Test #52: MPI_2D_EULER_INTEGRATION_CYLINDER_LONG ................................. Passed 224.97 sec + Start 53: MPI_2D_EULER_INTEGRATION_GAUSSIAN_BUMP_LONG +53/53 Test #53: MPI_2D_EULER_INTEGRATION_GAUSSIAN_BUMP_LONG ............................ Passed 707.42 sec + +96% tests passed, 2 tests failed out of 53 + +Total Test time (real) = 1576.04 sec + +The following tests FAILED: + 43 - MPI_3D_EULER_SPLIT_TAYLOR_GREEN (Failed) + 48 - 2D_EULER_LAXFRIEDRICHS_COLLOCATED_MANUFACTURED_SOLUTION (Failed) +Errors while running CTest + +2019/10/08 1a3972a + +Test project /home/ddong/Codes/PHiLiP/build_release + Start 1: 1D_numerical_flux_conservation + 1/50 Test #1: 1D_numerical_flux_conservation ......................................... Passed 0.57 sec + Start 2: 2D_numerical_flux_conservation + 2/50 Test #2: 2D_numerical_flux_conservation ......................................... Passed 0.42 sec + Start 3: 3D_numerical_flux_conservation + 3/50 Test #3: 3D_numerical_flux_conservation ......................................... Passed 0.43 sec + Start 4: 1D_jacobian_matrix_regression + 4/50 Test #4: 1D_jacobian_matrix_regression .......................................... Passed 0.48 sec + Start 5: 2D_jacobian_matrix_regression + 5/50 Test #5: 2D_jacobian_matrix_regression .......................................... Passed 0.53 sec + Start 6: 3D_jacobian_matrix_regression + 6/50 Test #6: 3D_jacobian_matrix_regression .......................................... Passed 2.77 sec + Start 7: 1D_euler_convert_primitive_conservative + 7/50 Test #7: 1D_euler_convert_primitive_conservative ................................ Passed 0.22 sec + Start 8: 2D_euler_convert_primitive_conservative + 8/50 Test #8: 2D_euler_convert_primitive_conservative ................................ Passed 0.21 sec + Start 9: 3D_euler_convert_primitive_conservative + 9/50 Test #9: 3D_euler_convert_primitive_conservative ................................ Passed 0.65 sec + Start 10: 1D_euler_manufactured_solution_source +10/50 Test #10: 1D_euler_manufactured_solution_source .................................. Passed 0.19 sec + Start 11: 2D_euler_manufactured_solution_source +11/50 Test #11: 2D_euler_manufactured_solution_source .................................. Passed 0.22 sec + Start 12: 3D_euler_manufactured_solution_source +12/50 Test #12: 3D_euler_manufactured_solution_source .................................. Passed 0.70 sec + Start 13: 1D_euler_convective_jacobian +13/50 Test #13: 1D_euler_convective_jacobian ........................................... Passed 0.20 sec + Start 14: 2D_euler_convective_jacobian +14/50 Test #14: 2D_euler_convective_jacobian ........................................... Passed 0.23 sec + Start 15: 3D_euler_convective_jacobian +15/50 Test #15: 3D_euler_convective_jacobian ........................................... Passed 1.44 sec + Start 16: 2D_HighOrder_MappingFEField +16/50 Test #16: 2D_HighOrder_MappingFEField ............................................ Passed 0.59 sec + Start 17: 3D_HighOrder_MappingFEField +17/50 Test #17: 3D_HighOrder_MappingFEField ............................................ Passed 2.70 sec + Start 18: 1D_ADVECTION_IMPLICIT_MANUFACTURED_SOLUTION +18/50 Test #18: 1D_ADVECTION_IMPLICIT_MANUFACTURED_SOLUTION ............................ Passed 0.65 sec + Start 19: 2D_ADVECTION_IMPLICIT_MANUFACTURED_SOLUTION +19/50 Test #19: 2D_ADVECTION_IMPLICIT_MANUFACTURED_SOLUTION ............................ Passed 3.63 sec + Start 20: MPI_2D_ADVECTION_EXPLICIT_MANUFACTURED_SOLUTION_LONG +20/50 Test #20: MPI_2D_ADVECTION_EXPLICIT_MANUFACTURED_SOLUTION_LONG ................... Passed 133.17 sec + Start 21: MPI_2D_ADVECTION_IMPLICIT_MANUFACTURED_SOLUTION +21/50 Test #21: MPI_2D_ADVECTION_IMPLICIT_MANUFACTURED_SOLUTION ........................ Passed 2.28 sec + Start 22: MPI_3D_ADVECTION_IMPLICIT_MANUFACTURED_SOLUTION +22/50 Test #22: MPI_3D_ADVECTION_IMPLICIT_MANUFACTURED_SOLUTION ........................ Passed 9.68 sec + Start 23: 1D_ADVECTION_IMPLICIT_MANUFACTURED_SOLUTION_STRONG +23/50 Test #23: 1D_ADVECTION_IMPLICIT_MANUFACTURED_SOLUTION_STRONG ..................... Passed 0.55 sec + Start 24: 2D_ADVECTION_IMPLICIT_MANUFACTURED_SOLUTION_STRONG +24/50 Test #24: 2D_ADVECTION_IMPLICIT_MANUFACTURED_SOLUTION_STRONG ..................... Passed 5.86 sec + Start 25: MPI_3D_ADVECTION_IMPLICIT_MANUFACTURED_SOLUTION_STRONG_MEDIUM +25/50 Test #25: MPI_3D_ADVECTION_IMPLICIT_MANUFACTURED_SOLUTION_STRONG_MEDIUM .......... Passed 29.90 sec + Start 26: 1D_ADVECTION_IMPLICIT_COLLOCATED_MANUFACTURED_SOLUTION +26/50 Test #26: 1D_ADVECTION_IMPLICIT_COLLOCATED_MANUFACTURED_SOLUTION ................. Passed 0.54 sec + Start 27: 1D_DIFFUSION_IMPLICIT_MANUFACTURED_SOLUTION +27/50 Test #27: 1D_DIFFUSION_IMPLICIT_MANUFACTURED_SOLUTION ............................ Passed 0.51 sec + Start 28: 2D_DIFFUSION_IMPLICIT_MANUFACTURED_SOLUTION +28/50 Test #28: 2D_DIFFUSION_IMPLICIT_MANUFACTURED_SOLUTION ............................ Passed 1.47 sec + Start 29: MPI_3D_DIFFUSION_IMPLICIT_MANUFACTURED_SOLUTION_MEDIUM +29/50 Test #29: MPI_3D_DIFFUSION_IMPLICIT_MANUFACTURED_SOLUTION_MEDIUM ................. Passed 31.92 sec + Start 30: 1D_CONVECTION_DIFFUSION_IMPLICIT_MANUFACTURED_SOLUTION +30/50 Test #30: 1D_CONVECTION_DIFFUSION_IMPLICIT_MANUFACTURED_SOLUTION ................. Passed 0.54 sec + Start 31: 2D_CONVECTION_DIFFUSION_IMPLICIT_MANUFACTURED_SOLUTION +31/50 Test #31: 2D_CONVECTION_DIFFUSION_IMPLICIT_MANUFACTURED_SOLUTION ................. Passed 1.52 sec + Start 32: 2D_CONVECTION_DIFFUSION_IMPLICIT_COLLOCATED_MANUFACTURED_SOLUTION +32/50 Test #32: 2D_CONVECTION_DIFFUSION_IMPLICIT_COLLOCATED_MANUFACTURED_SOLUTION ...... Passed 1.51 sec + Start 33: MPI_3D_CONVECTION_DIFFUSION_IMPLICIT_MANUFACTURED_SOLUTION_MEDIUM +33/50 Test #33: MPI_3D_CONVECTION_DIFFUSION_IMPLICIT_MANUFACTURED_SOLUTION_MEDIUM ...... Passed 30.37 sec + Start 34: 1D_ADVECTION_VECTOR_VALUED_IMPLICIT_MANUFACTURED_SOLUTION +34/50 Test #34: 1D_ADVECTION_VECTOR_VALUED_IMPLICIT_MANUFACTURED_SOLUTION .............. Passed 0.57 sec + Start 35: 2D_ADVECTION_VECTOR_VALUED_IMPLICIT_MANUFACTURED_SOLUTION +35/50 Test #35: 2D_ADVECTION_VECTOR_VALUED_IMPLICIT_MANUFACTURED_SOLUTION .............. Passed 7.53 sec + Start 36: MPI_3D_ADVECTION_VECTOR_VALUED_IMPLICIT_MANUFACTURED_SOLUTION_MEDIUM +36/50 Test #36: MPI_3D_ADVECTION_VECTOR_VALUED_IMPLICIT_MANUFACTURED_SOLUTION_MEDIUM ... Passed 13.56 sec + Start 37: 1D_BURGERS_INVISCID_IMPLICIT_MANUFACTURED_SOLUTION +37/50 Test #37: 1D_BURGERS_INVISCID_IMPLICIT_MANUFACTURED_SOLUTION ..................... Passed 0.60 sec + Start 38: 1D_BURGERS_INVISCID_IMPLICIT_COLLOCATED_MANUFACTURED_SOLUTION +38/50 Test #38: 1D_BURGERS_INVISCID_IMPLICIT_COLLOCATED_MANUFACTURED_SOLUTION .......... Passed 0.50 sec + Start 39: 1D_burgers_stability_LONG +39/50 Test #39: 1D_burgers_stability_LONG .............................................. Passed 186.27 sec + Start 40: MPI_3D_EULER_SPLIT_TAYLOR_GREEN +40/50 Test #40: MPI_3D_EULER_SPLIT_TAYLOR_GREEN ........................................***Failed 0.99 sec + Start 41: MPI_2D_ADVECTION_EXPLICIT_PERIODIC_LONG +41/50 Test #41: MPI_2D_ADVECTION_EXPLICIT_PERIODIC_LONG ................................ Passed 24.51 sec + Start 42: 1D_EULER_LAXFRIEDRICHS_MANUFACTURED_SOLUTION +42/50 Test #42: 1D_EULER_LAXFRIEDRICHS_MANUFACTURED_SOLUTION ........................... Passed 0.95 sec + Start 43: 1D_EULER_LAXFRIEDRICHS_COLLOCATED_MANUFACTURED_SOLUTION +43/50 Test #43: 1D_EULER_LAXFRIEDRICHS_COLLOCATED_MANUFACTURED_SOLUTION ................ Passed 0.91 sec + Start 44: 2D_EULER_LAXFRIEDRICHS_MANUFACTURED_SOLUTION_LONG +44/50 Test #44: 2D_EULER_LAXFRIEDRICHS_MANUFACTURED_SOLUTION_LONG ...................... Passed 27.54 sec + Start 45: 2D_EULER_LAXFRIEDRICHS_COLLOCATED_MANUFACTURED_SOLUTION +45/50 Test #45: 2D_EULER_LAXFRIEDRICHS_COLLOCATED_MANUFACTURED_SOLUTION ................***Failed 28.02 sec + Start 46: MPI_2D_EULER_LAXFRIEDRICHS_MANUFACTURED_SOLUTION_MEDIUM +46/50 Test #46: MPI_2D_EULER_LAXFRIEDRICHS_MANUFACTURED_SOLUTION_MEDIUM ................ Passed 12.20 sec + Start 47: 1D_EULER_ROE_MANUFACTURED_SOLUTION_LONG +47/50 Test #47: 1D_EULER_ROE_MANUFACTURED_SOLUTION_LONG ................................ Passed 0.97 sec + Start 48: MPI_2D_EULER_ROE_MANUFACTURED_SOLUTION_MEDIUM +48/50 Test #48: MPI_2D_EULER_ROE_MANUFACTURED_SOLUTION_MEDIUM .......................... Passed 8.93 sec + Start 49: MPI_2D_EULER_INTEGRATION_CYLINDER_LONG +49/50 Test #49: MPI_2D_EULER_INTEGRATION_CYLINDER_LONG ................................. Passed 213.84 sec + Start 50: MPI_2D_EULER_INTEGRATION_GAUSSIAN_BUMP_LONG +50/50 Test #50: MPI_2D_EULER_INTEGRATION_GAUSSIAN_BUMP_LONG ............................ Passed 653.68 sec + +96% tests passed, 2 tests failed out of 50 + +Total Test time (real) = 1448.91 sec + +The following tests FAILED: + 40 - MPI_3D_EULER_SPLIT_TAYLOR_GREEN (Failed) + 45 - 2D_EULER_LAXFRIEDRICHS_COLLOCATED_MANUFACTURED_SOLUTION (Failed) +Errors while running CTest + diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 60e30a173..be1ceab53 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -1,8 +1,6 @@ -#set(PHILIP_SRC -# grid_study.cpp) - set(MAIN_SRC main.cpp) +add_subdirectory(dummy) add_subdirectory(linear_solver) add_subdirectory(parameters) add_subdirectory(physics) @@ -67,46 +65,12 @@ foreach(dim RANGE 1 3) target_link_libraries(${MAIN_TARGET} ${ParameterLib}) target_link_libraries(${MAIN_TARGET} ${TestsLib}) target_link_libraries(${MAIN_TARGET} m mpi) - DEAL_II_SETUP_TARGET(${MAIN_TARGET}) + if(NOT DOC_ONLY) + DEAL_II_SETUP_TARGET(${MAIN_TARGET}) + endif() unset(MAIN_TARGET) unset(ParameterLib) unset(TestsLib) endforeach() - -set(SMALL_TEST testrun) -set(TARGET ${SMALL_TEST}) -set(SMALL_TEST_SRC - temporary_test.cpp) -add_executable(${SMALL_TEST} ${SMALL_TEST_SRC}) -target_link_libraries(${SMALL_TEST} mpi ) -DEAL_II_SETUP_TARGET(${SMALL_TEST}) -message("Adding executable " ${SMALL_TEST} " with files " ${SMALL_TEST_SRC} "\n") -unset(SMALL_TEST) -unset(TARGET ${SMALL_TEST}) - -#set(SMALL_TEST instantiate_vector_ad) -#set(TARGET ${SMALL_TEST}) -#set(SMALL_TEST_SRC -# instantiate_vector_ad.cpp) -#add_executable(${SMALL_TEST} ${SMALL_TEST_SRC}) -#target_link_libraries(${SMALL_TEST} mpi ) -#DEAL_II_SETUP_TARGET(${SMALL_TEST}) -#message("Adding executable " ${SMALL_TEST} " with files " ${SMALL_TEST_SRC} "\n") -#unset(SMALL_TEST) -#unset(TARGET ${SMALL_TEST}) - -set(SMALL_TEST manifoldrun) -set(TARGET ${SMALL_TEST}) -set(SMALL_TEST_SRC - test_manifold.cpp) - -add_executable(${SMALL_TEST} ${SMALL_TEST_SRC}) -# target_link_libraries(${SMALL_TEST} m mpi) - -DEAL_II_SETUP_TARGET(${SMALL_TEST}) -message("Adding executable " ${SMALL_TEST} " with files " ${SMALL_TEST_SRC} "\n") - -unset(SMALL_TEST) -unset(TARGET ${SMALL_TEST}) diff --git a/tests/unit_tests/grid/mapping_manifold.cpp b/src/deprecated/mapping_manifold.cpp similarity index 100% rename from tests/unit_tests/grid/mapping_manifold.cpp rename to src/deprecated/mapping_manifold.cpp diff --git a/src/dg/CMakeLists.txt b/src/dg/CMakeLists.txt index 0ceb8abf2..e7324e16d 100644 --- a/src/dg/CMakeLists.txt +++ b/src/dg/CMakeLists.txt @@ -1,6 +1,7 @@ set(DG_SOURCE high_order_grid.cpp dg.cpp + drdx.cpp weak_dg.cpp strong_dg.cpp ) @@ -20,7 +21,9 @@ foreach(dim RANGE 1 3) target_link_libraries(${DiscontinuousGalerkinLib} ${NumericalFluxLib}) target_link_libraries(${DiscontinuousGalerkinLib} ${PhysicsLib}) # Setup target with deal.II - DEAL_II_SETUP_TARGET(${DiscontinuousGalerkinLib}) + if(NOT DOC_ONLY) + DEAL_II_SETUP_TARGET(${DiscontinuousGalerkinLib}) + endif() unset(DiscontinuousGalerkinLib) unset(NumericalFluxLib) diff --git a/src/dg/dg.cpp b/src/dg/dg.cpp index fcc69a39e..91e25544b 100644 --- a/src/dg/dg.cpp +++ b/src/dg/dg.cpp @@ -295,398 +295,262 @@ DGBase::~DGBase () // //set_all_cells_fe_degree ( max_degree-min_degree); // } - template -void DGBase::assemble_residual (const bool compute_dRdW) +template +real DGBase::evaluate_penalty_scaling ( + const DoFCellAccessorType &cell, + const int iface, + const dealii::hp::FECollection fe_collection) const { - right_hand_side = 0; - if (compute_dRdW) system_matrix = 0; + const unsigned int fe_index = cell->active_fe_index(); + const unsigned int degree = fe_collection[fe_index].tensor_degree(); + const unsigned int degsq = (degree == 0) ? 1 : degree * (degree+1); - // For now assume same polynomial degree across domain - const unsigned int max_dofs_per_cell = dof_handler.get_fe_collection().max_dofs_per_cell(); - std::vector current_dofs_indices(max_dofs_per_cell); - std::vector neighbor_dofs_indices(max_dofs_per_cell); + const unsigned int normal_direction = dealii::GeometryInfo::unit_normal_direction[iface]; + const real vol_div_facearea = cell->extent_in_direction(normal_direction); - //dealii::hp::MappingCollection mapping_collection(*(high_order_grid.mapping_fe_field)); - //const dealii::MappingManifold mapping; - //const dealii::MappingQ mapping(max_degree+1); - const auto mapping = (*(high_order_grid.mapping_fe_field)); - dealii::hp::MappingCollection mapping_collection(mapping); - - dealii::hp::FEValues fe_values_collection_volume (mapping_collection, fe_collection, volume_quadrature_collection, this->volume_update_flags); ///< FEValues of volume. - dealii::hp::FEFaceValues fe_values_collection_face_int (mapping_collection, fe_collection, face_quadrature_collection, this->face_update_flags); ///< FEValues of interior face. - dealii::hp::FEFaceValues fe_values_collection_face_ext (mapping_collection, fe_collection, face_quadrature_collection, this->neighbor_face_update_flags); ///< FEValues of exterior face. - dealii::hp::FESubfaceValues fe_values_collection_subface (mapping_collection, fe_collection, face_quadrature_collection, this->face_update_flags); ///< FEValues of subface. - - dealii::hp::FEValues fe_values_collection_volume_lagrange (mapping_collection, fe_collection_lagrange, volume_quadrature_collection, this->volume_update_flags); - - unsigned int n_cell_visited = 0; - unsigned int n_face_visited = 0; - - solution.update_ghost_values(); - for (auto current_cell = dof_handler.begin_active(); current_cell != dof_handler.end(); ++current_cell) { - if (!current_cell->is_locally_owned()) continue; - n_cell_visited++; - - // Current reference element related to this physical cell - const unsigned int mapping_index = 0; - const unsigned int fe_index_curr_cell = current_cell->active_fe_index(); - const unsigned int quad_index = fe_index_curr_cell; - const dealii::FESystem ¤t_fe_ref = fe_collection[fe_index_curr_cell]; - const unsigned int curr_cell_degree = current_fe_ref.tensor_degree(); - const unsigned int n_dofs_curr_cell = current_fe_ref.n_dofs_per_cell(); - - // Local vector contribution from each cell - dealii::Vector current_cell_rhs (n_dofs_curr_cell); // Defaults to 0.0 initialization - - // Obtain the mapping from local dof indices to global dof indices - current_dofs_indices.resize(n_dofs_curr_cell); - current_cell->get_dof_indices (current_dofs_indices); - - // fe_values_collection.reinit(current_cell, quad_collection_index, mapping_collection_index, fe_collection_index) - fe_values_collection_volume.reinit (current_cell, quad_index, mapping_index, fe_index_curr_cell); - const dealii::FEValues &fe_values_volume = fe_values_collection_volume.get_present_fe_values(); + const real penalty = degsq / vol_div_facearea; + return penalty; +} - dealii::TriaIterator> cell_iterator = static_cast> > (current_cell); - //if (!(all_parameters->use_weak_form)) fe_values_collection_volume_lagrange.reinit (current_cell, quad_index, mapping_index, fe_index_curr_cell); - fe_values_collection_volume_lagrange.reinit (cell_iterator, quad_index, mapping_index, fe_index_curr_cell); - const dealii::FEValues &fe_values_lagrange = fe_values_collection_volume_lagrange.get_present_fe_values(); - if ( compute_dRdW ) { - assemble_volume_terms_implicit (fe_values_volume, current_dofs_indices, current_cell_rhs, fe_values_lagrange); - } else { - assemble_volume_terms_explicit (fe_values_volume, current_dofs_indices, current_cell_rhs, fe_values_lagrange); +template +template +bool DGBase::current_cell_should_do_the_work (const DoFCellAccessorType1 ¤t_cell, const DoFCellAccessorType2 &neighbor_cell) const +{ + if (neighbor_cell->is_ghost()) { + // In the case the neighbor is a ghost cell, we let the processor with the lower rank do the work on that face + // We cannot use the cell->index() because the index is relative to the distributed triangulation + // Therefore, the cell index of a ghost cell might be different to the physical cell index even if they refer to the same cell + return (current_cell->subdomain_id() < neighbor_cell->subdomain_id()); + } else { + // Locally owned neighbor cell + Assert(neighbor_cell->is_locally_owned(), dealii::ExcMessage("If not ghost, neighbor should be locally owned.")); + + if (current_cell->index() < neighbor_cell->index()) { + // Cell with lower index does work + return true; + } else if (neighbor_cell->index() == current_cell->index()) { + // If both cells have same index + // See https://www.dealii.org/developer/doxygen/deal.II/classTriaAccessorBase.html#a695efcbe84fefef3e4c93ee7bdb446ad + // then cell at the lower level does the work + return (current_cell->level() < neighbor_cell->level()); } + return false; + } + Assert(0==1, dealii::ExcMessage("Should not have reached here. Somehow another possible case has not been considered when two cells have the same coarseness.")); + return false; +} - for (unsigned int iface=0; iface < dealii::GeometryInfo::faces_per_cell; ++iface) { +template +template +void DGBase::assemble_cell_residual ( + const DoFCellAccessorType ¤t_cell, + const bool compute_dRdW, + dealii::hp::FEValues &fe_values_collection_volume, + dealii::hp::FEFaceValues &fe_values_collection_face_int, + dealii::hp::FEFaceValues &fe_values_collection_face_ext, + dealii::hp::FESubfaceValues &fe_values_collection_subface, + dealii::hp::FEValues &fe_values_collection_volume_lagrange, + dealii::LinearAlgebra::distributed::Vector &rhs) +{ + std::vector current_dofs_indices; + std::vector neighbor_dofs_indices; - auto current_face = current_cell->face(iface); - auto neighbor_cell = current_cell->neighbor(iface); + // Current reference element related to this physical cell + const int i_fele = current_cell->active_fe_index(); + const int i_quad = i_fele; + const int i_mapp = 0; - // See tutorial step-30 for breakdown of 4 face cases + const dealii::FESystem ¤t_fe_ref = fe_collection[i_fele]; + const unsigned int n_dofs_curr_cell = current_fe_ref.n_dofs_per_cell(); - // Case 1: - // Face at boundary - if (current_face->at_boundary() && !current_cell->has_periodic_neighbor(iface) ) { + // Local vector contribution from each cell + dealii::Vector current_cell_rhs (n_dofs_curr_cell); // Defaults to 0.0 initialization - n_face_visited++; + // Obtain the mapping from local dof indices to global dof indices + current_dofs_indices.resize(n_dofs_curr_cell); + current_cell->get_dof_indices (current_dofs_indices); - fe_values_collection_face_int.reinit (current_cell, iface, quad_index, mapping_index, fe_index_curr_cell); + fe_values_collection_volume.reinit (current_cell, i_quad, i_mapp, i_fele); + const dealii::FEValues &fe_values_volume = fe_values_collection_volume.get_present_fe_values(); - if(current_face->at_boundary() && all_parameters->use_periodic_bc == true && dim == 1) //using periodic BCs (for 1d) - { - int cell_index = current_cell->index(); - //int cell_index = current_cell->index(); - if (cell_index == 0 && iface == 0) - { - fe_values_collection_face_int.reinit(current_cell, iface, quad_index, mapping_index, fe_index_curr_cell); - neighbor_cell = dof_handler.begin_active(); - for (unsigned int i = 0 ; i < triangulation->n_active_cells() - 1; ++i) - { - ++neighbor_cell; - } - neighbor_cell->get_dof_indices(neighbor_dofs_indices); - const unsigned int fe_index_neigh_cell = neighbor_cell->active_fe_index(); - const unsigned int quad_index_neigh_cell = fe_index_neigh_cell; - const unsigned int mapping_index_neigh_cell = 0; + dealii::TriaIterator> cell_iterator = static_cast> > (current_cell); + //if (!(all_parameters->use_weak_form)) fe_values_collection_volume_lagrange.reinit (current_cell, i_quad, i_mapp, i_fele); + fe_values_collection_volume_lagrange.reinit (cell_iterator, i_quad, i_mapp, i_fele); + const dealii::FEValues &fe_values_lagrange = fe_values_collection_volume_lagrange.get_present_fe_values(); - fe_values_collection_face_ext.reinit(neighbor_cell,(iface == 1) ? 0 : 1,quad_index_neigh_cell,mapping_index_neigh_cell,fe_index_neigh_cell); + if ( compute_dRdW ) { + assemble_volume_terms_implicit (fe_values_volume, current_dofs_indices, current_cell_rhs, fe_values_lagrange); + } else { + assemble_volume_terms_explicit (fe_values_volume, current_dofs_indices, current_cell_rhs, fe_values_lagrange); + } - } - else if (cell_index == (int) triangulation->n_active_cells() - 1 && iface == 1) - { - fe_values_collection_face_int.reinit(current_cell, iface, quad_index, mapping_index, fe_index_curr_cell); - neighbor_cell = dof_handler.begin_active(); - neighbor_cell->get_dof_indices(neighbor_dofs_indices); - const unsigned int fe_index_neigh_cell = neighbor_cell->active_fe_index(); - const unsigned int quad_index_neigh_cell = fe_index_neigh_cell; - const unsigned int mapping_index_neigh_cell = 0; - fe_values_collection_face_ext.reinit(neighbor_cell,(iface == 1) ? 0 : 1, quad_index_neigh_cell, mapping_index_neigh_cell, fe_index_neigh_cell); //not sure how changing the face number would work in dim!=1-dimensions. - } + for (unsigned int iface=0; iface < dealii::GeometryInfo::faces_per_cell; ++iface) { - //std::cout << "cell " << current_cell->index() << "'s " << iface << "th face has neighbour: " << neighbor_cell->index() << std::endl; - const int neighbor_face_no = (iface ==1) ? 0:1; - const unsigned int fe_index_neigh_cell = neighbor_cell->active_fe_index(); + auto current_face = current_cell->face(iface); + auto neighbor_cell = current_cell->neighbor(iface); - const dealii::FEFaceValues &fe_values_face_int = fe_values_collection_face_int.get_present_fe_values(); - const dealii::FEFaceValues &fe_values_face_ext = fe_values_collection_face_ext.get_present_fe_values(); + // Case 1: Face at boundary + if (current_face->at_boundary() && !current_cell->has_periodic_neighbor(iface) ) { - const dealii::FESystem &neigh_fe_ref = fe_collection[fe_index_neigh_cell]; - const unsigned int neigh_cell_degree = neigh_fe_ref.tensor_degree(); - const unsigned int n_dofs_neigh_cell = neigh_fe_ref.n_dofs_per_cell(); + fe_values_collection_face_int.reinit(current_cell, iface, i_quad, i_mapp, i_fele); - dealii::Vector neighbor_cell_rhs (n_dofs_neigh_cell); // Defaults to 0.0 initialization + // Case 1.1: 1D Periodic boundary condition + if(current_face->at_boundary() && all_parameters->use_periodic_bc == true && dim == 1) { + const int neighbor_iface = (iface == 1) ? 0 : 1; - const unsigned int normal_direction1 = dealii::GeometryInfo::unit_normal_direction[iface]; - const unsigned int normal_direction2 = dealii::GeometryInfo::unit_normal_direction[neighbor_face_no]; - const unsigned int deg1sq = (curr_cell_degree == 0) ? 1 : curr_cell_degree * (curr_cell_degree+1); - const unsigned int deg2sq = (neigh_cell_degree == 0) ? 1 : neigh_cell_degree * (neigh_cell_degree+1); + int cell_index = current_cell->index(); + if (cell_index == 0 && iface == 0) { + // First cell of the domain, neighbor is the last. + neighbor_cell = dof_handler.begin_active(); + for (unsigned int i = 0 ; i < triangulation->n_active_cells() - 1; ++i) { + ++neighbor_cell; + } + } else if (cell_index == (int) triangulation->n_active_cells() - 1 && iface == 1) { + // Last cell of the domain, neighbor is the first. + neighbor_cell = dof_handler.begin_active(); + } - //const real vol_div_facearea1 = current_cell->extent_in_direction(normal_direction1) / current_face->number_of_children(); - const real vol_div_facearea1 = current_cell->extent_in_direction(normal_direction1); - const real vol_div_facearea2 = neighbor_cell->extent_in_direction(normal_direction2); + const int i_fele_n = neighbor_cell->active_fe_index(), i_quad_n = i_fele_n, i_mapp_n = 0; + const unsigned int n_dofs_neigh_cell = fe_collection[i_fele_n].n_dofs_per_cell(); + neighbor_dofs_indices.resize(n_dofs_neigh_cell); + fe_values_collection_face_ext.reinit(neighbor_cell, neighbor_iface, i_quad_n, i_mapp_n, i_fele_n); - const real penalty1 = deg1sq / vol_div_facearea1; - const real penalty2 = deg2sq / vol_div_facearea2; + const dealii::FEFaceValues &fe_values_face_int = fe_values_collection_face_int.get_present_fe_values(); + const dealii::FEFaceValues &fe_values_face_ext = fe_values_collection_face_ext.get_present_fe_values(); - real penalty = 0.5 * ( penalty1 + penalty2 ); + dealii::Vector neighbor_cell_rhs (n_dofs_neigh_cell); // Defaults to 0.0 initialization - if ( compute_dRdW ) { - assemble_face_term_implicit ( - fe_values_face_int, fe_values_face_ext, - penalty, - current_dofs_indices, neighbor_dofs_indices, - current_cell_rhs, neighbor_cell_rhs); - } else { - assemble_face_term_explicit ( - fe_values_face_int, fe_values_face_ext, - penalty, - current_dofs_indices, neighbor_dofs_indices, - current_cell_rhs, neighbor_cell_rhs); - } + const real penalty1 = evaluate_penalty_scaling (current_cell, iface, fe_collection); + const real penalty2 = evaluate_penalty_scaling (neighbor_cell, neighbor_iface, fe_collection); + const real penalty = 0.5 * (penalty1 + penalty2); + if ( compute_dRdW ) { + assemble_face_term_implicit ( + fe_values_face_int, fe_values_face_ext, + penalty, + current_dofs_indices, neighbor_dofs_indices, + current_cell_rhs, neighbor_cell_rhs); } else { - const dealii::FEFaceValues &fe_values_face_int = fe_values_collection_face_int.get_present_fe_values(); - const unsigned int deg1sq = (curr_cell_degree == 0) ? 1 : curr_cell_degree * (curr_cell_degree+1); - const unsigned int normal_direction = dealii::GeometryInfo::unit_normal_direction[iface]; - const real vol_div_facearea1 = current_cell->extent_in_direction(normal_direction); - - real penalty = deg1sq / vol_div_facearea1; - - const unsigned int boundary_id = current_face->boundary_id(); - // Need to somehow get boundary type from the mesh - if ( compute_dRdW ) { - assemble_boundary_term_implicit (boundary_id, fe_values_face_int, penalty, current_dofs_indices, current_cell_rhs); - } else { - assemble_boundary_term_explicit (boundary_id, fe_values_face_int, penalty, current_dofs_indices, current_cell_rhs); - } - } - - //CASE 1.5: periodic boundary conditions - //note that periodicity is not adapted for hp adaptivity yet. this needs to be figured out in the future - } else if (current_face->at_boundary() && current_cell->has_periodic_neighbor(iface)){ - - neighbor_cell = current_cell->periodic_neighbor(iface); - //std::cout << "cell " << current_cell->index() << " at boundary" <index() << std::endl; - - - if (!current_cell->periodic_neighbor_is_coarser(iface) && - (neighbor_cell->index() > current_cell->index() || - (neighbor_cell->index() == current_cell->index() && current_cell->level() < neighbor_cell->level()) - ) - ) - { - n_face_visited++; - Assert (current_cell->periodic_neighbor(iface).state() == dealii::IteratorState::valid, dealii::ExcInternalError()); - - - // Corresponding face of the neighbor. - // e.g. The 4th face of the current cell might correspond to the 3rd face of the neighbor - const unsigned int neighbor_face_no = current_cell->periodic_neighbor_of_periodic_neighbor(iface); - - // Get information about neighbor cell - const unsigned int fe_index_neigh_cell = neighbor_cell->active_fe_index(); - const unsigned int quad_index_neigh_cell = fe_index_neigh_cell; - const unsigned int mapping_index_neigh_cell = 0; - const dealii::FESystem &neigh_fe_ref = fe_collection[fe_index_neigh_cell]; - const unsigned int neigh_cell_degree = neigh_fe_ref.tensor_degree(); - const unsigned int n_dofs_neigh_cell = neigh_fe_ref.n_dofs_per_cell(); - - // Local rhs contribution from neighbor - dealii::Vector neighbor_cell_rhs (n_dofs_neigh_cell); // Defaults to 0.0 initialization - - // Obtain the mapping from local dof indices to global dof indices for neighbor cell - neighbor_dofs_indices.resize(n_dofs_neigh_cell); - neighbor_cell->get_dof_indices (neighbor_dofs_indices); - - fe_values_collection_face_int.reinit (current_cell, iface, quad_index, mapping_index, fe_index_curr_cell); - const dealii::FEFaceValues &fe_values_face_int = fe_values_collection_face_int.get_present_fe_values(); - fe_values_collection_face_ext.reinit (neighbor_cell, neighbor_face_no, quad_index_neigh_cell, mapping_index_neigh_cell, fe_index_neigh_cell); - const dealii::FEFaceValues &fe_values_face_ext = fe_values_collection_face_ext.get_present_fe_values(); - - const unsigned int normal_direction1 = dealii::GeometryInfo::unit_normal_direction[iface]; - const unsigned int normal_direction2 = dealii::GeometryInfo::unit_normal_direction[neighbor_face_no]; - const unsigned int deg1sq = (curr_cell_degree == 0) ? 1 : curr_cell_degree * (curr_cell_degree+1); - const unsigned int deg2sq = (neigh_cell_degree == 0) ? 1 : neigh_cell_degree * (neigh_cell_degree+1); - - //const real vol_div_facearea1 = current_cell->extent_in_direction(normal_direction1) / current_face->number_of_children(); - const real vol_div_facearea1 = current_cell->extent_in_direction(normal_direction1); - const real vol_div_facearea2 = neighbor_cell->extent_in_direction(normal_direction2); - - const real penalty1 = deg1sq / vol_div_facearea1; - const real penalty2 = deg2sq / vol_div_facearea2; - - real penalty = 0.5 * ( penalty1 + penalty2 ); - //penalty = 1;//99; - - if ( compute_dRdW ) { - assemble_face_term_implicit ( - fe_values_face_int, fe_values_face_ext, - penalty, - current_dofs_indices, neighbor_dofs_indices, - current_cell_rhs, neighbor_cell_rhs); - } else { - assemble_face_term_explicit ( - fe_values_face_int, fe_values_face_ext, - penalty, - current_dofs_indices, neighbor_dofs_indices, - current_cell_rhs, neighbor_cell_rhs); - } - - // Add local contribution from neighbor cell to global vector - for (unsigned int i=0; i &fe_values_face_int = fe_values_collection_face_int.get_present_fe_values(); - // Case 2: - // Neighbour is finer occurs if the face has children - // In this case, we loop over the current large face's subfaces and visit multiple neighbors - } else if (current_face->has_children()) { - - Assert (current_cell->neighbor(iface).state() == dealii::IteratorState::valid, dealii::ExcInternalError()); + const real penalty = evaluate_penalty_scaling (current_cell, iface, fe_collection); - // Obtain cell neighbour - const unsigned int neighbor_face_no = current_cell->neighbor_face_no(iface); + const unsigned int boundary_id = current_face->boundary_id(); + // Need to somehow get boundary type from the mesh + if ( compute_dRdW ) { + assemble_boundary_term_implicit (boundary_id, fe_values_face_int, penalty, current_dofs_indices, current_cell_rhs); + } else { + assemble_boundary_term_explicit (boundary_id, fe_values_face_int, penalty, current_dofs_indices, current_cell_rhs); + } + } - for (unsigned int subface_no=0; subface_no < current_face->number_of_children(); ++subface_no) { + //CASE 1.5: periodic boundary conditions + //note that periodicity is not adapted for hp adaptivity yet. this needs to be figured out in the future + } else if (current_face->at_boundary() && current_cell->has_periodic_neighbor(iface)){ - n_face_visited++; + neighbor_cell = current_cell->periodic_neighbor(iface); + //std::cout << "cell " << current_cell->index() << " at boundary" <index() << std::endl; - // Get neighbor on ith subface - auto neighbor_cell = current_cell->neighbor_child_on_subface (iface, subface_no); - // Since the neighbor cell is finer than the current cell, it should not have more children - Assert (!neighbor_cell->has_children(), dealii::ExcInternalError()); - // Get information about neighbor cell - const unsigned int fe_index_neigh_cell = neighbor_cell->active_fe_index(); - const unsigned int quad_index_neigh_cell = fe_index_neigh_cell; - const unsigned int mapping_index_neigh_cell = 0; - const dealii::FESystem &neigh_fe_ref = fe_collection[fe_index_neigh_cell]; - const unsigned int neigh_cell_degree = neigh_fe_ref.tensor_degree(); - const unsigned int n_dofs_neigh_cell = neigh_fe_ref.n_dofs_per_cell(); + if (!current_cell->periodic_neighbor_is_coarser(iface) && current_cell_should_do_the_work(current_cell, neighbor_cell)) { - dealii::Vector neighbor_cell_rhs (n_dofs_neigh_cell); // Defaults to 0.0 initialization + Assert (current_cell->periodic_neighbor(iface).state() == dealii::IteratorState::valid, dealii::ExcInternalError()); - // Obtain the mapping from local dof indices to global dof indices for neighbor cell - neighbor_dofs_indices.resize(n_dofs_neigh_cell); - neighbor_cell->get_dof_indices (neighbor_dofs_indices); + const unsigned int n_dofs_neigh_cell = fe_collection[neighbor_cell->active_fe_index()].n_dofs_per_cell(); + dealii::Vector neighbor_cell_rhs (n_dofs_neigh_cell); // Defaults to 0.0 initialization - fe_values_collection_subface.reinit (current_cell, iface, subface_no, quad_index, mapping_index, fe_index_curr_cell); - const dealii::FESubfaceValues &fe_values_face_int = fe_values_collection_subface.get_present_fe_values(); + // Obtain the mapping from local dof indices to global dof indices for neighbor cell + neighbor_dofs_indices.resize(n_dofs_neigh_cell); + neighbor_cell->get_dof_indices (neighbor_dofs_indices); - fe_values_collection_face_ext.reinit (neighbor_cell, neighbor_face_no, quad_index_neigh_cell, mapping_index_neigh_cell, fe_index_neigh_cell); - const dealii::FEFaceValues &fe_values_face_ext = fe_values_collection_face_ext.get_present_fe_values(); + fe_values_collection_face_int.reinit (current_cell, iface, i_quad, i_mapp, i_fele); + const dealii::FEFaceValues &fe_values_face_int = fe_values_collection_face_int.get_present_fe_values(); - const unsigned int normal_direction1 = dealii::GeometryInfo::unit_normal_direction[iface]; - const unsigned int normal_direction2 = dealii::GeometryInfo::unit_normal_direction[neighbor_face_no]; - const unsigned int deg1sq = (curr_cell_degree == 0) ? 1 : curr_cell_degree * (curr_cell_degree+1); - const unsigned int deg2sq = (neigh_cell_degree == 0) ? 1 : neigh_cell_degree * (neigh_cell_degree+1); + // Corresponding face of the neighbor. + const unsigned int neighbor_iface = current_cell->periodic_neighbor_of_periodic_neighbor(iface); - const real vol_div_facearea1 = current_cell->extent_in_direction(normal_direction1); - const real vol_div_facearea2 = neighbor_cell->extent_in_direction(normal_direction2); + const int i_fele_n = neighbor_cell->active_fe_index(), i_quad_n = i_fele_n, i_mapp_n = 0; + fe_values_collection_face_ext.reinit (neighbor_cell, neighbor_iface, i_quad_n, i_mapp_n, i_fele_n); + const dealii::FEFaceValues &fe_values_face_ext = fe_values_collection_face_ext.get_present_fe_values(); - const real penalty1 = deg1sq / vol_div_facearea1; - const real penalty2 = deg2sq / vol_div_facearea2; - - real penalty = 0.5 * ( penalty1 + penalty2 ); + const real penalty1 = evaluate_penalty_scaling (current_cell, iface, fe_collection); + const real penalty2 = evaluate_penalty_scaling (neighbor_cell, neighbor_iface, fe_collection); + const real penalty = 0.5 * (penalty1 + penalty2); - if ( compute_dRdW ) { - assemble_face_term_implicit ( - fe_values_face_int, fe_values_face_ext, - penalty, - current_dofs_indices, neighbor_dofs_indices, - current_cell_rhs, neighbor_cell_rhs); - } else { - assemble_face_term_explicit ( + if ( compute_dRdW ) { + assemble_face_term_implicit ( + fe_values_face_int, fe_values_face_ext, + penalty, + current_dofs_indices, neighbor_dofs_indices, + current_cell_rhs, neighbor_cell_rhs); + } else { + assemble_face_term_explicit ( fe_values_face_int, fe_values_face_ext, penalty, current_dofs_indices, neighbor_dofs_indices, current_cell_rhs, neighbor_cell_rhs); - } - // Add local contribution from neighbor cell to global vector - for (unsigned int i=0; ineighbor_is_coarser(iface)) - // In the case the neighbor is a ghost cell, we let the processor with the lower rank do the work on that face - // We cannot use the cell->index() because the index is relative to the distributed triangulation - // Therefore, the cell index of a ghost cell might be different to the physical cell index even if they refer to the same cell - && neighbor_cell->is_ghost() - && current_cell->subdomain_id() < neighbor_cell->subdomain_id() - ) - || - ( !(current_cell->neighbor_is_coarser(iface)) - // In the case the neighbor is a local cell, we let the cell with the lower index do the work on that face - && neighbor_cell->is_locally_owned() - && - ( // Cell with lower index does work - current_cell->index() < neighbor_cell->index() - || - // If both cells have same index - // See https://www.dealii.org/developer/doxygen/deal.II/classTriaAccessorBase.html#a695efcbe84fefef3e4c93ee7bdb446ad - // then cell at the lower level does the work - (neighbor_cell->index() == current_cell->index() && current_cell->level() < neighbor_cell->level()) - ) - ) - ) - { - n_face_visited++; - Assert (current_cell->neighbor(iface).state() == dealii::IteratorState::valid, dealii::ExcInternalError()); + // Add local contribution from neighbor cell to global vector + for (unsigned int i=0; ineighbor_or_periodic_neighbor(iface); - // Corresponding face of the neighbor. - // e.g. The 4th face of the current cell might correspond to the 3rd face of the neighbor - const unsigned int neighbor_face_no = current_cell->neighbor_of_neighbor(iface); - // Get information about neighbor cell - const unsigned int fe_index_neigh_cell = neighbor_cell->active_fe_index(); - const unsigned int quad_index_neigh_cell = fe_index_neigh_cell; - const unsigned int mapping_index_neigh_cell = 0; - const dealii::FESystem &neigh_fe_ref = fe_collection[fe_index_neigh_cell]; - const unsigned int neigh_cell_degree = neigh_fe_ref.tensor_degree(); - const unsigned int n_dofs_neigh_cell = neigh_fe_ref.n_dofs_per_cell(); + // Case 2: + // Neighbour is finer occurs if the face has children + // In this case, we loop over the current large face's subfaces and visit multiple neighbors + } else if (current_face->has_children()) { + + Assert (current_cell->neighbor(iface).state() == dealii::IteratorState::valid, dealii::ExcInternalError()); + Assert (current_cell->neighbor(iface)->has_children(), dealii::ExcInternalError()); + + // Obtain cell neighbour + const unsigned int neighbor_iface = current_cell->neighbor_face_no(iface); + + for (unsigned int subface_no=0; subface_no < current_face->number_of_children(); ++subface_no) { + + // Get neighbor on ith subface + auto neighbor_cell = current_cell->neighbor_child_on_subface (iface, subface_no); + // Since the neighbor cell is finer than the current cell, it should not have more children + Assert (!neighbor_cell->has_children(), dealii::ExcInternalError()); + Assert (neighbor_cell->neighbor(neighbor_iface) == current_cell, dealii::ExcInternalError()); - // Local rhs contribution from neighbor - dealii::Vector neighbor_cell_rhs (n_dofs_neigh_cell); // Defaults to 0.0 initialization + const int i_fele_n = neighbor_cell->active_fe_index(), i_quad_n = i_fele_n, i_mapp_n = 0; + + const unsigned int n_dofs_neigh_cell = fe_collection[i_fele_n].n_dofs_per_cell(); + dealii::Vector neighbor_cell_rhs (n_dofs_neigh_cell); // Defaults to 0.0 initialization // Obtain the mapping from local dof indices to global dof indices for neighbor cell neighbor_dofs_indices.resize(n_dofs_neigh_cell); neighbor_cell->get_dof_indices (neighbor_dofs_indices); - fe_values_collection_face_int.reinit (current_cell, iface, quad_index, mapping_index, fe_index_curr_cell); - const dealii::FEFaceValues &fe_values_face_int = fe_values_collection_face_int.get_present_fe_values(); - fe_values_collection_face_ext.reinit (neighbor_cell, neighbor_face_no, quad_index_neigh_cell, mapping_index_neigh_cell, fe_index_neigh_cell); - const dealii::FEFaceValues &fe_values_face_ext = fe_values_collection_face_ext.get_present_fe_values(); + fe_values_collection_subface.reinit (current_cell, iface, subface_no, i_quad, i_mapp, i_fele); + const dealii::FESubfaceValues &fe_values_face_int = fe_values_collection_subface.get_present_fe_values(); - const unsigned int normal_direction1 = dealii::GeometryInfo::unit_normal_direction[iface]; - const unsigned int normal_direction2 = dealii::GeometryInfo::unit_normal_direction[neighbor_face_no]; - const unsigned int deg1sq = (curr_cell_degree == 0) ? 1 : curr_cell_degree * (curr_cell_degree+1); - const unsigned int deg2sq = (neigh_cell_degree == 0) ? 1 : neigh_cell_degree * (neigh_cell_degree+1); - - //const real vol_div_facearea1 = current_cell->extent_in_direction(normal_direction1) / current_face->number_of_children(); - const real vol_div_facearea1 = current_cell->extent_in_direction(normal_direction1); - const real vol_div_facearea2 = neighbor_cell->extent_in_direction(normal_direction2); + fe_values_collection_face_ext.reinit (neighbor_cell, neighbor_iface, i_quad_n, i_mapp_n, i_fele_n); + const dealii::FEFaceValues &fe_values_face_ext = fe_values_collection_face_ext.get_present_fe_values(); - const real penalty1 = deg1sq / vol_div_facearea1; - const real penalty2 = deg2sq / vol_div_facearea2; - - real penalty = 0.5 * ( penalty1 + penalty2 ); - //penalty = 1;//99; + const real penalty1 = evaluate_penalty_scaling (current_cell, iface, fe_collection); + const real penalty2 = evaluate_penalty_scaling (neighbor_cell, neighbor_iface, fe_collection); + const real penalty = 0.5 * (penalty1 + penalty2); if ( compute_dRdW ) { assemble_face_term_implicit ( @@ -696,29 +560,117 @@ void DGBase::assemble_residual (const bool compute_dRdW) current_cell_rhs, neighbor_cell_rhs); } else { assemble_face_term_explicit ( - fe_values_face_int, fe_values_face_ext, - penalty, - current_dofs_indices, neighbor_dofs_indices, - current_cell_rhs, neighbor_cell_rhs); + fe_values_face_int, fe_values_face_ext, + penalty, + current_dofs_indices, neighbor_dofs_indices, + current_cell_rhs, neighbor_cell_rhs); } - // Add local contribution from neighbor cell to global vector for (unsigned int i=0; ineighbor_is_coarser(iface)) && current_cell_should_do_the_work(current_cell, neighbor_cell) ) { + Assert (current_cell->neighbor(iface).state() == dealii::IteratorState::valid, dealii::ExcInternalError()); + + auto neighbor_cell = current_cell->neighbor_or_periodic_neighbor(iface); + // Corresponding face of the neighbor. + // e.g. The 4th face of the current cell might correspond to the 3rd face of the neighbor + const unsigned int neighbor_iface = current_cell->neighbor_of_neighbor(iface); + + // Get information about neighbor cell + const unsigned int n_dofs_neigh_cell = fe_collection[neighbor_cell->active_fe_index()].n_dofs_per_cell(); + + // Local rhs contribution from neighbor + dealii::Vector neighbor_cell_rhs (n_dofs_neigh_cell); // Defaults to 0.0 initialization + + // Obtain the mapping from local dof indices to global dof indices for neighbor cell + neighbor_dofs_indices.resize(n_dofs_neigh_cell); + neighbor_cell->get_dof_indices (neighbor_dofs_indices); + + fe_values_collection_face_int.reinit (current_cell, iface, i_quad, i_mapp, i_fele); + const dealii::FEFaceValues &fe_values_face_int = fe_values_collection_face_int.get_present_fe_values(); + + const int i_fele_n = neighbor_cell->active_fe_index(), i_quad_n = i_fele_n, i_mapp_n = 0; + fe_values_collection_face_ext.reinit (neighbor_cell, neighbor_iface, i_quad_n, i_mapp_n, i_fele_n); + const dealii::FEFaceValues &fe_values_face_ext = fe_values_collection_face_ext.get_present_fe_values(); + + const real penalty1 = evaluate_penalty_scaling (current_cell, iface, fe_collection); + const real penalty2 = evaluate_penalty_scaling (neighbor_cell, neighbor_iface, fe_collection); + const real penalty = 0.5 * (penalty1 + penalty2); - for (unsigned int i=0; i +void DGBase::assemble_residual (const bool compute_dRdW) +{ + right_hand_side = 0; + + if (compute_dRdW) system_matrix = 0; + + //dealii::hp::MappingCollection mapping_collection(*(high_order_grid.mapping_fe_field)); + //const dealii::MappingManifold mapping; + //const dealii::MappingQ mapping(max_degree+1); + const auto mapping = (*(high_order_grid.mapping_fe_field)); + dealii::hp::MappingCollection mapping_collection(mapping); + + dealii::hp::FEValues fe_values_collection_volume (mapping_collection, fe_collection, volume_quadrature_collection, this->volume_update_flags); ///< FEValues of volume. + dealii::hp::FEFaceValues fe_values_collection_face_int (mapping_collection, fe_collection, face_quadrature_collection, this->face_update_flags); ///< FEValues of interior face. + dealii::hp::FEFaceValues fe_values_collection_face_ext (mapping_collection, fe_collection, face_quadrature_collection, this->neighbor_face_update_flags); ///< FEValues of exterior face. + dealii::hp::FESubfaceValues fe_values_collection_subface (mapping_collection, fe_collection, face_quadrature_collection, this->face_update_flags); ///< FEValues of subface. + + dealii::hp::FEValues fe_values_collection_volume_lagrange (mapping_collection, fe_collection_lagrange, volume_quadrature_collection, this->volume_update_flags); + + solution.update_ghost_values(); + + for (auto current_cell = dof_handler.begin_active(); current_cell != dof_handler.end(); ++current_cell) { + if (!current_cell->is_locally_owned()) continue; + + // Add right-hand side contributions this cell can compute + assemble_cell_residual (current_cell, compute_dRdW, + fe_values_collection_volume, + fe_values_collection_face_int, + fe_values_collection_face_ext, + fe_values_collection_subface, + fe_values_collection_volume_lagrange, + right_hand_side); } // end of cell loop + right_hand_side.compress(dealii::VectorOperation::add); if ( compute_dRdW ) system_matrix.compress(dealii::VectorOperation::add); @@ -789,9 +741,9 @@ void DGBase::output_results_vtk (const unsigned int cycle)// const const int iproc = dealii::Utilities::MPI::this_mpi_process(mpi_communicator); //data_out.build_patches (mapping_collection[mapping_collection.size()-1]); - data_out.build_patches(*(high_order_grid.mapping_fe_field), max_degree, dealii::DataOut>::CurvedCellRegion::curved_inner_cells); + data_out.build_patches(*(high_order_grid.mapping_fe_field), max_degree, dealii::DataOut>::CurvedCellRegion::no_curved_cells); //data_out.build_patches(*(high_order_grid.mapping_fe_field), fe_collection.size(), dealii::DataOut::CurvedCellRegion::curved_inner_cells); - std::string filename = "solution-" + dealii::Utilities::int_to_string(dim, 1) +"D-"; + std::string filename = "solution-" + dealii::Utilities::int_to_string(dim, 1) +"D_maxpoly"+dealii::Utilities::int_to_string(max_degree, 2)+"-"; filename += dealii::Utilities::int_to_string(cycle, 4) + "."; filename += dealii::Utilities::int_to_string(iproc, 4); filename += ".vtu"; @@ -801,13 +753,13 @@ void DGBase::output_results_vtk (const unsigned int cycle)// const if (iproc == 0) { std::vector filenames; for (unsigned int iproc = 0; iproc < dealii::Utilities::MPI::n_mpi_processes(mpi_communicator); ++iproc) { - std::string fn = "solution-" + dealii::Utilities::int_to_string(dim, 1) +"D-"; + std::string fn = "solution-" + dealii::Utilities::int_to_string(dim, 1) +"D_maxpoly"+dealii::Utilities::int_to_string(max_degree, 2)+"-"; fn += dealii::Utilities::int_to_string(cycle, 4) + "."; fn += dealii::Utilities::int_to_string(iproc, 4); fn += ".vtu"; filenames.push_back(fn); } - std::string master_fn = "solution-" + dealii::Utilities::int_to_string(dim, 1) +"D-"; + std::string master_fn = "solution-" + dealii::Utilities::int_to_string(dim, 1) +"D_maxpoly"+dealii::Utilities::int_to_string(max_degree, 2)+"-"; master_fn += dealii::Utilities::int_to_string(cycle, 4) + ".pvtu"; std::ofstream master_output(master_fn); data_out.write_pvtu_record(master_output, filenames); diff --git a/src/dg/dg.h b/src/dg/dg.h index aeecb0719..528cafa49 100644 --- a/src/dg/dg.h +++ b/src/dg/dg.h @@ -67,8 +67,16 @@ template class DGBase { #if PHILIP_DIM==1 // dealii::parallel::distributed::Triangulation does not work for 1D + /** Triangulation to store the grid. + * In 1D, dealii::Triangulation is used. + * In 2D, 3D, dealii::parallel::distributed::Triangulation is used. + */ using Triangulation = dealii::Triangulation; #else + /** Triangulation to store the grid. + * In 1D, dealii::Triangulation is used. + * In 2D, 3D, dealii::parallel::distributed::Triangulation is used. + */ using Triangulation = dealii::parallel::distributed::Triangulation; #endif public: @@ -217,6 +225,16 @@ class DGBase */ dealii::LinearAlgebra::distributed::Vector solution; + /// Evaluate SparsityPattern of dRdX + /* Where R represents the residual and X represents the grid degrees of freedom stored as high_order_grid.nodes. + */ + dealii::SparsityPattern get_dRdX_sparsity_pattern (); + + /// Evaluate dRdX using finite-differences + /* Where R represents the residual and X represents the grid degrees of freedom stored as high_order_grid.nodes. + */ + dealii::TrilinosWrappers::SparseMatrix get_dRdX_finite_differences (dealii::SparsityPattern dRdX_sparsity_pattern); + void initialize_manufactured_solution (); ///< Virtual function defined in DG void output_results_vtk (const unsigned int ith_grid); ///< Output solution @@ -256,6 +274,22 @@ class DGBase //void assemble_residual_dRdW (); void assemble_residual (const bool compute_dRdW=false); + /// Used in assemble_residual(). + /** IMPORTANT: This does not fully compute the cell residual since it might not + * perform the work on all the faces. + * All the active cells must be traversed to ensure that the right hand side is correct. + */ + template + void assemble_cell_residual ( + const DoFCellAccessorType ¤t_cell, + const bool compute_dRdW, + dealii::hp::FEValues &fe_values_collection_volume, + dealii::hp::FEFaceValues &fe_values_collection_face_int, + dealii::hp::FEFaceValues &fe_values_collection_face_ext, + dealii::hp::FESubfaceValues &fe_values_collection_subface, + dealii::hp::FEValues &fe_values_collection_volume_lagrange, + dealii::LinearAlgebra::distributed::Vector &rhs); + /// Finite Element Collection for p-finite-element to represent the solution /** This is a collection of FESystems */ const dealii::hp::FECollection fe_collection; @@ -268,8 +302,11 @@ class DGBase //const dealii::hp::FECollection fe_collection_grid; //const dealii::FESystem fe_grid; + /// Quadrature used to evaluate volume integrals. dealii::hp::QCollection volume_quadrature_collection; + /// Quadrature used to evaluate face integrals. dealii::hp::QCollection face_quadrature_collection; + /// 1D quadrature to generate Lagrange polynomials for the sake of flux interpolation. dealii::hp::QCollection<1> oned_quadrature_collection; protected: @@ -379,7 +416,7 @@ class DGBase const dealii::UpdateFlags face_update_flags = dealii::update_values | dealii::update_gradients | dealii::update_quadrature_points | dealii::update_JxW_values | dealii::update_normal_vectors; /// Update flags needed at neighbor' face points. /** NOTE: With hp-adaptation, might need to query neighbor's quadrature points depending on the order of the cells. */ - const dealii::UpdateFlags neighbor_face_update_flags = dealii::update_values | dealii::update_gradients; + const dealii::UpdateFlags neighbor_face_update_flags = dealii::update_values | dealii::update_gradients | dealii::update_quadrature_points | dealii::update_JxW_values; @@ -388,6 +425,28 @@ class DGBase dealii::ConditionalOStream pcout; ///< Parallel std::cout that only outputs on mpi_rank==0 private: + /// Evaluate the average penalty term at the face + /** For a cell with solution of degree p, and Hausdorff measure h, + * which represents the element dimension orthogonal to the face, + * the penalty term is given by p*(p+1)/h . + */ + template + real evaluate_penalty_scaling ( + const DoFCellAccessorType &cell, + const int iface, + const dealii::hp::FECollection fe_collection) const; + + /// In the case that two cells have the same coarseness, this function decides if the current cell should perform the work. + /** In the case the neighbor is a ghost cell, we let the processor with the lower rank do the work on that face. + * We cannot use the cell->index() because the index is relative to the distributed triangulation. + * Therefore, the cell index of a ghost cell might be different to the physical cell index even if they refer to the same cell. + * + * For a locally owned neighbor cell, cell with lower index does work or if both cells have same index, then cell at the lower level does the work + * See https://www.dealii.org/developer/doxygen/deal.II/classTriaAccessorBase.html#a695efcbe84fefef3e4c93ee7bdb446ad + */ + template + bool current_cell_should_do_the_work (const DoFCellAccessorType1 ¤t_cell, const DoFCellAccessorType2 &neighbor_cell) const; + /// Used in the delegated constructor /** The main reason we use this weird function is because all of the above objects * need to be looped with the various p-orders. This function allows us to do this in a @@ -404,22 +463,29 @@ template class DGWeak : public DGBase { #if PHILIP_DIM==1 // dealii::parallel::distributed::Triangulation does not work for 1D + /** Triangulation to store the grid. + * In 1D, dealii::Triangulation is used. + * In 2D, 3D, dealii::parallel::distributed::Triangulation is used. + */ using Triangulation = dealii::Triangulation; #else + /** Triangulation to store the grid. + * In 1D, dealii::Triangulation is used. + * In 2D, 3D, dealii::parallel::distributed::Triangulation is used. + */ using Triangulation = dealii::parallel::distributed::Triangulation; #endif public: - /// Constructor + /// Constructor. DGWeak( const Parameters::AllParameters *const parameters_input, const unsigned int degree, Triangulation *const triangulation_input); - ~DGWeak(); ///< Destructor + ~DGWeak(); ///< Destructor. private: - /// Contains the physics of the PDE std::shared_ptr < Physics::PhysicsBase > > pde_physics; /// Convective numerical flux @@ -497,8 +563,16 @@ template class DGStrong : public DGBase { #if PHILIP_DIM==1 // dealii::parallel::distributed::Triangulation does not work for 1D + /** Triangulation to store the grid. + * In 1D, dealii::Triangulation is used. + * In 2D, 3D, dealii::parallel::distributed::Triangulation is used. + */ using Triangulation = dealii::Triangulation; #else + /** Triangulation to store the grid. + * In 1D, dealii::Triangulation is used. + * In 2D, 3D, dealii::parallel::distributed::Triangulation is used. + */ using Triangulation = dealii::parallel::distributed::Triangulation; #endif public: @@ -591,8 +665,16 @@ template class DGFactory { #if PHILIP_DIM==1 // dealii::parallel::distributed::Triangulation does not work for 1D + /** Triangulation to store the grid. + * In 1D, dealii::Triangulation is used. + * In 2D, 3D, dealii::parallel::distributed::Triangulation is used. + */ using Triangulation = dealii::Triangulation; #else + /** Triangulation to store the grid. + * In 1D, dealii::Triangulation is used. + * In 2D, 3D, dealii::parallel::distributed::Triangulation is used. + */ using Triangulation = dealii::parallel::distributed::Triangulation; #endif public: diff --git a/src/dg/drdx.cpp b/src/dg/drdx.cpp new file mode 100644 index 000000000..2ba3922f1 --- /dev/null +++ b/src/dg/drdx.cpp @@ -0,0 +1,447 @@ +#include "dg.h" + +namespace PHiLiP { + +template +dealii::SparsityPattern DGBase::get_dRdX_sparsity_pattern () { + + const unsigned n_residuals = dof_handler.n_dofs(); + const unsigned n_nodes_coeff = high_order_grid.dof_handler_grid.n_dofs(); + const unsigned n_nodes_per_cell = high_order_grid.dof_handler_grid.get_fe_collection().max_dofs_per_cell(); + dealii::SparsityPattern sparsity_pattern(n_residuals, n_nodes_coeff, n_nodes_per_cell); + + std::vector resi_indices; + std::vector node_indices; + for (auto cell = dof_handler.begin_active(); cell != dof_handler.end(); ++cell) { + if (!cell->is_locally_owned()) continue; + + const unsigned int level = cell->level(); + const unsigned int index = cell->index(); + + const unsigned int n_resi_cell = fe_collection[cell->active_fe_index()].n_dofs_per_cell(); + resi_indices.resize(n_resi_cell); + cell->get_dof_indices (resi_indices); + + for (auto cell_grid = high_order_grid.dof_handler_grid.begin_active(); cell_grid != dof_handler.end(); ++cell_grid) { + // Brute force search for the same cell. + // There must be a better way + if (!cell_grid->is_locally_owned()) continue; + if (cell_grid->level() != level) continue; + if (cell_grid->index() != index) continue; + + const unsigned int n_node_cell = high_order_grid.fe_system.n_dofs_per_cell(); + node_indices.resize(n_node_cell); + cell_grid->get_dof_indices (node_indices); + + for (auto resi_row = resi_indices.begin(); resi_row!=resi_indices.end(); ++resi_row) { + sparsity_pattern.add_entries(*resi_row, node_indices.begin(), node_indices.end()); + } + + break; + } + + } // end of cell loop + sparsity_pattern.compress(); + + return sparsity_pattern; +} + +//template +//dealii::TrilinosWrappers::SparseMatrix DGBase::get_dRdX_finite_differences (dealii::SparsityPattern dRdX_sparsity_pattern) { +// +// const double pertubation = 1e-8; +// +// dealii::TrilinosWrappers::SparseMatrix dRdX; +// dRdX.reinit(locally_owned_dofs, dRdX_sparsity_pattern, mpi_communicator); +// +// // For now assume same polynomial degree across domain +// const unsigned int max_dofs_per_cell = dof_handler.get_fe_collection().max_dofs_per_cell(); +// std::vector current_dofs_indices(max_dofs_per_cell); +// std::vector neighbor_dofs_indices(max_dofs_per_cell); +// +// for (auto current_cell = dof_handler.begin_active(); current_cell != dof_handler.end(); ++current_cell) { +// if (!current_cell->is_locally_owned()) continue; +// +// +// // Current reference element related to this physical cell +// const unsigned int mapping_index = 0; +// const unsigned int fe_index_curr_cell = current_cell->active_fe_index(); +// const unsigned int quad_index = fe_index_curr_cell; +// const dealii::FESystem ¤t_fe_ref = fe_collection[fe_index_curr_cell]; +// const unsigned int curr_cell_degree = current_fe_ref.tensor_degree(); +// const unsigned int n_dofs_curr_cell = current_fe_ref.n_dofs_per_cell(); +// +// // Local vector contribution from each cell +// dealii::Vector current_cell_rhs (n_dofs_curr_cell); // Defaults to 0.0 initialization +// +// // Obtain the mapping from local dof indices to global dof indices +// current_dofs_indices.resize(n_dofs_curr_cell); +// current_cell->get_dof_indices (current_dofs_indices); +// +// +// // Cell diameter used to scale perturbation +// double cell_diameter = current_cell.diameter(); +// for (auto resi_row = current_dofs_indices.begin(); resi_row != current_dofs_indices.end(); ++resi_row) { +// const int irow_glob = *resi_row; +// const int n_cols = dRdX_sparsity_pattern.row_length(irow_glob); +// for (int icol=0; icol < n_cols; ++icol) { +// const int icol_glob = dRdX_sparsity_pattern.column_number(irow_glob, icol); +// +// double old_node = high_order_grid.nodes[icol_glob]; +// double dx = cell_diameter * perturbation; +// double high_order_grid.nodes[icol_glob] += dx; +// high_order_grid.nodes.update_ghost_values(); +// +// double high_order_grid.nodes[icol_glob] = old_node; +// high_order_grid.nodes.update_ghost_values(); +// } +// } +// +// //dealii::hp::MappingCollection mapping_collection(*(high_order_grid.mapping_fe_field)); +// //const dealii::MappingManifold mapping; +// //const dealii::MappingQ mapping(max_degree+1); +// const auto mapping = (*(high_order_grid.mapping_fe_field)); +// dealii::hp::MappingCollection mapping_collection(mapping); +// +// dealii::hp::FEValues fe_values_collection_volume (mapping_collection, fe_collection, volume_quadrature_collection, this->volume_update_flags); ///< FEValues of volume. +// dealii::hp::FEFaceValues fe_values_collection_face_int (mapping_collection, fe_collection, face_quadrature_collection, this->face_update_flags); ///< FEValues of interior face. +// dealii::hp::FEFaceValues fe_values_collection_face_ext (mapping_collection, fe_collection, face_quadrature_collection, this->neighbor_face_update_flags); ///< FEValues of exterior face. +// dealii::hp::FESubfaceValues fe_values_collection_subface (mapping_collection, fe_collection, face_quadrature_collection, this->face_update_flags); ///< FEValues of subface. +// +// dealii::hp::FEValues fe_values_collection_volume_lagrange (mapping_collection, fe_collection_lagrange, volume_quadrature_collection, this->volume_update_flags); +// +// // fe_values_collection.reinit(current_cell, quad_collection_index, mapping_collection_index, fe_collection_index) +// fe_values_collection_volume.reinit (current_cell, quad_index, mapping_index, fe_index_curr_cell); +// const dealii::FEValues &fe_values_volume = fe_values_collection_volume.get_present_fe_values(); +// +// +// dealii::TriaIterator> cell_iterator = static_cast> > (current_cell); +// //if (!(all_parameters->use_weak_form)) fe_values_collection_volume_lagrange.reinit (current_cell, quad_index, mapping_index, fe_index_curr_cell); +// fe_values_collection_volume_lagrange.reinit (cell_iterator, quad_index, mapping_index, fe_index_curr_cell); +// const dealii::FEValues &fe_values_lagrange = fe_values_collection_volume_lagrange.get_present_fe_values(); +// assemble_volume_terms_explicit (fe_values_volume, current_dofs_indices, current_cell_rhs, fe_values_lagrange); +// +// for (unsigned int iface=0; iface < dealii::GeometryInfo::faces_per_cell; ++iface) { +// +// auto current_face = current_cell->face(iface); +// auto neighbor_cell = current_cell->neighbor(iface); +// +// // See tutorial step-30 for breakdown of 4 face cases +// +// // Case 1: +// // Face at boundary +// if (current_face->at_boundary() && !current_cell->has_periodic_neighbor(iface) ) { +// +// fe_values_collection_face_int.reinit (current_cell, iface, quad_index, mapping_index, fe_index_curr_cell); +// +// if(current_face->at_boundary() && all_parameters->use_periodic_bc == true && dim == 1) //using periodic BCs (for 1d) +// { +// int cell_index = current_cell->index(); +// //int cell_index = current_cell->index(); +// if (cell_index == 0 && iface == 0) +// { +// fe_values_collection_face_int.reinit(current_cell, iface, quad_index, mapping_index, fe_index_curr_cell); +// neighbor_cell = dof_handler.begin_active(); +// for (unsigned int i = 0 ; i < triangulation->n_active_cells() - 1; ++i) +// { +// ++neighbor_cell; +// } +// neighbor_cell->get_dof_indices(neighbor_dofs_indices); +// const unsigned int fe_index_neigh_cell = neighbor_cell->active_fe_index(); +// const unsigned int quad_index_neigh_cell = fe_index_neigh_cell; +// const unsigned int mapping_index_neigh_cell = 0; +// +// fe_values_collection_face_ext.reinit(neighbor_cell,(iface == 1) ? 0 : 1,quad_index_neigh_cell,mapping_index_neigh_cell,fe_index_neigh_cell); +// +// } +// else if (cell_index == (int) triangulation->n_active_cells() - 1 && iface == 1) +// { +// fe_values_collection_face_int.reinit(current_cell, iface, quad_index, mapping_index, fe_index_curr_cell); +// neighbor_cell = dof_handler.begin_active(); +// neighbor_cell->get_dof_indices(neighbor_dofs_indices); +// const unsigned int fe_index_neigh_cell = neighbor_cell->active_fe_index(); +// const unsigned int quad_index_neigh_cell = fe_index_neigh_cell; +// const unsigned int mapping_index_neigh_cell = 0; +// fe_values_collection_face_ext.reinit(neighbor_cell,(iface == 1) ? 0 : 1, quad_index_neigh_cell, mapping_index_neigh_cell, fe_index_neigh_cell); //not sure how changing the face number would work in dim!=1-dimensions. +// } +// +// //std::cout << "cell " << current_cell->index() << "'s " << iface << "th face has neighbour: " << neighbor_cell->index() << std::endl; +// const int neighbor_face_no = (iface ==1) ? 0:1; +// const unsigned int fe_index_neigh_cell = neighbor_cell->active_fe_index(); +// +// const dealii::FEFaceValues &fe_values_face_int = fe_values_collection_face_int.get_present_fe_values(); +// const dealii::FEFaceValues &fe_values_face_ext = fe_values_collection_face_ext.get_present_fe_values(); +// +// const dealii::FESystem &neigh_fe_ref = fe_collection[fe_index_neigh_cell]; +// const unsigned int neigh_cell_degree = neigh_fe_ref.tensor_degree(); +// const unsigned int n_dofs_neigh_cell = neigh_fe_ref.n_dofs_per_cell(); +// +// dealii::Vector neighbor_cell_rhs (n_dofs_neigh_cell); // Defaults to 0.0 initialization +// +// +// const unsigned int normal_direction1 = dealii::GeometryInfo::unit_normal_direction[iface]; +// const unsigned int normal_direction2 = dealii::GeometryInfo::unit_normal_direction[neighbor_face_no]; +// const unsigned int deg1sq = (curr_cell_degree == 0) ? 1 : curr_cell_degree * (curr_cell_degree+1); +// const unsigned int deg2sq = (neigh_cell_degree == 0) ? 1 : neigh_cell_degree * (neigh_cell_degree+1); +// +// //const real vol_div_facearea1 = current_cell->extent_in_direction(normal_direction1) / current_face->number_of_children(); +// const real vol_div_facearea1 = current_cell->extent_in_direction(normal_direction1); +// const real vol_div_facearea2 = neighbor_cell->extent_in_direction(normal_direction2); +// +// const real penalty1 = deg1sq / vol_div_facearea1; +// const real penalty2 = deg2sq / vol_div_facearea2; +// +// real penalty = 0.5 * ( penalty1 + penalty2 ); +// +// assemble_face_term_explicit ( +// fe_values_face_int, fe_values_face_ext, +// penalty, +// current_dofs_indices, neighbor_dofs_indices, +// current_cell_rhs, neighbor_cell_rhs); +// +// } else { +// const dealii::FEFaceValues &fe_values_face_int = fe_values_collection_face_int.get_present_fe_values(); +// const unsigned int deg1sq = (curr_cell_degree == 0) ? 1 : curr_cell_degree * (curr_cell_degree+1); +// const unsigned int normal_direction = dealii::GeometryInfo::unit_normal_direction[iface]; +// const real vol_div_facearea1 = current_cell->extent_in_direction(normal_direction); +// +// real penalty = deg1sq / vol_div_facearea1; +// +// const unsigned int boundary_id = current_face->boundary_id(); +// // Need to somehow get boundary type from the mesh +// assemble_boundary_term_explicit (boundary_id, fe_values_face_int, penalty, current_dofs_indices, current_cell_rhs); +// } +// +// //CASE 1.5: periodic boundary conditions +// //note that periodicity is not adapted for hp adaptivity yet. this needs to be figured out in the future +// } else if (current_face->at_boundary() && current_cell->has_periodic_neighbor(iface)){ +// +// neighbor_cell = current_cell->periodic_neighbor(iface); +// //std::cout << "cell " << current_cell->index() << " at boundary" <index() << std::endl; +// +// +// if (!current_cell->periodic_neighbor_is_coarser(iface) && +// (neighbor_cell->index() > current_cell->index() || +// (neighbor_cell->index() == current_cell->index() && current_cell->level() < neighbor_cell->level()) +// ) +// ) +// { +// Assert (current_cell->periodic_neighbor(iface).state() == dealii::IteratorState::valid, dealii::ExcInternalError()); +// +// +// // Corresponding face of the neighbor. +// // e.g. The 4th face of the current cell might correspond to the 3rd face of the neighbor +// const unsigned int neighbor_face_no = current_cell->periodic_neighbor_of_periodic_neighbor(iface); +// +// // Get information about neighbor cell +// const unsigned int fe_index_neigh_cell = neighbor_cell->active_fe_index(); +// const unsigned int quad_index_neigh_cell = fe_index_neigh_cell; +// const unsigned int mapping_index_neigh_cell = 0; +// const dealii::FESystem &neigh_fe_ref = fe_collection[fe_index_neigh_cell]; +// const unsigned int neigh_cell_degree = neigh_fe_ref.tensor_degree(); +// const unsigned int n_dofs_neigh_cell = neigh_fe_ref.n_dofs_per_cell(); +// +// // Local rhs contribution from neighbor +// dealii::Vector neighbor_cell_rhs (n_dofs_neigh_cell); // Defaults to 0.0 initialization +// +// // Obtain the mapping from local dof indices to global dof indices for neighbor cell +// neighbor_dofs_indices.resize(n_dofs_neigh_cell); +// neighbor_cell->get_dof_indices (neighbor_dofs_indices); +// +// fe_values_collection_face_int.reinit (current_cell, iface, quad_index, mapping_index, fe_index_curr_cell); +// const dealii::FEFaceValues &fe_values_face_int = fe_values_collection_face_int.get_present_fe_values(); +// fe_values_collection_face_ext.reinit (neighbor_cell, neighbor_face_no, quad_index_neigh_cell, mapping_index_neigh_cell, fe_index_neigh_cell); +// const dealii::FEFaceValues &fe_values_face_ext = fe_values_collection_face_ext.get_present_fe_values(); +// +// const unsigned int normal_direction1 = dealii::GeometryInfo::unit_normal_direction[iface]; +// const unsigned int normal_direction2 = dealii::GeometryInfo::unit_normal_direction[neighbor_face_no]; +// const unsigned int deg1sq = (curr_cell_degree == 0) ? 1 : curr_cell_degree * (curr_cell_degree+1); +// const unsigned int deg2sq = (neigh_cell_degree == 0) ? 1 : neigh_cell_degree * (neigh_cell_degree+1); +// +// //const real vol_div_facearea1 = current_cell->extent_in_direction(normal_direction1) / current_face->number_of_children(); +// const real vol_div_facearea1 = current_cell->extent_in_direction(normal_direction1); +// const real vol_div_facearea2 = neighbor_cell->extent_in_direction(normal_direction2); +// +// const real penalty1 = deg1sq / vol_div_facearea1; +// const real penalty2 = deg2sq / vol_div_facearea2; +// +// real penalty = 0.5 * ( penalty1 + penalty2 ); +// //penalty = 1;//99; +// +// assemble_face_term_explicit ( +// fe_values_face_int, fe_values_face_ext, +// penalty, +// current_dofs_indices, neighbor_dofs_indices, +// current_cell_rhs, neighbor_cell_rhs); +// +// // Add local contribution from neighbor cell to global vector +// for (unsigned int i=0; ihas_children()) { +// +// Assert (current_cell->neighbor(iface).state() == dealii::IteratorState::valid, dealii::ExcInternalError()); +// +// // Obtain cell neighbour +// const unsigned int neighbor_face_no = current_cell->neighbor_face_no(iface); +// +// for (unsigned int subface_no=0; subface_no < current_face->number_of_children(); ++subface_no) { +// +// // Get neighbor on ith subface +// auto neighbor_cell = current_cell->neighbor_child_on_subface (iface, subface_no); +// // Since the neighbor cell is finer than the current cell, it should not have more children +// Assert (!neighbor_cell->has_children(), dealii::ExcInternalError()); +// +// // Get information about neighbor cell +// const unsigned int fe_index_neigh_cell = neighbor_cell->active_fe_index(); +// const unsigned int quad_index_neigh_cell = fe_index_neigh_cell; +// const unsigned int mapping_index_neigh_cell = 0; +// const dealii::FESystem &neigh_fe_ref = fe_collection[fe_index_neigh_cell]; +// const unsigned int neigh_cell_degree = neigh_fe_ref.tensor_degree(); +// const unsigned int n_dofs_neigh_cell = neigh_fe_ref.n_dofs_per_cell(); +// +// dealii::Vector neighbor_cell_rhs (n_dofs_neigh_cell); // Defaults to 0.0 initialization +// +// // Obtain the mapping from local dof indices to global dof indices for neighbor cell +// neighbor_dofs_indices.resize(n_dofs_neigh_cell); +// neighbor_cell->get_dof_indices (neighbor_dofs_indices); +// +// fe_values_collection_subface.reinit (current_cell, iface, subface_no, quad_index, mapping_index, fe_index_curr_cell); +// const dealii::FESubfaceValues &fe_values_face_int = fe_values_collection_subface.get_present_fe_values(); +// +// fe_values_collection_face_ext.reinit (neighbor_cell, neighbor_face_no, quad_index_neigh_cell, mapping_index_neigh_cell, fe_index_neigh_cell); +// const dealii::FEFaceValues &fe_values_face_ext = fe_values_collection_face_ext.get_present_fe_values(); +// +// const unsigned int normal_direction1 = dealii::GeometryInfo::unit_normal_direction[iface]; +// const unsigned int normal_direction2 = dealii::GeometryInfo::unit_normal_direction[neighbor_face_no]; +// const unsigned int deg1sq = (curr_cell_degree == 0) ? 1 : curr_cell_degree * (curr_cell_degree+1); +// const unsigned int deg2sq = (neigh_cell_degree == 0) ? 1 : neigh_cell_degree * (neigh_cell_degree+1); +// +// const real vol_div_facearea1 = current_cell->extent_in_direction(normal_direction1); +// const real vol_div_facearea2 = neighbor_cell->extent_in_direction(normal_direction2); +// +// const real penalty1 = deg1sq / vol_div_facearea1; +// const real penalty2 = deg2sq / vol_div_facearea2; +// +// real penalty = 0.5 * ( penalty1 + penalty2 ); +// +// assemble_face_term_explicit ( +// fe_values_face_int, fe_values_face_ext, +// penalty, +// current_dofs_indices, neighbor_dofs_indices, +// current_cell_rhs, neighbor_cell_rhs); +// // Add local contribution from neighbor cell to global vector +// for (unsigned int i=0; ineighbor_is_coarser(iface)) +// // In the case the neighbor is a ghost cell, we let the processor with the lower rank do the work on that face +// // We cannot use the cell->index() because the index is relative to the distributed triangulation +// // Therefore, the cell index of a ghost cell might be different to the physical cell index even if they refer to the same cell +// && neighbor_cell->is_ghost() +// && current_cell->subdomain_id() < neighbor_cell->subdomain_id() +// ) +// || +// ( !(current_cell->neighbor_is_coarser(iface)) +// // In the case the neighbor is a local cell, we let the cell with the lower index do the work on that face +// && neighbor_cell->is_locally_owned() +// && +// ( // Cell with lower index does work +// current_cell->index() < neighbor_cell->index() +// || +// // If both cells have same index +// // See https://www.dealii.org/developer/doxygen/deal.II/classTriaAccessorBase.html#a695efcbe84fefef3e4c93ee7bdb446ad +// // then cell at the lower level does the work +// (neighbor_cell->index() == current_cell->index() && current_cell->level() < neighbor_cell->level()) +// ) +// ) +// ) +// { +// Assert (current_cell->neighbor(iface).state() == dealii::IteratorState::valid, dealii::ExcInternalError()); +// +// auto neighbor_cell = current_cell->neighbor_or_periodic_neighbor(iface); +// // Corresponding face of the neighbor. +// // e.g. The 4th face of the current cell might correspond to the 3rd face of the neighbor +// const unsigned int neighbor_face_no = current_cell->neighbor_of_neighbor(iface); +// +// // Get information about neighbor cell +// const unsigned int fe_index_neigh_cell = neighbor_cell->active_fe_index(); +// const unsigned int quad_index_neigh_cell = fe_index_neigh_cell; +// const unsigned int mapping_index_neigh_cell = 0; +// const dealii::FESystem &neigh_fe_ref = fe_collection[fe_index_neigh_cell]; +// const unsigned int neigh_cell_degree = neigh_fe_ref.tensor_degree(); +// const unsigned int n_dofs_neigh_cell = neigh_fe_ref.n_dofs_per_cell(); +// +// // Local rhs contribution from neighbor +// dealii::Vector neighbor_cell_rhs (n_dofs_neigh_cell); // Defaults to 0.0 initialization +// +// // Obtain the mapping from local dof indices to global dof indices for neighbor cell +// neighbor_dofs_indices.resize(n_dofs_neigh_cell); +// neighbor_cell->get_dof_indices (neighbor_dofs_indices); +// +// fe_values_collection_face_int.reinit (current_cell, iface, quad_index, mapping_index, fe_index_curr_cell); +// const dealii::FEFaceValues &fe_values_face_int = fe_values_collection_face_int.get_present_fe_values(); +// fe_values_collection_face_ext.reinit (neighbor_cell, neighbor_face_no, quad_index_neigh_cell, mapping_index_neigh_cell, fe_index_neigh_cell); +// const dealii::FEFaceValues &fe_values_face_ext = fe_values_collection_face_ext.get_present_fe_values(); +// +// const unsigned int normal_direction1 = dealii::GeometryInfo::unit_normal_direction[iface]; +// const unsigned int normal_direction2 = dealii::GeometryInfo::unit_normal_direction[neighbor_face_no]; +// const unsigned int deg1sq = (curr_cell_degree == 0) ? 1 : curr_cell_degree * (curr_cell_degree+1); +// const unsigned int deg2sq = (neigh_cell_degree == 0) ? 1 : neigh_cell_degree * (neigh_cell_degree+1); +// +// //const real vol_div_facearea1 = current_cell->extent_in_direction(normal_direction1) / current_face->number_of_children(); +// const real vol_div_facearea1 = current_cell->extent_in_direction(normal_direction1); +// const real vol_div_facearea2 = neighbor_cell->extent_in_direction(normal_direction2); +// +// const real penalty1 = deg1sq / vol_div_facearea1; +// const real penalty2 = deg2sq / vol_div_facearea2; +// +// real penalty = 0.5 * ( penalty1 + penalty2 ); +// //penalty = 1;//99; +// +// assemble_face_term_explicit ( +// fe_values_face_int, fe_values_face_ext, +// penalty, +// current_dofs_indices, neighbor_dofs_indices, +// current_cell_rhs, neighbor_cell_rhs); +// +// // Add local contribution from neighbor cell to global vector +// for (unsigned int i=0; i + +#include + +// For metric Jacobian testing +#include +#include +#include +#include +// ***************** + #include +#include #include +#include +#include + +#include + #include +#include +#include + +#include + #include "high_order_grid.h" namespace PHiLiP { @@ -22,11 +44,38 @@ HighOrderGrid::HighOrderGrid( , mpi_communicator(MPI_COMM_WORLD) , pcout(std::cout, dealii::Utilities::MPI::this_mpi_process(mpi_communicator)==0) { + Assert(max_degree > 0, dealii::ExcMessage("Grid must be at least order 1.")); + nth_refinement = 0; allocate(); const dealii::ComponentMask mask(dim, true); - dealii::VectorTools::get_position_vector(dof_handler_grid, nodes, mask); + get_position_vector(dof_handler_grid, nodes, mask); nodes.update_ghost_values(); + update_surface_indices(); + update_surface_nodes(); mapping_fe_field = std::make_shared< dealii::MappingFEField > (dof_handler_grid,nodes,mask); + output_results_vtk(nth_refinement++); + + // Used to check Jacobian validity + const unsigned int exact_jacobian_order = (max_degree-1) * dim; + const unsigned int min_jacobian_order = 1; + const unsigned int used_jacobian_order = std::max(exact_jacobian_order, min_jacobian_order); + evaluate_lagrange_to_bernstein_operator(used_jacobian_order); + + auto cell = dof_handler_grid.begin_active(); + auto endcell = dof_handler_grid.end(); + std::cout << "Disabled check_valid_cells. Took too much time due to shape_grad()." << std::endl; + for (; cell!=endcell; ++cell) { + if (!cell->is_locally_owned()) continue; + const bool is_invalid_cell = check_valid_cell(cell); + + if ( !is_invalid_cell ) { + std::cout << " Poly: " << max_degree + << " Grid: " << nth_refinement + << " Cell: " << cell->active_cell_index() << " has an invalid Jacobian." << std::endl; + // bool fixed_invalid_cell = fix_invalid_cell(cell); + // if (fixed_invalid_cell) std::cout << "Fixed it." << std::endl; + } + } } template @@ -60,6 +109,625 @@ HighOrderGrid::allocate() // +template +void HighOrderGrid +::get_position_vector(const DoFHandlerType &dh, VectorType &vector, const dealii::ComponentMask &mask) +{ + AssertDimension(vector.size(), dh.n_dofs()); + const dealii::FESystem &fe = dh.get_fe(); + + // Construct default fe_mask; + const dealii::ComponentMask fe_mask(mask.size() ? mask : dealii::ComponentMask(fe.get_nonzero_components(0).size(), true)); + + AssertDimension(fe_mask.size(), fe.get_nonzero_components(0).size()); + + std::vector fe_to_real(fe_mask.size(), dealii::numbers::invalid_unsigned_int); + unsigned int size = 0; + for (unsigned int i = 0; i < fe_mask.size(); ++i) { + if (fe_mask[i]) fe_to_real[i] = size++; + } + Assert(size == dim, dealii::ExcMessage( + "The Component Mask you provided is invalid. It has to select exactly dim entries.")); + + const dealii::Quadrature quad(fe.get_unit_support_points()); + + dealii::MappingQ map_q(fe.degree); + dealii::FEValues fe_v(map_q, fe, quad, dealii::update_quadrature_points); + std::vector dofs(fe.dofs_per_cell); + + AssertDimension(fe.dofs_per_cell, fe.get_unit_support_points().size()); + Assert(fe.is_primitive(), dealii::ExcMessage("FE is not Primitive! This won't work.")); + + for (const auto &cell : dh.active_cell_iterators()) { + if (cell->is_locally_owned()) { + fe_v.reinit(cell); + cell->get_dof_indices(dofs); + const std::vector> &points = fe_v.get_quadrature_points(); + for (unsigned int q = 0; q < points.size(); ++q) { + const unsigned int comp = fe.system_to_component_index(q).first; + if (fe_mask[comp]) ::dealii::internal::ElementAccess::set(points[q][fe_to_real[comp]], dofs[q], vector); + } + } + } +} + + +template +template +inline real2 HighOrderGrid +::determinant(const std::array< dealii::Tensor<1,dim,real2>, dim > jacobian) const +{ + if constexpr(dim == 1) return jacobian[0][0]; + if constexpr(dim == 2) + return jacobian[0][0] * jacobian[1][1] - jacobian[1][0] * jacobian[0][1]; + if constexpr(dim == 3) + return (jacobian[0][0] * jacobian[1][1] * jacobian[2][2] - + jacobian[0][0] * jacobian[1][2] * jacobian[2][1] - + jacobian[1][0] * jacobian[0][1] * jacobian[2][2] + + jacobian[1][0] * jacobian[0][2] * jacobian[2][1] + + jacobian[2][0] * jacobian[0][1] * jacobian[1][2] - + jacobian[2][0] * jacobian[0][2] * jacobian[1][1]); +} + + +template +std::vector HighOrderGrid::evaluate_jacobian_at_points( + const VectorType &solution, + const typename DoFHandlerType::cell_iterator &cell, + const std::vector> &points) const +{ + const dealii::FESystem &fe = cell->get_fe(); + const unsigned int n_quad_pts = points.size(); + const unsigned int n_dofs_cell = fe.n_dofs_per_cell(); + const unsigned int n_axis = dim; + + std::vector dofs_indices(n_dofs_cell); + cell->get_dof_indices (dofs_indices); + + std::array coords; + std::array< dealii::Tensor<1,dim,real>, n_axis > coords_grad; + std::vector jac_det(n_quad_pts); + + for (unsigned int iquad=0; iquad phys_point; + for (int d=0;d +template +real2 HighOrderGrid::evaluate_jacobian_at_point( + const std::vector &dofs, + const dealii::FESystem &fe, + const dealii::Point &point) const +{ + AssertDimension(dim, fe.n_components()); + + const unsigned int n_dofs_coords = fe.n_dofs_per_cell(); + const unsigned int n_axis = dim; + + std::array< dealii::Tensor<1,dim,real2>, n_axis > coords_grad; // Tensor initialize with zeros + for (unsigned int idof=0; idof +template +void HighOrderGrid::evaluate_jacobian_at_points( + const std::vector &dofs, + const dealii::FESystem &fe, + const std::vector> &points, + std::vector &jacobian_determinants) const +{ + AssertDimension(dim, fe.n_components()); + AssertDimension(jacobian_determinants.size(), points.size()); + + const unsigned int n_points = points.size(); + + for (unsigned int i=0; i +std::vector matrix_stdvector_mult(const dealii::FullMatrix &matrix, const std::vector &vector_in) +{ + const unsigned int m = matrix.m(), n = matrix.n(); + AssertDimension(vector_in.size(),n); + std::vector vector_out(m,0.0); + for (unsigned int row=0; row +void HighOrderGrid::evaluate_lagrange_to_bernstein_operator(const unsigned int order) +{ + const dealii::FE_Q lagrange_basis(order); + const std::vector< dealii::Point > &lagrange_pts = lagrange_basis.get_unit_support_points(); + const unsigned int n_lagrange_pts = lagrange_pts.size(); + + const dealii::FE_Bernstein bernstein_basis(order); + const unsigned int n_bernstein = bernstein_basis.n_dofs_per_cell(); + AssertDimension(n_bernstein, n_lagrange_pts); + // Evaluate Vandermonde matrix such that V u_bernstein = u_lagrange + // where the matrix's rows correspond to the different Bernstein polynomials + // and the matrix's column correspond to the unit support points of the Lagrage polynomials + dealii::FullMatrix bernstein_to_lagrange(n_bernstein, n_lagrange_pts); + for (unsigned int ibernstein=0; ibernstein &point = lagrange_pts[ijacpts]; + bernstein_to_lagrange[ibernstein][ijacpts] = bernstein_basis.shape_value(ibernstein, point); + } + } + lagrange_to_bernstein_operator.reinit(n_lagrange_pts, n_bernstein); + std::cout << "Careful, about to invert a " << n_lagrange_pts << " x " << n_lagrange_pts << " dense matrix..." << std::endl; + lagrange_to_bernstein_operator.invert(bernstein_to_lagrange); + std::cout << "Done inverting a " << n_lagrange_pts << " x " << n_lagrange_pts << " dense matrix..." << std::endl; +} + + +template +bool HighOrderGrid::check_valid_cell(const typename DoFHandlerType::cell_iterator &cell) const +{ + return true; + const unsigned int exact_jacobian_order = (max_degree-1) * dim, min_jacobian_order = 1; + const unsigned int used_jacobian_order = std::max(exact_jacobian_order, min_jacobian_order); + + // Evaluate Jacobian at Lagrange interpolation points + const dealii::FESystem &fe_coords = cell->get_fe(); + const unsigned int n_dofs_coords = fe_coords.n_dofs_per_cell(); + std::vector dofs_indices(n_dofs_coords); + cell->get_dof_indices (dofs_indices); + + std::vector< real > cell_nodes(n_dofs_coords); + for (unsigned int idof = 0; idof < n_dofs_coords; ++idof) { + cell_nodes[idof] = nodes(dofs_indices[idof]); + } + + const dealii::FE_Q lagrange_basis(used_jacobian_order); + const std::vector< dealii::Point > &lagrange_pts = lagrange_basis.get_unit_support_points(); + const unsigned int n_lagrange_pts = lagrange_pts.size(); + const unsigned int n_bernstein = n_lagrange_pts; + std::vector lagrange_coeff(n_lagrange_pts); + evaluate_jacobian_at_points(cell_nodes, fe_coords, lagrange_pts, lagrange_coeff); + std::vector bernstein_coeff = matrix_stdvector_mult(lagrange_to_bernstein_operator, lagrange_coeff); + + const real tol = 1e-12; + for (unsigned int i=0; i lagrange_to_bernstein_operator( +// const dealii::FE_Q &lagrange_basis, +// const dealii::FE_Bernstein &bernstein_basis, +// const typename DoFHandlerType::cell_iterator &cell) +// { +// const dealii::FE_Bernstein bernstein_basis(jacobian_order); +// const unsigned int n_bernstein = bernstein_basis.n_dofs_per_cell(); +// const unsigned int n_lagrange_pts = lagrange_pts.size(); +// AssertDimension(n_bernstein, n_lagrange_pts); +// // Evaluate Vandermonde matrix such that V u_bernstein = u_lagrange +// // where the matrix's rows correspond to the different Bernstein polynomials +// // and the matrix's column correspond to the unit support points of the Lagrage polynomials +// dealii::FullMatrix bernstein_to_lagrange(n_bernstein, n_lagrange_pts); +// for (unsigned int ibernstein=0; ibernstein &point = lagrange_pts[ijacpts]; +// bernstein_to_lagrange[ibernstein][ijacpts] = bernstein_basis.shape_value(ibernstein, point); +// } +// } +// } + +template +bool HighOrderGrid::fix_invalid_cell(const typename DoFHandlerType::cell_iterator &cell) +{ + // This will be the target ratio between the current minimum (estimated) cell Jacobian + // and the value of the straight-sided Jacobian. This estimates depends on the Bernstein + // coefficients that serve as a convex hull + const double target_ratio = 0.1; + + // Maximum number of times we will move the barrier + const int max_barrier_iterations = 100; + + const unsigned int exact_jacobian_order = (max_degree-1) * dim, min_jacobian_order = 1; + const unsigned int used_jacobian_order = std::max(exact_jacobian_order, min_jacobian_order); + const dealii::FE_Q lagrange_basis(used_jacobian_order); + const std::vector< dealii::Point > &lagrange_pts = lagrange_basis.get_unit_support_points(); + const unsigned int n_lagrange_pts = lagrange_pts.size(), n_bernstein = n_lagrange_pts; + + const dealii::FESystem &fe_coords = cell->get_fe(); + const unsigned int n_dofs_coords = fe_coords.n_dofs_per_cell(); + // Evaluate Jacobian at Lagrange interpolation points + std::vector dofs_indices(n_dofs_coords); + cell->get_dof_indices (dofs_indices); + + // Use reverse mode for more efficiency + using ADtype = Sacado::Fad::DFad; + std::vector cell_nodes(n_dofs_coords); + std::vector lagrange_coeff(n_lagrange_pts); + std::vector bernstein_coeff(n_bernstein); + + // Count and tag movable nodes + std::vector< bool > movable(n_dofs_coords); + unsigned int n_movable_nodes = 0; + // for (unsigned int idof = 0; idof < n_dofs_coords; ++idof) { + // bool is_interior_node = true; + // for (unsigned int iface = 0; iface::faces_per_cell; ++iface) { + // is_interior_node = is_interior_node && !(fe_coords.has_support_on_face(idof, iface)); + // } + // movable[idof] = is_interior_node; + // if (is_interior_node) n_movable_nodes++; + // } + for (unsigned int idof = 0; idof < n_dofs_coords; ++idof) { + const bool is_movable = (idof/dim > 2*dim); + movable[idof] = is_movable; + if (is_movable) n_movable_nodes++; + } + + // Evaluate straight sided cell Jacobian. + // Note that we can't simply use the Triangulation's cell's vertices since + // MappingFEField does not preserves_vertex_locations() + const unsigned int n_vertices = dealii::GeometryInfo::vertices_per_cell; + const dealii::FE_Q straight_sided_elem_fe(1); + const dealii::FESystem straight_sided_elem_fesystem(straight_sided_elem_fe, dim); + std::vector straight_sided_dofs(straight_sided_elem_fesystem.n_dofs_per_cell()); + for (unsigned int ivertex = 0; ivertex < n_vertices; ++ivertex) { + const dealii::Point unit_vertex = dealii::GeometryInfo::unit_cell_vertex(ivertex); + const dealii::Point phys_vertex = mapping_fe_field->transform_unit_to_real_cell(cell, unit_vertex); + for (int d=0;d unit_cell_center; + unit_cell_center[0] = 0.5; + if constexpr (dim>=2) unit_cell_center[1] = 0.5; + if constexpr (dim>=3) unit_cell_center[2] = 0.5; + + //const double straight_sided_jacobian = cell->measure(); + const double straight_sided_jacobian = evaluate_jacobian_at_point(straight_sided_dofs, straight_sided_elem_fesystem, unit_cell_center); + + + // Initialize movable nodes, which are our design variables + dealii::Vector movable_nodes(n_movable_nodes); + unsigned int idesign = 0; + for (unsigned int idof = 0; idof < n_dofs_coords; ++idof) { + cell_nodes[idof] = nodes(dofs_indices[idof]); + if (movable[idof]) movable_nodes[idesign++] = cell_nodes[idof].val(); + } + evaluate_jacobian_at_points(cell_nodes, fe_coords, lagrange_pts, lagrange_coeff); + bernstein_coeff = matrix_stdvector_mult(lagrange_to_bernstein_operator, lagrange_coeff); + real min_ratio = -1.0; + for (int barrier_iterations = 0; barrier_iterations < max_barrier_iterations && min_ratio < target_ratio; ++barrier_iterations) { + min_ratio = bernstein_coeff[0].val(); + for (unsigned int i=0; i &, dealii::Vector &)> func = + [&](const dealii::Vector &movable_nodes, dealii::Vector &gradient) { + + unsigned int idesign = 0; + for (unsigned int idof = 0; idof < n_dofs_coords; ++idof) { + if (movable[idof]) { + cell_nodes[idof] = movable_nodes[idesign]; + cell_nodes[idof].diff(idesign, n_movable_nodes); + idesign++; + } + } + evaluate_jacobian_at_points(cell_nodes, fe_coords, lagrange_pts, lagrange_coeff); + bernstein_coeff = matrix_stdvector_mult(lagrange_to_bernstein_operator, lagrange_coeff); + + ADtype functional = 0.0; + min_ratio = bernstein_coeff[0].val(); + for (unsigned int i=0; i>::AdditionalData data(max_history, false); + // dealii::SolverBFGS> solver(solver_control, data); + // solver.connect_preconditioner_slot(preconditioner); + // solver.solve(func, movable_nodes); + const unsigned int max_opt_iterations = 1000; + const unsigned int n_line_search = 40; + const double initial_step_length = 1.0;//1e-3 * cell->minimum_vertex_distance(); + dealii::Vector old_movable_nodes(n_movable_nodes); + dealii::Vector search_direction(n_movable_nodes); + dealii::Vector grad(n_movable_nodes); + dealii::Vector old_grad(n_movable_nodes); + real functional = func(movable_nodes, grad); + const real initial_grad_norm = grad.l2_norm(); + double grad_norm = grad.l2_norm(); + dealii::FullMatrix hessian_inverse(n_movable_nodes); + dealii::FullMatrix outer_product_term(n_movable_nodes); + dealii::Vector dg(n_movable_nodes); + dealii::Vector dx(n_movable_nodes); + dealii::Vector B_dg(n_movable_nodes); + hessian_inverse = 0; + for (unsigned int inode=0; inode gradient_drop;++i) { + + old_movable_nodes = movable_nodes; + old_grad = grad; + + hessian_inverse.vmult(search_direction, grad); + search_direction *= -1.0; + double step_length = initial_step_length; + real old_functional = functional; + unsigned int i_line_search; + for (i_line_search = 0; i_line_search +// bool HighOrderGrid::make_valid_cell( +// const typename DoFHandlerType::cell_iterator &cell); +// { +// const dealii::FESystem &fe_coords = cell->get_fe(); +// const order = fe_coords.tensor_degree(); +// const jacobian_order = std::pow(fe_coords.tensor_degree()-1, dim); +// const dealii::FE_Q lagrange_basis(jacobian_order); +// const std::vector< dealii::Point > &jacobian_points = lagrange_basis.get_unit_support_points(); +// +// const dealii::FE_Bernstein bernstein_basis(jacobian_order); +// +// const unsigned int n_lagrange_pts = jacobian_points.size(); +// const unsigned int n_dofs_coords = fe_coords.n_dofs_per_cell(); +// const unsigned int n_axis = dim; +// +// // Evaluate Jacobian at Lagrange interpolation points +// std::vector dofs_indices(n_dofs_cell); +// cell->get_dof_indices (dofs_indices); +// +// dealii::Vector lagrange_coeff(n_lagrange_pts); +// +// for (unsigned int i_lagrange=0; i_lagrange &point = jacobian_points[i_lagrange]; +// +// std::array< dealii::Tensor<1,dim,real>, n_axis > > coords_grad; // Tensor initialize with zeros +// for (unsigned int idof=0; idof jacobian; +// for (unsigned int a=0;a lagrange_coeff(i_lagrange); +// } +// +// const unsigned int n_bernstein = bernstein_basis.n_dofs_per_cell(); +// AssertDimension(n_bernstein == n_lagrange_pts); +// // Evaluate Vandermonde matrix such that V u_bernstein = u_lagrange +// // where the matrix's rows correspond to the different Bernstein polynomials +// // and the matrix's column correspond to the unit support points of the Lagrage polynomials +// dealii::FullMatrix bernstein_to_lagrange(n_bernstein, n_lagrange_pts); +// for (unsigned int ibernstein=0; ibernstein &point = jacobian_points[ijacpts]; +// bernstein_to_lagrange[ibernstein][ijacpts] = bernstein_basis.shape_value(ibernstein, point); +// } +// } +// dealii::FullMatrix lagrange_to_bernstein; +// lagrange_to_bernstein.invert(bernstein_to_lagrange); +// +// dealii::Vector bernstein_coeff(n_bernstein); +// lagrange_to_bernstein.vmult(bernstein_coeff, lagrange_coeff); +// +// return false; +// } + +template +void HighOrderGrid::test_jacobian() +{ + // // Setup a dummy system + // const unsigned int solution_degree = max_degree-1; + // const unsigned int dummy_n_state = 5; + // const dealii::FE_DGQ fe_dgq(solution_degree); + // const dealii::FESystem fe_system_dgq(fe_dgq, dummy_n_state); + // const dealii::QGauss qgauss(solution_degree+1); + // const dealii::QGaussLobatto qgauss_lobatto(solution_degree+20); + + // dealii::DoFHandler dof_handler; + // dof_handler.initialize(*triangulation, fe_system_dgq); + + // const dealii::Mapping &mapping = (*(mapping_fe_field)); + // const dealii::FESystem &fe = fe_system_dgq; + // const dealii::Quadrature &quadrature = qgauss; + // const dealii::UpdateFlags volume_update_flags = + // dealii::update_values + // | dealii::update_gradients + // | dealii::update_quadrature_points + // | dealii::update_JxW_values + // | dealii::update_jacobians + // ; + + // { + // const dealii::Quadrature &overquadrature = qgauss_lobatto; + // const std::vector< dealii::Point< dim > > points = overquadrature.get_points(); + // std::vector jac_det; + // for (auto cell = dof_handler_grid.begin_active(); cell!=dof_handler_grid.end(); ++cell) { + // if (!cell->is_locally_owned()) continue; + // jac_det = evaluate_jacobian_at_points(nodes, cell, points); + // } + // } + + // dealii::FEValues fe_values(mapping, fe, quadrature, volume_update_flags); + + // for (auto cell = dof_handler.begin_active(); cell!=dof_handler.end(); ++cell) { + + // if (!cell->is_locally_owned()) continue; + + // // const unsigned int n_quad_pts = quadrature.size(); + + // // std::cout << " Cell " << cell->active_cell_index(); + // // fe_values.reinit(cell); + // // for (unsigned int iquad=0; iquad> &points = fe_values.get_quadrature_points(); + // // std::cout << " Point: " << points[iquad] + // // << " JxW: " << fe_values.JxW(iquad) + // // << " J: " << (fe_values.jacobian(iquad)).determinant() + // // << std::endl; + // // } + + // // for (unsigned int itrial=0; itrialget_dof_indices (dofs_indices); + // } +} + + template void HighOrderGrid::prepare_for_coarsening_and_refinement() { @@ -73,17 +741,488 @@ void HighOrderGrid::prepare_for_coarsening_a } template -void HighOrderGrid::execute_coarsening_and_refinement() { +void HighOrderGrid::execute_coarsening_and_refinement(const bool output_mesh) { allocate(); -#if PHILIP_DIM==1 // dealii::parallel::distributed::Triangulation does not work for 1D +#if PHILIP_DIM==1 solution_transfer.interpolate(old_nodes, nodes); #else solution_transfer.interpolate(nodes); #endif nodes.update_ghost_values(); + + dealii::AffineConstraints hanging_node_constraints; + hanging_node_constraints.clear(); + dealii::DoFTools::make_hanging_node_constraints(dof_handler_grid, hanging_node_constraints); + hanging_node_constraints.close(); + hanging_node_constraints.distribute(nodes); + + nodes.update_ghost_values(); + + update_surface_indices(); + update_surface_nodes(); mapping_fe_field = std::make_shared< dealii::MappingFEField > (dof_handler_grid,nodes); + if (output_mesh) output_results_vtk(nth_refinement++); + + auto cell = dof_handler_grid.begin_active(); + auto endcell = dof_handler_grid.end(); + std::cout << "Disabled check_valid_cells. Took too much time due to shape_grad()." << std::endl; + for (; cell!=endcell; ++cell) { + if (!cell->is_locally_owned()) continue; + const bool is_invalid_cell = check_valid_cell(cell); + + if ( !is_invalid_cell ) { + std::cout << " Poly: " << max_degree + << " Grid: " << nth_refinement + << " Cell: " << cell->active_cell_index() << " has an invalid Jacobian." << std::endl; + //bool fixed_invalid_cell = fix_invalid_cell(cell); + //if (fixed_invalid_cell) std::cout << "Fixed it." << std::endl; + } + } +} + + +template +void HighOrderGrid::deform_mesh(std::vector local_surface_displacements) { + + int n_mpi; + int rank; + MPI_Comm_rank(MPI_COMM_WORLD, &rank); + MPI_Comm_size(MPI_COMM_WORLD, &n_mpi); + + (void) local_surface_displacements; + // const unsigned int n_local_surface_disp = local_surface_displacements.size(); + // Assert(n_local_surface_disp==local_surface_nodes.size(), dealii::ExcDimensionMismatch(n_local_surface_disp,local_surface_nodes.size())); + + const unsigned int n_surface_nodes = all_surface_nodes.size(); + dealii::Vector surface_displacements(n_surface_nodes); + Assert(surface_displacements.size() == all_surface_nodes.size(), dealii::ExcDimensionMismatch(surface_displacements.size(),all_surface_nodes.size())); + + const unsigned int n_surface_points = n_surface_nodes / dim; + Assert( n_surface_nodes % dim == 0, dealii::ExcMessage("Surface nodes has incorrect size.")); + + (void) surface_displacements; + const double support_radius = 1.0; + const double support_radius2 = support_radius*support_radius; + + dealii::SparsityPattern sparsity_pattern_M(n_surface_points, n_surface_points, n_surface_points); + int row=0, col=0; + for (auto node1 = all_surface_nodes.begin(); node1 != all_surface_nodes.end(); node1+=dim) { + for (auto node2 = node1; node2 != all_surface_nodes.end(); node2+=dim) { + double distance2 = 0; + // Evaluate the squared distance + for (int d=0;d +std::vector flatten(const std::vector>& v) { + std::size_t total_size = 0; + for (const auto& sub : v) { + total_size += sub.size(); + } + std::vector result; + result.reserve(total_size); + for (const auto& sub : v) { + result.insert(result.end(), sub.begin(), sub.end()); + } + return result; +} + +//std::vector concatenate_vectors_across_mpi(std::vector vector) { +// +// const unsigned int n_local_entries = locally_owned_entries_indices.size(); +// std::vector n_local_entries_per_mpi(n_mpi); +// MPI_Allgather(&n_local_entries, 1, MPI::UNSIGNED, &(n_local_entries_per_mpi[0]), 1, MPI::UNSIGNED, MPI_COMM_WORLD); +// +// std::vector> all_local_entries(n_mpi); +// std::vector request(n_mpi); +// for (int i_mpi=0; i_mpi +void HighOrderGrid::update_surface_nodes() { + int n_mpi; + int rank; + MPI_Comm_rank(MPI_COMM_WORLD, &rank); + MPI_Comm_size(MPI_COMM_WORLD, &n_mpi); + + const unsigned int n_local_surface_nodes = locally_owned_surface_nodes_indices.size(); + // Copy surface node locations + local_surface_nodes.clear(); + local_surface_nodes.resize(n_local_surface_nodes); + unsigned int i = 0; + for (auto index = locally_owned_surface_nodes_indices.begin(); index != locally_owned_surface_nodes_indices.end(); index++) { + local_surface_nodes[i++] = nodes[*index]; + } + std::vector n_local_surface_nodes_per_mpi(n_mpi); + MPI_Allgather(&n_local_surface_nodes, 1, MPI::UNSIGNED, &(n_local_surface_nodes_per_mpi[0]), 1, MPI::UNSIGNED, MPI_COMM_WORLD); + n_surface_nodes = 0; + for (int i_mpi=0; i_mpi> all_local_surface_nodes(n_mpi); + std::vector request(n_mpi); + for (int i_mpi=0; i_mpi +void HighOrderGrid::update_surface_indices() { + locally_owned_surface_nodes_indices.clear(); + + const unsigned int dofs_per_cell = fe_system.n_dofs_per_cell(); + const unsigned int dofs_per_face = fe_system.n_dofs_per_face(); + std::vector< dealii::types::global_dof_index > dof_indices(dofs_per_cell); + + // Use unordered sets which uses hashmap. + // Has an average complexity of O(1) (worst case O(n)) for finding and inserting + // Overall algorithm average will be O(n), worst case O(n^2) + std::unordered_set locally_owned_dofs(locally_owned_dofs_grid.begin(), locally_owned_dofs_grid.end()); + std::unordered_set surface_dof_indices_temp; + for (auto cell = dof_handler_grid.begin_active(); cell!=dof_handler_grid.end(); ++cell) { + + if (!cell->is_locally_owned()) continue; + if (!cell->at_boundary()) continue; + + cell->get_dof_indices(dof_indices); + + for (unsigned int iface=0; iface < dealii::GeometryInfo::faces_per_cell; ++iface) { + const auto face = cell->face(iface); + if (!face->at_boundary()) continue; + + const unsigned int boundary_id = face->boundary_id(); + + for (unsigned int idof_face=0; idof_face +void HighOrderGrid::output_results_vtk (const unsigned int cycle) const +{ + +// { +// dealii::DoFHandler dof_handler_fine(*triangulation); +// const unsigned int jacobian_degree = std::pow((max_degree-1),dim); +// +// unsigned int output_degree = std::max(max_degree, jacobian_degree); +// output_degree = max_degree; +// +// const dealii::FE_Q fe_q_fine(output_degree); +// const dealii::FESystem fe_system_fine(fe_q_fine, dim); +// dof_handler_fine.initialize(*triangulation, fe_system_fine); +// dof_handler_fine.distribute_dofs(fe_system_fine); +// +// VectorType nodes_fine; +// #if PHILIP_DIM==1 // dealii::parallel::distributed::Triangulation does not work for 1D +// nodes_fine.reinit(dof_handler_fine.n_dofs()); +// #else +// dealii::IndexSet locally_owned_dofs_fine = dof_handler_fine.locally_owned_dofs(); +// dealii::IndexSet ghost_dofs_fine; +// dealii::DoFTools::extract_locally_relevant_dofs(dof_handler_fine, ghost_dofs_fine); +// dealii::IndexSet locally_relevant_dofs_fine = ghost_dofs_fine; +// ghost_dofs_fine.subtract_set(locally_owned_dofs_fine); +// nodes_fine.reinit(locally_owned_dofs_fine, ghost_dofs_fine, mpi_communicator); +// #endif +// dealii::FullMatrix interpolation_matrix(fe_system_fine.dofs_per_cell, fe_system.dofs_per_cell); +// fe_system_fine.get_interpolation_matrix(fe_system, interpolation_matrix); +// //dealii::FullMatrix interpolation_matrix(fe_system.dofs_per_cell, fe_system_fine.dofs_per_cell); +// //fe_system.get_interpolation_matrix(fe_system_fine, interpolation_matrix); +// dealii::VectorTools::interpolate(dof_handler_grid, dof_handler_fine, interpolation_matrix, nodes, nodes_fine); +// +// dealii::DataOut> data_out; +// data_out.attach_dof_handler (dof_handler_fine); +// +// std::vector solution_names; +// for(int d=0;d data_component_interpretation(dim, dealii::DataComponentInterpretation::component_is_scalar); +// data_out.add_data_vector (nodes_fine, solution_names, dealii::DataOut::type_dof_data, data_component_interpretation); +// +// dealii::Vector subdomain(triangulation->n_active_cells()); +// for (unsigned int i = 0; i < subdomain.size(); ++i) { +// subdomain[i] = triangulation->locally_owned_subdomain(); +// } +// data_out.add_data_vector(subdomain, "subdomain", dealii::DataOut_DoFData,dim>::DataVectorType::type_cell_data); +// +// VectorType jacobian_determinant; +// #if PHILIP_DIM==1 // dealii::parallel::distributed::Triangulation does not work for 1D +// jacobian_determinant.reinit(dof_handler_fine.n_dofs()); +// #else +// jacobian_determinant.reinit(locally_owned_dofs_fine, ghost_dofs_fine, mpi_communicator); +// #endif +// const unsigned int n_dofs_per_cell = fe_system_fine.n_dofs_per_cell(); +// std::vector dofs_indices(n_dofs_per_cell); +// +// const std::vector< dealii::Point > &points = fe_system_fine.get_unit_support_points(); +// std::vector jac_det; +// for (auto cell = dof_handler_fine.begin_active(); cell!=dof_handler_fine.end(); ++cell) { +// if (!cell->is_locally_owned()) continue; +// jac_det = evaluate_jacobian_at_points(nodes_fine, cell, points); +// cell->get_dof_indices (dofs_indices); +// for (unsigned int i=0; i> data_out; + data_out.attach_dof_handler (dof_handler_grid); + + std::vector solution_names; + for(int d=0;d data_component_interpretation(dim, dealii::DataComponentInterpretation::component_is_scalar); + data_out.add_data_vector (nodes, solution_names, dealii::DataOut::type_dof_data, data_component_interpretation); + + dealii::Vector subdomain(triangulation->n_active_cells()); + for (unsigned int i = 0; i < subdomain.size(); ++i) { + subdomain[i] = triangulation->locally_owned_subdomain(); + } + data_out.add_data_vector(subdomain, "subdomain", dealii::DataOut_DoFData,dim>::DataVectorType::type_cell_data); + + // const GridPostprocessor grid_post_processor; + // data_out.add_data_vector (nodes, grid_post_processor); + + VectorType jacobian_determinant; +#if PHILIP_DIM==1 // dealii::parallel::distributed::Triangulation does not work for 1D + jacobian_determinant.reinit(dof_handler_grid.n_dofs()); +#else + jacobian_determinant.reinit(locally_owned_dofs_grid, ghost_dofs_grid, mpi_communicator); +#endif + // const unsigned int n_dofs_per_cell = fe_system.n_dofs_per_cell(); + // std::vector dofs_indices(n_dofs_per_cell); + // const std::vector< dealii::Point > &points = fe_system.get_unit_support_points(); + // std::vector jac_det; + // for (auto cell = dof_handler_grid.begin_active(); cell!=dof_handler_grid.end(); ++cell) { + // if (!cell->is_locally_owned()) continue; + // jac_det = evaluate_jacobian_at_points(nodes, cell, points); + // cell->get_dof_indices (dofs_indices); + // for (unsigned int i=0; i jacobian_name; + for (int d=0;d jacobian_component_interpretation(dim, dealii::DataComponentInterpretation::component_is_scalar); + data_out.add_data_vector (jacobian_determinant, jacobian_name, dealii::DataOut::type_dof_data, jacobian_component_interpretation); + + typename dealii::DataOut::CurvedCellRegion curved = dealii::DataOut::CurvedCellRegion::curved_inner_cells; + //typename dealii::DataOut::CurvedCellRegion curved = dealii::DataOut::CurvedCellRegion::curved_boundary; + //typename dealii::DataOut::CurvedCellRegion curved = dealii::DataOut::CurvedCellRegion::no_curved_cells; + + const dealii::Mapping &mapping = (*(mapping_fe_field)); + const int n_subdivisions = max_degree;;//+30; // if write_higher_order_cells, n_subdivisions represents the order of the cell + data_out.build_patches(mapping, n_subdivisions, curved); + const bool write_higher_order_cells = (dim>1) ? true : false; + dealii::DataOutBase::VtkFlags vtkflags(0.0,cycle,true,dealii::DataOutBase::VtkFlags::ZlibCompressionLevel::best_compression,write_higher_order_cells); + data_out.set_flags(vtkflags); + + const int iproc = dealii::Utilities::MPI::this_mpi_process(mpi_communicator); + std::string filename = "Mesh-" + dealii::Utilities::int_to_string(dim, 1) +"D_GridP"+dealii::Utilities::int_to_string(max_degree, 2)+"-"; + filename += dealii::Utilities::int_to_string(cycle, 4) + "." + dealii::Utilities::int_to_string(iproc, 4); + filename += ".vtu"; + std::ofstream output(filename); + data_out.write_vtu(output); + if (iproc == 0) { + std::vector filenames; + for (unsigned int iproc = 0; iproc < dealii::Utilities::MPI::n_mpi_processes(mpi_communicator); ++iproc) { + std::string fn = "Mesh-" + dealii::Utilities::int_to_string(dim, 1) +"D_GridP"+dealii::Utilities::int_to_string(max_degree, 2)+"-"; + fn += dealii::Utilities::int_to_string(cycle, 4) + "." + dealii::Utilities::int_to_string(iproc, 4); + fn += ".vtu"; + filenames.push_back(fn); + } + std::ofstream master_output(master_fn); + data_out.write_pvtu_record(master_output, filenames); + } + +} + +template +void GridPostprocessor::evaluate_vector_field ( + const dealii::DataPostprocessorInputs::Vector &inputs, + std::vector> &computed_quantities) const +{ + const unsigned int n_solution_points = inputs.solution_values.size(); + Assert (computed_quantities.size() == n_solution_points, dealii::ExcInternalError()); + Assert (inputs.solution_values[0].size() == dim, dealii::ExcInternalError()); + std::vector names = get_names (); + for (unsigned int q=0; q &uh = inputs.solution_values[q]; + const std::vector > &duh = inputs.solution_gradients[q]; + // const std::vector > &dduh = inputs.solution_hessians[q]; + // const dealii::Tensor<1,dim> &normals = inputs.normals[q]; + // const dealii::Point &evaluation_points = inputs.evaluation_points[q]; + + unsigned int current_data_index = 0; + + dealii::Tensor<2,dim,double> jacobian; + for (unsigned int a=0;a +std::vector GridPostprocessor +::get_data_component_interpretation () const +{ + namespace DCI = dealii::DataComponentInterpretation; + std::vector interpretation; + interpretation.push_back (DCI::component_is_scalar); // Jacobian + + std::vector names = get_names(); + if (names.size() != interpretation.size()) { + std::cout << "Number of DataComponentInterpretation is not the same as number of names for output file" << std::endl; + } + return interpretation; +} + + +template +std::vector GridPostprocessor::get_names () const +{ + std::vector names; + names.push_back ("JacobianDeterminant"); + return names; +} + +template +dealii::UpdateFlags GridPostprocessor::get_needed_update_flags () const +{ + //return update_values | update_gradients; + return dealii::update_values | dealii::update_gradients; +} //template class HighOrderGrid; #if PHILIP_DIM==1 // dealii::parallel::distributed::Triangulation does not work for 1D diff --git a/src/dg/high_order_grid.h b/src/dg/high_order_grid.h index 30ea31638..28e4058ac 100644 --- a/src/dg/high_order_grid.h +++ b/src/dg/high_order_grid.h @@ -15,6 +15,10 @@ #include "parameters/all_parameters.h" +#include +#include +#include + namespace PHiLiP { /** This HighOrderGrid class basically contains all the different part necessary to generate @@ -35,10 +39,25 @@ template does not work for 1D + /** Vector class used to store the solution. + * In 1D dealii::Vector is used. + * In 2D, 3D, dealii::LinearAlgebra::distributed::Vector is used. + */ using Vector = dealii::Vector; + /** In 1D SolutionTransfer = dealii::SolutionTransfer, dealii::DoFHandler> is used. + * In 2D, 3D, dealii::parallel::distributed::SolutionTransfer, dealii::DoFHandler> is used. + */ using SolutionTransfer = dealii::SolutionTransfer, dealii::DoFHandler>; #else + /** Vector class used to store the solution. + * In 1D dealii::Vector is used. + * In 2D, 3D, dealii::LinearAlgebra::distributed::Vector is used. + */ using Vector = dealii::LinearAlgebra::distributed::Vector; + /// SolutionTransfer class used during refinement. + /** In 1D SolutionTransfer = dealii::SolutionTransfer, dealii::DoFHandler> is used. + * In 2D, 3D, dealii::parallel::distributed::SolutionTransfer, dealii::DoFHandler> is used. + */ using SolutionTransfer = dealii::parallel::distributed::SolutionTransfer, dealii::DoFHandler>; #endif public: @@ -65,37 +84,102 @@ class HighOrderGrid dealii::DoFHandler dof_handler_grid; /// Current nodal coefficients of the high-order grid. + /** Note that this contains all \ directions. + * By convention, the DoF representing the z-direction follows the DoF representing + * the y-direction, which follows the one representing the x-direction such that + * the integer division "idof_index / dim" gives the coordinates related to the same + * point. + */ Vector nodes; - /// List of surface points - Vector surface_nodes; + /// List of surface nodes. + /** Note that this contains all \ directions. + * By convention, the DoF representing the z-direction follows the DoF representing + * the y-direction, which follows the one representing the x-direction such that + * the integer division "idof_index / dim" gives the coordinates related to the same + * point. + */ + std::vector local_surface_nodes; + /// List of surface node indices + std::vector locally_owned_surface_nodes_indices; + /// List of surface node boundary IDs, corresponding to locally_owned_surface_nodes_indices + std::vector locally_owned_surface_nodes_boundary_id; + + void test_jacobian(); ///< Test metric Jacobian + + /// Evaluates Jacobian of a cell given some points using a global solution vector + std::vector evaluate_jacobian_at_points( + const VectorType &solution, + const typename DoFHandlerType::cell_iterator &cell, + const std::vector> &points) const; + + /// Evaluates Jacobian given some DoF, associated FE, and some points. + /** Calls evaluate_jacobian_at_point() on the vector of points. + * Can be used in conjunction with AD since the type is templated. + */ + template + void evaluate_jacobian_at_points( + const std::vector &dofs, + const dealii::FESystem &fe, + const std::vector> &points, + std::vector &jacobian_determinants) const; + /// Evaluates Jacobian given some DoF, associated FE, and some point. + /** Can be used in conjunction with AD since the type is templated. + */ + template + real2 evaluate_jacobian_at_point( + const std::vector &dofs, + const dealii::FESystem &fe, + const dealii::Point &point) const; + + /// Evaluate exact Jacobian determinant polynomial and uses Bernstein polynomials to determine positivity + bool check_valid_cell(const typename DoFHandlerType::cell_iterator &cell) const; + + /// Evaluate exact Jacobian determinant polynomial and uses Bernstein polynomials to determine positivity + bool fix_invalid_cell(const typename DoFHandlerType::cell_iterator &cell); + + /// Used to transform coefficients from a Lagrange basis to a Bernstein basis + dealii::FullMatrix lagrange_to_bernstein_operator; + /// Evaluates the operator to obtain Bernstein coefficients from a set of Lagrange coefficients + /** This is used in the evaluation of the Jacobian positivity by checking the convex hull of the + * resulting Bezier curve. + */ + void evaluate_lagrange_to_bernstein_operator(const unsigned int order); + + /// Update surface indices + void update_surface_indices(); + /// Update surface indices + void update_surface_nodes(); + + void output_results_vtk (const unsigned int cycle) const; ///< Output mesh with metric informations + + /// List of all surface nodes + std::vector all_surface_nodes; + unsigned int n_surface_nodes; ///< Total number of surface nodes /// RBF mesh deformation - To be done - void deform_mesh(Vector surface_displacements); - - /// Evaluate cell metric Jacobian - /** The metric Jacobian is given by the gradient of the physical location - * with respect to the reference locations - * \f[ J_{ij} = \frac{ - * ( \mathbf{F}_{conv}( u ) - * + \mathbf{F}_{diss}( u, \boldsymbol{\nabla}(u) ) - * = s(\mathbf{x}) - * \f] - */ - //dealii::Tensor<2,dim,real> cell_jacobian (const typename dealii::Triangulation::cell_iterator &cell, const dealii::Point &point) const override - //{ - //} - // - - - /** Prepares the solution transfer such that the curved refined grid is on top of the curved coarse grid. - * This function needs to be called before Triangulation::execute_coarsening_and_refinement() or Triangulation::refine_global() + //void deform_mesh(Vector surface_displacements); + void deform_mesh(std::vector local_surface_displacements); + + // /// Evaluate cell metric Jacobian + // /** The metric Jacobian is given by the gradient of the physical location + // * with respect to the reference locations + // */ + // dealii::Tensor<2,dim,real> cell_jacobian (const typename dealii::Triangulation::cell_iterator &cell, const dealii::Point &point) const override + // { + // } + + + /// Prepares the solution transfer such that the curved refined grid is on top of the curved coarse grid. + /** This function needs to be called before dealii::Triangulation::execute_coarsening_and_refinement() or dealii::Triangulation::refine_global() + * and this->execute_coarsening_and_refinement(). */ void prepare_for_coarsening_and_refinement(); - /** Transfers the coarse curved curve onto the fine curved grid. - * This function needs to be called after Triangulation::execute_coarsening_and_refinement() or Triangulation::refine_global() + /// Executes the solution transfer such that the curved refined grid is on top of the curved coarse grid. + /** This function needs to be after this->prepare_for_coarsening_and_refinement() and + * dealii::Triangulation::execute_coarsening_and_refinement() or dealii::Triangulation::refine_global(). */ - void execute_coarsening_and_refinement(); + void execute_coarsening_and_refinement(const bool output_mesh = false); /// Use Lagrange polynomial to represent the spatial location. const dealii::FE_Q fe_q; @@ -103,10 +187,10 @@ class HighOrderGrid const dealii::FESystem fe_system; - /** MappingFEField that will provide the polynomial-based grid. - * It is a shared smart pointer because the constructor requires the dof_handler_grid to be properly initialized. - * See discussion in the following thread: - * https://stackoverflow.com/questions/7557153/defining-an-object-without-calling-its-constructor-in-c + /// MappingFEField that will provide the polynomial-based grid. + /** It is a shared smart pointer because the constructor requires the dof_handler_grid to be properly initialized. + * See discussion in the following + * thread. */ std::shared_ptr> mapping_fe_field; @@ -114,6 +198,7 @@ class HighOrderGrid dealii::IndexSet ghost_dofs_grid; ///< Locally relevant ghost degrees of freedom for the grid dealii::IndexSet locally_relevant_dofs_grid; ///< Union of locally owned degrees of freedom and relevant ghost degrees of freedom for the grid protected: + int nth_refinement; ///< Used to name the various files outputted. /// Used for the SolutionTransfer when performing grid adaptation. Vector old_nodes; @@ -126,6 +211,34 @@ class HighOrderGrid MPI_Comm mpi_communicator; ///< MPI communicator dealii::ConditionalOStream pcout; ///< Parallel std::cout that only outputs on mpi_rank==0 + /// Evaluate the determinant of a matrix given in the format of a std::array,dim>. + /** The indices of the array represent the matrix rows, and the indices of the Tensor represents its columns. + */ + template + real2 determinant(const std::array< dealii::Tensor<1,dim,real2>, dim > jacobian) const; + + /// A stripped down copy of dealii::VectorTools::get_position_vector() + void get_position_vector(const DoFHandlerType &dh, VectorType &vector, const dealii::ComponentMask &mask); + + +}; + +/// Postprocessor used to output the grid. +template +class GridPostprocessor : public dealii::DataPostprocessor +{ +public: + // /// Constructor + // GridPostprocessor(); + + /// Evaluates the values of interest to output. + virtual void evaluate_vector_field (const dealii::DataPostprocessorInputs::Vector &inputs, std::vector> &computed_quantities) const override; + /// Returns the names associated with the output data. + virtual std::vector get_names () const override; + /// Returns the DCI associated with the output data. + virtual std::vector get_data_component_interpretation () const override; + /// Returns the update flags required to evaluate the output data. + virtual dealii::UpdateFlags get_needed_update_flags () const override; }; } // namespace PHiLiP diff --git a/src/dg/strong_dg.cpp b/src/dg/strong_dg.cpp index 9bfb24c87..d7313f07e 100644 --- a/src/dg/strong_dg.cpp +++ b/src/dg/strong_dg.cpp @@ -49,6 +49,8 @@ DGStrong::~DGStrong () pcout << "Destructing DGStrong..." << std::endl; delete conv_num_flux; delete diss_num_flux; + delete conv_num_flux_double; + delete diss_num_flux_double; } template diff --git a/src/dg/weak_dg.cpp b/src/dg/weak_dg.cpp index bda005340..2503d3ade 100644 --- a/src/dg/weak_dg.cpp +++ b/src/dg/weak_dg.cpp @@ -46,6 +46,9 @@ DGWeak::~DGWeak () pcout << "Destructing DGWeak..." << std::endl; delete conv_num_flux; delete diss_num_flux; + + delete conv_num_flux_double; + delete diss_num_flux_double; } template @@ -293,8 +296,10 @@ void DGWeak::assemble_face_term_implicit( AssertDimension (n_dofs_ext, dof_indices_ext.size()); // Jacobian and normal should always be consistent between two elements - // even for non-conforming meshes? - const std::vector &JxW_int = fe_values_int.get_JxW_values (); + // In the case of the non-conforming mesh, we should be using the Jacobian + // of the smaller face since it would be "half" of the larger one. + // However, their curvature should match. + const std::vector &JxW_ext = fe_values_ext.get_JxW_values (); const std::vector > &normals_int = fe_values_int.get_normal_vectors (); // AD variable @@ -387,11 +392,16 @@ void DGWeak::assemble_face_term_implicit( const unsigned int istate = fe_values_int.get_fe().system_to_component_index(itest_int).first; for (unsigned int iquad=0; iquad::assemble_face_term_implicit( for (unsigned int iquad=0; iquad::assemble_face_term_explicit( AssertDimension (n_dofs_ext, dof_indices_ext.size()); // Jacobian and normal should always be consistent between two elements - // even for non-conforming meshes? - const std::vector &JxW_int = fe_values_int.get_JxW_values (); + // In the case of the non-conforming mesh, we should be using the Jacobian + // of the smaller face since it would be "half" of the larger one. + // However, their curvature should match. + const std::vector &JxW_ext = fe_values_ext.get_JxW_values (); const std::vector > &normals_int = fe_values_int.get_normal_vectors (); // AD variable @@ -734,10 +746,10 @@ void DGWeak::assemble_face_term_explicit( for (unsigned int iquad=0; iquad::assemble_face_term_explicit( for (unsigned int iquad=0; iquad +namespace dealii{ +namespace numbers{ + bool is_finite(const Sacado::Fad::DFad &x) { + (void) x; + return true; + } +} +} + +#include +#include + + +int main (int /*argc*/, char * /*argv*/[]) +{ + using namespace dealii; + dealii::LinearAlgebra::distributed::Vector vector_double; + + using ADtype = Sacado::Fad::DFad; + dealii::LinearAlgebra::distributed::Vector vector_ad; + + return 0; +} + diff --git a/src/integral_hypercube.m b/src/dummy/integral_hypercube.m similarity index 100% rename from src/integral_hypercube.m rename to src/dummy/integral_hypercube.m diff --git a/src/template_instantiator.h b/src/dummy/template_instantiator.h similarity index 100% rename from src/template_instantiator.h rename to src/dummy/template_instantiator.h diff --git a/src/temporary_test.cpp b/src/dummy/temporary_test.cpp similarity index 100% rename from src/temporary_test.cpp rename to src/dummy/temporary_test.cpp diff --git a/src/test_manifold.cpp b/src/dummy/test_manifold.cpp similarity index 95% rename from src/test_manifold.cpp rename to src/dummy/test_manifold.cpp index 6b275507d..a64580fd1 100644 --- a/src/test_manifold.cpp +++ b/src/dummy/test_manifold.cpp @@ -48,12 +48,9 @@ namespace Step10 // Description of the upper face of the square class UpperManifold: public ChartManifold<2,2,2> { public: - virtual Point<2> - pull_back(const Point<2> &space_point) const; - virtual Point<2> - push_forward(const Point<2> &chart_point) const; - - virtual std::unique_ptr > clone() const; + virtual Point<2> pull_back(const Point<2> &space_point) const; ///< See dealii::Manifold::pull_back(). + virtual Point<2> push_forward(const Point<2> &chart_point) const; ///< See dealii::Manifold::pull_forward(). + virtual std::unique_ptr > clone() const; ///< See dealii::Manifold::clone(). }; Point<2> UpperManifold::pull_back(const Point<2> &space_point) const { diff --git a/src/linear_solver/CMakeLists.txt b/src/linear_solver/CMakeLists.txt index 9d6663511..54e9e7bfc 100644 --- a/src/linear_solver/CMakeLists.txt +++ b/src/linear_solver/CMakeLists.txt @@ -9,7 +9,9 @@ add_library(${LinearSolverLib} STATIC ${SOURCE}) target_compile_definitions(${LinearSolverLib} PRIVATE PHILIP_DIM=${dim}) # Setup target with deal.II -DEAL_II_SETUP_TARGET(${LinearSolverLib}) +if(NOT DOC_ONLY) + DEAL_II_SETUP_TARGET(${LinearSolverLib}) +endif() unset(LinearSolverLib) diff --git a/src/main.cpp b/src/main.cpp index ce244f178..3f1b2666f 100644 --- a/src/main.cpp +++ b/src/main.cpp @@ -14,9 +14,9 @@ int main (int argc, char *argv[]) { -#if !defined(__APPLE__) - feenableexcept(FE_INVALID | FE_OVERFLOW); // catch nan -#endif +//#if !defined(__APPLE__) +// feenableexcept(FE_INVALID | FE_OVERFLOW); // catch nan +//#endif dealii::deallog.depth_console(99); dealii::Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, 1); const int n_mpi = dealii::Utilities::MPI::n_mpi_processes(MPI_COMM_WORLD); diff --git a/src/numerical_flux/CMakeLists.txt b/src/numerical_flux/CMakeLists.txt index f105464d2..027aaca00 100644 --- a/src/numerical_flux/CMakeLists.txt +++ b/src/numerical_flux/CMakeLists.txt @@ -15,7 +15,9 @@ foreach(dim RANGE 1 3) string(CONCAT PhysicsLib Physics_${dim}D) target_link_libraries(${NumericalFluxLib} ${PhysicsLib}) # Setup target with deal.II - DEAL_II_SETUP_TARGET(${NumericalFluxLib}) + if(NOT DOC_ONLY) + DEAL_II_SETUP_TARGET(${NumericalFluxLib}) + endif() unset(NumericalFluxLib) diff --git a/src/numerical_flux/viscous_numerical_flux.h b/src/numerical_flux/viscous_numerical_flux.h index 717c48cdd..69627c0b1 100644 --- a/src/numerical_flux/viscous_numerical_flux.h +++ b/src/numerical_flux/viscous_numerical_flux.h @@ -30,10 +30,10 @@ virtual std::array evaluate_auxiliary_flux ( const real &penalty, const bool on_boundary = false) const = 0; -dealii::Tensor<1,dim, dealii::Tensor<1,nstate,real>> diffusion_matrix_int; -dealii::Tensor<1,dim, dealii::Tensor<1,nstate,real>> diffusion_matrix_int_transpose; -dealii::Tensor<1,dim, dealii::Tensor<1,nstate,real>> diffusion_matrix_ext; -dealii::Tensor<1,dim, dealii::Tensor<1,nstate,real>> diffusion_matrix_ext_transpose; +// dealii::Tensor<1,dim, dealii::Tensor<1,nstate,real>> diffusion_matrix_int; +// dealii::Tensor<1,dim, dealii::Tensor<1,nstate,real>> diffusion_matrix_int_transpose; +// dealii::Tensor<1,dim, dealii::Tensor<1,nstate,real>> diffusion_matrix_ext; +// dealii::Tensor<1,dim, dealii::Tensor<1,nstate,real>> diffusion_matrix_ext_transpose; }; @@ -74,7 +74,7 @@ std::array evaluate_auxiliary_flux ( const bool on_boundary = false) const; protected: -const std::shared_ptr < Physics::PhysicsBase > pde_physics; +const std::shared_ptr < Physics::PhysicsBase > pde_physics; ///< Associated physics. }; diff --git a/src/ode_solver/CMakeLists.txt b/src/ode_solver/CMakeLists.txt index 07a4ee3b0..19f03fb44 100644 --- a/src/ode_solver/CMakeLists.txt +++ b/src/ode_solver/CMakeLists.txt @@ -15,7 +15,9 @@ foreach(dim RANGE 1 3) target_link_libraries(${ODESolverLib} ${DiscontinuousGalerkinLib}) target_link_libraries(${ODESolverLib} ${LinearSolverLib}) # Setup target with deal.II - DEAL_II_SETUP_TARGET(${ODESolverLib}) + if(NOT DOC_ONLY) + DEAL_II_SETUP_TARGET(${ODESolverLib}) + endif() unset(ODESolverLib) diff --git a/src/ode_solver/ode_solver.cpp b/src/ode_solver/ode_solver.cpp index 3e07e4ab9..efada4930 100644 --- a/src/ode_solver/ode_solver.cpp +++ b/src/ode_solver/ode_solver.cpp @@ -61,6 +61,7 @@ template int ODESolver::steady_state () { Parameters::ODESolverParam ode_param = ODESolver::all_parameters->ode_solver_param; + if (ode_param.output_solution_every_x_steps >= 0) this->dg->output_results_vtk(this->current_iteration); pcout << " Performing steady state analysis... " << std::endl; allocate_ode_system (); @@ -69,8 +70,6 @@ int ODESolver::steady_state () update_norm = 1; // Always do at least 1 iteration this->current_iteration = 0; - if (ode_param.output_solution_every_x_steps >= 0) this->dg->output_results_vtk(this->current_iteration); - pcout << " Evaluating right-hand side and setting system_matrix to Jacobian before starting iterations... " << std::endl; this->dg->assemble_residual (); initial_residual_norm = this->dg->get_residual_l2norm(); diff --git a/src/ode_solver/ode_solver.h b/src/ode_solver/ode_solver.h index 61398ba2d..11df52d05 100644 --- a/src/ode_solver/ode_solver.h +++ b/src/ode_solver/ode_solver.h @@ -22,17 +22,24 @@ template class ODESolver { public: - ODESolver(int ode_solver_type); ///< Constructor - ODESolver(std::shared_ptr< DGBase > dg_input); ///< Constructor - virtual ~ODESolver() {}; ///< Destructor + ODESolver(int ode_solver_type); ///< Constructor. + ODESolver(std::shared_ptr< DGBase > dg_input); ///< Constructor. + virtual ~ODESolver() {}; ///< Destructor. /// Useful for accurate time-stepping. /** This variable will change when advance_solution_time() or step_in_time() is called. */ double current_time; - /// Evaluate steady state solution + /// Evaluate steady state solution. int steady_state (); + /// Ramps up the solution from p0 all the way up to the given global_final_poly_degree. + /** This first interpolates the current solution to the P0 space as an initial solution. + * The P0 is then solved, interpolated to P1, and the process is repeated until the + * desired polynomial is solved. + * + * This is mainly usely for grid studies. + */ void initialize_steady_polynomial_ramping (const unsigned int global_final_poly_degree); @@ -48,8 +55,8 @@ class ODESolver unsigned int current_iteration; ///< Current iteration. protected: - double update_norm; - double initial_residual_norm; + double update_norm; ///< Norm of the solution update. + double initial_residual_norm; ///< Initial residual norm. /// Virtual function to evaluate solution update virtual void step_in_time(real dt) = 0; @@ -61,15 +68,19 @@ class ODESolver /// Solution update given by the ODE solver dealii::LinearAlgebra::distributed::Vector solution_update; + /// Stores the various RK stages. + /** Currently hard-coded to RK4. + */ std::vector> rk_stage; /// Smart pointer to DGBase std::shared_ptr> dg; + /// Input parameters. const Parameters::AllParameters *const all_parameters; - const MPI_Comm mpi_communicator; - dealii::ConditionalOStream pcout; + const MPI_Comm mpi_communicator; ///< MPI communicator. + dealii::ConditionalOStream pcout; ///< Parallel std::cout that only outputs on mpi_rank==0 }; // end of ODESolver class @@ -104,10 +115,15 @@ class Implicit_ODESolver ODESolver::ODESolver(dg_input) {}; ~Implicit_ODESolver() {}; ///< Destructor. + /// Allocates ODE system based on given DGBase. + /** Basically allocates solution vector and asks DGBase to evaluate the mass matrix. + */ void allocate_ode_system (); protected: + /// Advances the solution in time by \p dt. void step_in_time(real dt); - using ODESolver::pcout; + + using ODESolver::pcout; ///< Parallel std::cout that only outputs on mpi_rank==0 }; // end of Implicit_ODESolver class @@ -121,16 +137,20 @@ class Explicit_ODESolver : public ODESolver { public: - Explicit_ODESolver() = delete; + /// Constructor. Explicit_ODESolver(std::shared_ptr> dg_input) : ODESolver::ODESolver(dg_input) {}; + /// Destructor. ~Explicit_ODESolver() {}; + /// Allocates ODE system based on given DGBase. + /** Basically allocates rk-stages and solution vector and asks DGBase to evaluate the inverse mass matrix. + */ void allocate_ode_system (); protected: - void step_in_time(real dt); - using ODESolver::pcout; + void step_in_time(real dt); ///< Advances the solution in time by \p dt. + using ODESolver::pcout; ///< Parallel std::cout that only outputs on mpi_rank==0 }; // end of Explicit_ODESolver class /// Creates and assemble Explicit_ODESolver or Implicit_ODESolver as ODESolver based on input. @@ -138,8 +158,11 @@ template class ODESolverFactory { public: + /// Creates the ODE solver given a DGBase. + /** The input parameters are copied from the DGBase since they should be consistent + */ static std::shared_ptr> create_ODESolver(std::shared_ptr< DGBase > dg_input); - static std::shared_ptr> create_ODESolver(Parameters::ODESolverParam::ODESolverEnum ode_solver_type); + // static std::shared_ptr> create_ODESolver(Parameters::ODESolverParam::ODESolverEnum ode_solver_type); }; diff --git a/src/parameters/CMakeLists.txt b/src/parameters/CMakeLists.txt index 167a19f8e..c1d68dffe 100644 --- a/src/parameters/CMakeLists.txt +++ b/src/parameters/CMakeLists.txt @@ -12,6 +12,8 @@ set(ParameterLib ParametersLibrary) add_library(${ParameterLib} STATIC ${SOURCE}) # Setup target with deal.II -DEAL_II_SETUP_TARGET(${ParameterLib}) +if(NOT DOC_ONLY) + DEAL_II_SETUP_TARGET(${ParameterLib}) +endif() unset(ParameterLib) diff --git a/src/parameters/all_parameters.h b/src/parameters/all_parameters.h index f66f11f9a..f21dfbcd8 100644 --- a/src/parameters/all_parameters.h +++ b/src/parameters/all_parameters.h @@ -41,6 +41,9 @@ class AllParameters /// Flag to use split form. bool use_split_form; + /// Flag to use periodic BC. + /** Not fully tested. + */ bool use_periodic_bc; /// Number of state variables. Will depend on PDE @@ -60,7 +63,7 @@ class AllParameters burgers_split_form, advection_periodicity, }; - TestType test_type; + TestType test_type; ///< Selected TestType from the input file. /// Currently allows to solve advection, diffusion, convection-diffusion enum PartialDifferentialEquation { @@ -115,7 +118,7 @@ class AllParameters //FunctionParser initial_conditions; //BoundaryConditions boundary_conditions[max_n_boundaries]; protected: - dealii::ConditionalOStream pcout; + dealii::ConditionalOStream pcout; ///< Parallel std::cout that only outputs on mpi_rank==0 }; diff --git a/src/parameters/parameters_euler.h b/src/parameters/parameters_euler.h index d6dcc08f3..c2ed5e010 100644 --- a/src/parameters/parameters_euler.h +++ b/src/parameters/parameters_euler.h @@ -9,9 +9,9 @@ namespace Parameters { class EulerParam { public: - double ref_length; - double mach_inf; - double gamma_gas; + double ref_length; ///< Reference length. + double mach_inf; ///< Mach number at infinity. + double gamma_gas; ///< Adiabatic index of the fluid. /// Input file provides in degrees, but the value stored here is in radians double angle_of_attack; /// Input file provides in degrees, but the value stored here is in radians diff --git a/src/parameters/parameters_linear_solver.h b/src/parameters/parameters_linear_solver.h index ee4e50545..627ede01c 100644 --- a/src/parameters/parameters_linear_solver.h +++ b/src/parameters/parameters_linear_solver.h @@ -13,24 +13,30 @@ namespace Parameters { class LinearSolverParam { public: - LinearSolverParam (); ///< Constructor - enum LinearSolverEnum { direct, gmres }; ///< Types of linear solvers available + LinearSolverParam (); ///< Constructor. + /// Types of linear solvers available. + enum LinearSolverEnum { + direct, /// LU. + gmres /// GMRES. + }; /// Can either be verbose or quiet. /** Verbose will print the full dense matrix. Will not work for large matrices */ - OutputEnum linear_solver_output; - LinearSolverEnum linear_solver_type; + OutputEnum linear_solver_output; ///< quiet or verbose. + LinearSolverEnum linear_solver_type; ///< direct or gmres. // GMRES options double ilut_drop; ///< Threshold to drop terms close to zero. + //@{ /// Add to the diagonal. - /** For some reason it helps some problems. + /** Based on deal.II doc it helps some problems? * From what I have seen, it doesn't help ours. */ double ilut_rtol, ilut_atol; - /// ILU fill-in - int ilut_fill; + //@} + + int ilut_fill; ///< ILU fill-in double linear_residual; ///< Tolerance for linear residual. int max_iterations; ///< Maximum number of linear iteration. diff --git a/src/parameters/parameters_ode_solver.h b/src/parameters/parameters_ode_solver.h index 8910e379a..37d32685a 100644 --- a/src/parameters/parameters_ode_solver.h +++ b/src/parameters/parameters_ode_solver.h @@ -12,7 +12,11 @@ class ODESolverParam { public: ODESolverParam (); ///< Constructor. - enum ODESolverEnum { explicit_solver, implicit_solver }; ///< Types of ODE solver. Will determine memory allocation. + /// Types of ODE solver + enum ODESolverEnum { + explicit_solver, /// RK4 + implicit_solver /// Backward-Euler + }; OutputEnum ode_output; ///< verbose or quiet. ODESolverEnum ode_solver_type; ///< ODE solver type. Note that only implicit has been fully tested for now. @@ -28,10 +32,8 @@ class ODESolverParam double time_step_factor_residual; ///< Multiplies initial time-step by time_step_factor_residual*(-log10(residual_norm_decrease)) double time_step_factor_residual_exp; ///< Scales initial time step by pow(time_step_factor_residual*(-log10(residual_norm_decrease)),time_step_factor_residual_exp) - /// Declares the possible variables and sets the defaults. - static void declare_parameters (dealii::ParameterHandler &prm); - /// Parses input file and sets the variables. - void parse_parameters (dealii::ParameterHandler &prm); + static void declare_parameters (dealii::ParameterHandler &prm); ///< Declares the possible variables and sets the defaults. + void parse_parameters (dealii::ParameterHandler &prm); ///< Parses input file and sets the variables. }; } // Parameters namespace diff --git a/src/physics/CMakeLists.txt b/src/physics/CMakeLists.txt index b3a588035..d2758fd21 100644 --- a/src/physics/CMakeLists.txt +++ b/src/physics/CMakeLists.txt @@ -14,7 +14,9 @@ foreach(dim RANGE 1 3) add_library(${PhysicsLib} STATIC ${SOURCE}) # Setup target with deal.II - DEAL_II_SETUP_TARGET(${PhysicsLib}) + if(NOT DOC_ONLY) + DEAL_II_SETUP_TARGET(${PhysicsLib}) + endif() target_compile_definitions(${PhysicsLib} PRIVATE PHILIP_DIM=${dim}) diff --git a/src/physics/convection_diffusion.cpp b/src/physics/convection_diffusion.cpp index 0a464fb52..03c5180bd 100644 --- a/src/physics/convection_diffusion.cpp +++ b/src/physics/convection_diffusion.cpp @@ -76,6 +76,19 @@ ::convective_flux (const std::array &solution) const return conv_flux; } +template +std::array,nstate> ConvectionDiffusion +::convective_numerical_split_flux ( + const std::array &soln1, + const std::array &soln2) const +{ + std::array arr_avg; + for (int i = 0 ; i < nstate; ++i) { + arr_avg[i] = (soln1[i] + soln2[i])/2.; + } + return convective_flux(arr_avg); +} + template dealii::Tensor<1,dim,real> ConvectionDiffusion ::advection_speed () const diff --git a/src/physics/convection_diffusion.h b/src/physics/convection_diffusion.h index 02c7e4ae3..1ca7fbefc 100644 --- a/src/physics/convection_diffusion.h +++ b/src/physics/convection_diffusion.h @@ -28,8 +28,10 @@ template class ConvectionDiffusion : public PhysicsBase { public: - const bool hasConvection; - const bool hasDiffusion; + const bool hasConvection; ///< Turns ON/OFF convection term. + + const bool hasDiffusion; ///< Turns ON/OFF diffusion term. + /// Constructor ConvectionDiffusion (const bool convection = true, const bool diffusion = true) : hasConvection(convection), hasDiffusion(diffusion) @@ -42,17 +44,9 @@ class ConvectionDiffusion : public PhysicsBase /// Convective flux: \f$ \mathbf{F}_{conv} = u \f$ std::array,nstate> convective_flux (const std::array &solution) const; - std::array,nstate> convective_numerical_split_flux ( - const std::array &soln_const, const std::array &soln_loop) const - { - std::array arr_avg; - for (int i = 0 ; i < nstate; ++i) - { - arr_avg[i] = (soln_const[i] + soln_loop[i])/2.; - } - return convective_flux(arr_avg); - }; + const std::array &soln1, + const std::array &soln2) const; /// Spectral radius of convective term Jacobian is 'c' std::array convective_eigenvalues ( diff --git a/src/physics/euler.cpp b/src/physics/euler.cpp index c078cbccf..498e72761 100644 --- a/src/physics/euler.cpp +++ b/src/physics/euler.cpp @@ -257,40 +257,68 @@ ::compute_mach_number ( const std::array &conservative_soln ) const // Split form functions: +template +std::array,nstate> Euler +::convective_numerical_split_flux(const std::array &conservative_soln1, + const std::array &conservative_soln2) const +{ + std::array,nstate> conv_num_split_flux; + const real mean_density = compute_mean_density(conservative_soln1, conservative_soln2); + const real mean_pressure = compute_mean_pressure(conservative_soln1, conservative_soln2); + const dealii::Tensor<1,dim,real> mean_velocities = compute_mean_velocities(conservative_soln1,conservative_soln2); + const real mean_specific_energy = compute_mean_specific_energy(conservative_soln1, conservative_soln2); + + for (int flux_dim = 0; flux_dim < dim; ++flux_dim) + { + // Density equation + conv_num_split_flux[0][flux_dim] = mean_density * mean_velocities[flux_dim];//conservative_soln[1+flux_dim]; + // Momentum equation + for (int velocity_dim=0; velocity_dim inline real Euler:: -compute_mean_density(const std::array &soln_const, - const std::array &soln_loop) const +compute_mean_density(const std::array &conservative_soln1, + const std::array &conservative_soln2) const { - return (soln_const[0] + soln_loop[0])/2.; + return (conservative_soln1[0] + conservative_soln2[0])/2.; } template inline real Euler:: -compute_mean_pressure(const std::array &soln_const, - const std::array &soln_loop) const +compute_mean_pressure(const std::array &conservative_soln1, + const std::array &conservative_soln2) const { - real pressure_const = compute_pressure(soln_const); - real pressure_loop = compute_pressure(soln_loop); - return (pressure_const + pressure_loop)/2.; + real pressure_1 = compute_pressure(conservative_soln1); + real pressure_2 = compute_pressure(conservative_soln2); + return (pressure_1 + pressure_2)/2.; } template inline dealii::Tensor<1,dim,real> Euler:: -compute_mean_velocities(const std::array &soln_const, - const std::array &soln_loop) const +compute_mean_velocities(const std::array &conservative_soln1, + const std::array &conservative_soln2) const { - dealii::Tensor<1,dim,real> vel_const = compute_velocities(soln_const); - dealii::Tensor<1,dim,real> vel_loop = compute_velocities(soln_loop); - return (vel_const + vel_loop)/2.; + dealii::Tensor<1,dim,real> vel_1 = compute_velocities(conservative_soln1); + dealii::Tensor<1,dim,real> vel_2 = compute_velocities(conservative_soln2); + return (vel_1 + vel_2)/2.; } template inline real Euler:: -compute_mean_specific_energy(const std::array &soln_const, - const std::array &soln_loop) const +compute_mean_specific_energy(const std::array &conservative_soln1, + const std::array &conservative_soln2) const { - return ((soln_const[nstate-1]/soln_const[0]) + (soln_loop[nstate-1]/soln_loop[0]))/2.; + return ((conservative_soln1[nstate-1]/conservative_soln1[0]) + (conservative_soln2[nstate-1]/conservative_soln2[0]))/2.; } @@ -319,33 +347,6 @@ ::convective_flux (const std::array &conservative_soln) const return conv_flux; } -template -std::array,nstate> Euler -::convective_numerical_split_flux(const std::array &soln_const, - const std::array &soln_loop) const -{ - std::array,nstate> conv_num_split_flux; - const real mean_density = compute_mean_density(soln_const, soln_loop); - const real mean_pressure = compute_mean_pressure(soln_const, soln_loop); - const dealii::Tensor<1,dim,real> mean_velocities = compute_mean_velocities(soln_const,soln_loop); - const real mean_specific_energy = compute_mean_specific_energy(soln_const, soln_loop); - - for (int flux_dim = 0; flux_dim < dim; ++flux_dim) - { - // Density equation - conv_num_split_flux[0][flux_dim] = mean_density * mean_velocities[flux_dim];//conservative_soln[1+flux_dim]; - // Momentum equation - for (int velocity_dim=0; velocity_dim std::array Euler ::convective_normal_flux (const std::array &conservative_soln, const dealii::Tensor<1,dim,real> &normal) const diff --git a/src/physics/euler.h b/src/physics/euler.h index 5a6a3e0e1..ec52411e5 100644 --- a/src/physics/euler.h +++ b/src/physics/euler.h @@ -119,46 +119,49 @@ class Euler : public PhysicsBase ~Euler () {}; - const double ref_length; - /// Constant heat capacity ratio of air - const double gam; - /// Gamma-1.0 used often - const double gamm1; + const double ref_length; ///< Reference length. + const double gam; ///< Constant heat capacity ratio of fluid. + const double gamm1; ///< Constant heat capacity ratio (Gamma-1.0) used often. /// Non-dimensionalized density* at infinity. density* = density/density_ref /// Choose density_ref = density(inf) /// density*(inf) = density(inf) / density_ref = density(inf)/density(inf) = 1.0 const double density_inf; - const double mach_inf; - const double mach_inf_sqr; + const double mach_inf; ///< Farfield Mach number. + const double mach_inf_sqr; ///< Farfield Mach number squared. + /// Angle of attack. + /** Mandatory for 2D simulations. + */ const double angle_of_attack; + /// Sideslip angle. + /** Mandatory for 2D and 3D simulations. + */ const double side_slip_angle; - const double sound_inf; /// Non-dimensionalized sound* at infinity - const double pressure_inf; /// Non-dimensionalized pressure* at infinity - const double entropy_inf; /// Entropy measure at infinity - double temperature_inf; /// Non-dimensionalized temperature* at infinity. Should equal 1/density*(inf) + const double sound_inf; ///< Non-dimensionalized sound* at infinity + const double pressure_inf; ///< Non-dimensionalized pressure* at infinity + const double entropy_inf; ///< Entropy measure at infinity + double temperature_inf; ///< Non-dimensionalized temperature* at infinity. Should equal 1/density*(inf) //const double internal_energy_inf; + /// Non-dimensionalized Velocity vector at farfield + /** Evaluated using mach_number, angle_of_attack, and side_slip_angle. + */ dealii::Tensor<1,dim,double> velocities_inf; // should be const - dealii::Tensor<1,dim,double> compute_velocities_inf() const; - + // dealii::Tensor<1,dim,double> compute_velocities_inf() const; - - std::array manufactured_solution (const dealii::Point &pos) const; + // std::array manufactured_solution (const dealii::Point &pos) const; /// Convective flux: \f$ \mathbf{F}_{conv} \f$ std::array,nstate> convective_flux ( const std::array &conservative_soln) const; - std::array,nstate> convective_numerical_split_flux ( - const std::array &soln_const, const std::array &soln_loop) const; - + /// Convective normal flux: \f$ \mathbf{F}_{conv} \cdot \hat{n} \f$ std::array convective_normal_flux (const std::array &conservative_soln, const dealii::Tensor<1,dim,real> &normal) const; /// Convective flux Jacobian: \f$ \frac{\partial \mathbf{F}_{conv}}{\partial w} \cdot \mathbf{n} \f$ @@ -218,6 +221,9 @@ class Euler : public PhysicsBase /// Given primitive variables, returns velocities. dealii::Tensor<1,dim,real> extract_velocities_from_primitive ( const std::array &primitive_soln ) const; /// Given primitive variables, returns total energy + /** @param[in] primitive_soln Primitive solution (density, momentum, energy) + * \return Entropy measure + */ real compute_total_energy ( const std::array &primitive_soln ) const; /// Evaluate entropy from conservative variables @@ -225,6 +231,9 @@ class Euler : public PhysicsBase * Used to check entropy convergence * See discussion in * https://physics.stackexchange.com/questions/116779/entropy-is-constant-how-to-express-this-equation-in-terms-of-pressure-and-densi?answertab=votes#tab-top + * + * @param[in] conservative_soln Conservative solution (density, momentum, energy) + * \return Entropy measure */ real compute_entropy_measure ( const std::array &conservative_soln ) const; @@ -246,19 +255,39 @@ class Euler : public PhysicsBase /** See the book I do like CFD, sec 4.14.2 */ real compute_temperature_from_density_pressure ( const real density, const real pressure ) const; - ///These functions are only relevant to the split form. The Euler split form is that of Kennedy & Gruber. - /** Refer to Gassner's paper for more information: */ - real compute_mean_density(const std::array &soln_const, - const std::array &soln_loop) const; + /// The Euler split form is that of Kennedy & Gruber. + /** Refer to Gassner's paper (2016) Eq. 3.10 for more information: */ + std::array,nstate> convective_numerical_split_flux ( + const std::array &conservative_soln1, + const std::array &conservative_soln2) const; - real compute_mean_pressure(const std::array &soln_const, - const std::array &soln_loop) const; + /// Mean density given two sets of conservative solutions. + /** Used in the implementation of the split form. + */ + real compute_mean_density( + const std::array &conservative_soln1, + const std::array &convervative_soln2) const; - dealii::Tensor<1,dim,real> compute_mean_velocities(const std::array &soln_const, - const std::array &soln_loop) const; + /// Mean pressure given two sets of conservative solutions. + /** Used in the implementation of the split form. + */ + real compute_mean_pressure( + const std::array &conservative_soln1, + const std::array &convervative_soln2) const; - real compute_mean_specific_energy(const std::array &soln_const, - const std::array &soln_loop) const; + /// Mean velocities given two sets of conservative solutions. + /** Used in the implementation of the split form. + */ + dealii::Tensor<1,dim,real> compute_mean_velocities( + const std::array &conservative_soln1, + const std::array &convervative_soln2) const; + + /// Mean specific energy given two sets of conservative solutions. + /** Used in the implementation of the split form. + */ + real compute_mean_specific_energy( + const std::array &conservative_soln1, + const std::array &convervative_soln2) const; void boundary_face_values ( const int /*boundary_type*/, diff --git a/src/physics/manufactured_solution.h b/src/physics/manufactured_solution.h index 7bc048756..03e7d9a10 100644 --- a/src/physics/manufactured_solution.h +++ b/src/physics/manufactured_solution.h @@ -15,6 +15,9 @@ namespace PHiLiP { +/// Manufactured solution used for grid studies to check convergence orders. +/** This class also provides derivatives necessary to evaluate source terms. + */ template class ManufacturedSolutionFunction : public dealii::Function { @@ -70,6 +73,7 @@ class ManufacturedSolutionFunction : public dealii::Function /// Uses finite-difference to evaluate the hessian dealii::SymmetricTensor<2,dim,real> hessian_fd (const dealii::Point &point, const unsigned int istate = 0) const; + /// Same as Function::values() except it returns it into a std::vector format. std::vector stdvector_values (const dealii::Point &point) const; // Virtual functions inherited from dealii::Function @@ -107,9 +111,13 @@ class ManufacturedSolutionFunction : public dealii::Function // std::vector > > &gradients) const; private: + //@{ + /** Constants used to manufactured solution. + */ std::vector base_values; std::vector amplitudes; std::vector> frequencies; + //@} }; } diff --git a/src/physics/mhd.h b/src/physics/mhd.h index 7236b78c0..7e45f0904 100644 --- a/src/physics/mhd.h +++ b/src/physics/mhd.h @@ -95,16 +95,13 @@ class MHD : public PhysicsBase // double mach_inf_sqr = 1; - std::array manufactured_solution (const dealii::Point &pos) const; + //std::array manufactured_solution (const dealii::Point &pos) const; /// Convective flux: \f$ \mathbf{F}_{conv} \f$ std::array,nstate> convective_flux ( const std::array &conservative_soln) const; - std::array,nstate> convective_numerical_split_flux ( - const std::array &soln_const, const std::array &soln_loop) const; - - + /// Convective flux: \f$ \mathbf{F}_{conv} \hat{n} \f$ std::array convective_normal_flux (const std::array &conservative_soln, const dealii::Tensor<1,dim,real> &normal) const; /// Convective flux Jacobian: \f$ \frac{\partial \mathbf{F}_{conv}}{\partial w} \cdot \mathbf{n} \f$ @@ -195,19 +192,39 @@ class MHD : public PhysicsBase /** See the book I do like CFD, sec 4.14.2 */ real compute_temperature_from_density_pressure ( const real density, const real pressure ) const; - ///These functions are only relevant to the split form. The Euler split form is that of Kennedy & Gruber. - /** Refer to Gassner's paper for more information: */ - real compute_mean_density(const std::array &soln_const, - const std::array &soln_loop) const; + /// The Euler split form is that of Kennedy & Gruber. + /** Refer to Gassner's paper (2016) Eq. 3.10 for more information: */ + std::array,nstate> convective_numerical_split_flux ( + const std::array &conservative_soln1, + const std::array &conservative_soln2) const; - real compute_mean_pressure(const std::array &soln_const, - const std::array &soln_loop) const; + /// Mean density given two sets of conservative solutions. + /** Used in the implementation of the split form. + */ + real compute_mean_density( + const std::array &conservative_soln1, + const std::array &convervative_soln2) const; - dealii::Tensor<1,dim,real> compute_mean_velocities(const std::array &soln_const, - const std::array &soln_loop) const; + /// Mean pressure given two sets of conservative solutions. + /** Used in the implementation of the split form. + */ + real compute_mean_pressure( + const std::array &conservative_soln1, + const std::array &convervative_soln2) const; + + /// Mean velocities given two sets of conservative solutions. + /** Used in the implementation of the split form. + */ + dealii::Tensor<1,dim,real> compute_mean_velocities( + const std::array &conservative_soln1, + const std::array &convervative_soln2) const; - real compute_mean_specific_energy(const std::array &soln_const, - const std::array &soln_loop) const; + /// Mean specific energy given two sets of conservative solutions. + /** Used in the implementation of the split form. + */ + real compute_mean_specific_energy( + const std::array &conservative_soln1, + const std::array &convervative_soln2) const; // void boundary_face_values ( // const int /*boundary_type*/, diff --git a/src/physics/physics.h b/src/physics/physics.h index 9470735ff..72868ff27 100644 --- a/src/physics/physics.h +++ b/src/physics/physics.h @@ -49,7 +49,8 @@ class PhysicsBase const std::array &soln_const, const std::array &soln_loop) const = 0; /// Spectral radius of convective term Jacobian. - /// Used for scalar dissipation + /** Used for scalar dissipation + */ virtual std::array convective_eigenvalues ( const std::array &/*solution*/, const dealii::Tensor<1,dim,real> &/*normal*/) const = 0; @@ -63,7 +64,8 @@ class PhysicsBase // const std::array,nstate> &solution_grad) const = 0; /// Dissipative fluxes that will be differentiated once in space. - /// Evaluates the dissipative flux through the linearization F = A(u)*grad(u). + /** Evaluates the dissipative flux through the linearization F = A(u)*grad(u). + */ std::array,nstate> dissipative_flux_A_gradu ( const real scaling, const std::array &solution, @@ -90,27 +92,45 @@ class PhysicsBase std::array &/*soln_bc*/, std::array,nstate> &/*soln_grad_bc*/) const; + /// Returns current vector solution to be used by PhysicsPostprocessor to output current solution. + /** The implementation in this Physics base class simply returns the stored solution. + */ virtual dealii::Vector post_compute_derived_quantities_vector ( - const dealii::Vector &uh, - const std::vector > &duh, - const std::vector > &dduh, - const dealii::Tensor<1,dim> &normals, - const dealii::Point &evaluation_points) const; - + const dealii::Vector &uh, + const std::vector > &/*duh*/, + const std::vector > &/*dduh*/, + const dealii::Tensor<1,dim> &/*normals*/, + const dealii::Point &/*evaluation_points*/) const; + + /// Returns current scalar solution to be used by PhysicsPostprocessor to output current solution. + /** The implementation in this Physics base class simply returns the stored solution. + */ virtual dealii::Vector post_compute_derived_quantities_scalar ( const double &uh, const dealii::Tensor<1,dim> &/*duh*/, const dealii::Tensor<2,dim> &/*dduh*/, const dealii::Tensor<1,dim> &/*normals*/, const dealii::Point &/*evaluation_points*/) const; + /// Returns names of the solution to be used by PhysicsPostprocessor to output current solution. + /** The implementation in this Physics base class simply returns "state0, state1, etc.". + */ virtual std::vector post_get_names () const; + /// Returns DataComponentInterpretation of the solution to be used by PhysicsPostprocessor to output current solution. + /** Treats every solution state as an independent scalar. + */ virtual std::vector post_get_data_component_interpretation () const; + /// Returns required update flags of the solution to be used by PhysicsPostprocessor to output current solution. + /** Only update the solution at the output points. + */ virtual dealii::UpdateFlags post_get_needed_update_flags () const; protected: - /// Some constants used to define manufactured solution + //@{ + /** Velocity constants used to define convection-diffusion type of manufactured solution + */ double velo_x, velo_y, velo_z; - double diff_coeff; + //@} + double diff_coeff; ///< Diffusion constant used to define convection-diffusion type of manufactured solution /// Anisotropic diffusion matrix /** As long as the diagonal components are positive and diagonally dominant diff --git a/src/physics/physics_factory.h b/src/physics/physics_factory.h index 30185118d..1fa10596f 100644 --- a/src/physics/physics_factory.h +++ b/src/physics/physics_factory.h @@ -13,6 +13,7 @@ template class PhysicsFactory { public: + /// Factory to return the correct physics given input file. static std::shared_ptr< PhysicsBase > create_Physics(const Parameters::AllParameters *const parameters_input); }; diff --git a/src/post_processor/CMakeLists.txt b/src/post_processor/CMakeLists.txt index cb69fa72b..1952c8545 100644 --- a/src/post_processor/CMakeLists.txt +++ b/src/post_processor/CMakeLists.txt @@ -8,7 +8,9 @@ foreach(dim RANGE 1 3) add_library(${PostprocessingLib} STATIC ${SOURCE}) # Setup target with deal.II - DEAL_II_SETUP_TARGET(${PostprocessingLib}) + if(NOT DOC_ONLY) + DEAL_II_SETUP_TARGET(${PostprocessingLib}) + endif() target_compile_definitions(${PostprocessingLib} PRIVATE PHILIP_DIM=${dim}) # target_link_libraries(${PostprocessingLib} m mpi mpi_cxx) diff --git a/src/post_processor/physics_post_processor.h b/src/post_processor/physics_post_processor.h index 599170390..f2ba40738 100644 --- a/src/post_processor/physics_post_processor.h +++ b/src/post_processor/physics_post_processor.h @@ -8,22 +8,26 @@ namespace PHiLiP { namespace Postprocess { -/// Create the post-processor with the correct template parameter +/// Postprocessor factory. template class PostprocessorFactory { public: + /// Create the post-processor with the correct template parameters. static std::unique_ptr< dealii::DataPostprocessor > create_Postprocessor(const Parameters::AllParameters *const parameters_input); }; +/// Postprocessor to output solution and other values computed by associated physics. +/** The functions in this class will call the Physics functions to query data. + */ template class PhysicsPostprocessor : public dealii::DataPostprocessor { public: - /// Constructor + /// Constructor. PhysicsPostprocessor (const Parameters::AllParameters *const parameters_input); - /// Physics that the post-processor will use to evaluate derived data types + /// Physics that the post-processor will use to evaluate derived data types. std::shared_ptr < Physics::PhysicsBase > physics; /// Queries the Physics to output data of a vector-valued problem. diff --git a/src/testing/CMakeLists.txt b/src/testing/CMakeLists.txt index ca9a7bd3b..309908508 100644 --- a/src/testing/CMakeLists.txt +++ b/src/testing/CMakeLists.txt @@ -27,7 +27,9 @@ foreach(dim RANGE 1 3) target_link_libraries(${TestsLib} ${DiscontinuousGalerkinLib}) target_link_libraries(${TestsLib} ${ODESolverLib}) # Setup target with deal.II - DEAL_II_SETUP_TARGET(${TestsLib}) + if(NOT DOC_ONLY) + DEAL_II_SETUP_TARGET(${TestsLib}) + endif() unset(TestsLib) unset(DiscontinuousGalerkinLib) diff --git a/src/testing/advection_explicit_periodic.cpp b/src/testing/advection_explicit_periodic.cpp index 229472561..00b04cc96 100644 --- a/src/testing/advection_explicit_periodic.cpp +++ b/src/testing/advection_explicit_periodic.cpp @@ -28,8 +28,6 @@ template AdvectionPeriodic::AdvectionPeriodic(const PHiLiP::Parameters::AllParameters *const parameters_input) : TestsBase::TestsBase(parameters_input) -, mpi_communicator(MPI_COMM_WORLD) -, pcout(std::cout, dealii::Utilities::MPI::this_mpi_process(mpi_communicator)==0) {} diff --git a/src/testing/advection_explicit_periodic.h b/src/testing/advection_explicit_periodic.h index c75766878..349523ecb 100644 --- a/src/testing/advection_explicit_periodic.h +++ b/src/testing/advection_explicit_periodic.h @@ -8,17 +8,20 @@ namespace PHiLiP { namespace Tests { +/// Linear advection through periodic boundary conditions. template class AdvectionPeriodic: public TestsBase { public: - AdvectionPeriodic() = delete; + /// Constructor. AdvectionPeriodic(const Parameters::AllParameters *const parameters_input); + + /// Currently passes no matter what. + /** Since it is linear advection, the exact solution about time T is known. Convergence orders can/should be checked. + * TO BE FIXED. + */ int run_test () const override; private: - //const Parameters::AllParameters *const all_parameters; ///< Pointer to all parameters - const MPI_Comm mpi_communicator; - dealii::ConditionalOStream pcout; }; } // End of Tests namespace diff --git a/src/testing/burgers_stability.h b/src/testing/burgers_stability.h index 4b3bfd703..98cde6d86 100644 --- a/src/testing/burgers_stability.h +++ b/src/testing/burgers_stability.h @@ -8,13 +8,25 @@ namespace PHiLiP { namespace Tests { +/// Burgers sine wave to shock. +/** Ensure that the kinetic energy is bounded. + * Gassner 2017. + */ template class BurgersEnergyStability: public TestsBase { public: + /// Constructor. BurgersEnergyStability(const Parameters::AllParameters *const parameters_input); + /// Ensure that the kinetic energy is bounded. + /** If the kinetic energy increases about its initial value, then the test should fail. + * Gassner 2017. + */ int run_test () const override; private: + /// Computes an integral of the kinetic energy (solution squared) in the entire domain. + /** Uses Inverse of inverse mass matrix (?) to evaluate integral of u^2. + */ double compute_energy(std::shared_ptr < PHiLiP::DGBase > &dg) const; }; diff --git a/src/testing/euler_cylinder.cpp b/src/testing/euler_cylinder.cpp index f3203a251..5380ba112 100644 --- a/src/testing/euler_cylinder.cpp +++ b/src/testing/euler_cylinder.cpp @@ -28,8 +28,6 @@ #include "dg/dg.h" #include "ode_solver/ode_solver.h" -//#include "template_instantiator.h" - namespace PHiLiP { namespace Tests { @@ -99,12 +97,17 @@ void half_cylinder(dealii::parallel::distributed::Triangulation<2> & tria, } +/// Function used to evaluate farfield conservative solution template class FreeStreamInitialConditions : public dealii::Function { public: - std::array far_field_conservative; + /// Farfield conservative solution + std::array farfield_conservative; + /// Constructor. + /** Evaluates the primary farfield solution and converts it into the store farfield_conservative solution + */ FreeStreamInitialConditions (const Physics::Euler euler_physics) : dealii::Function(nstate) { @@ -114,14 +117,13 @@ class FreeStreamInitialConditions : public dealii::Function primitive_boundary_values[0] = density_bc; for (int d=0;d &/*point*/, const unsigned int istate) const { - return far_field_conservative[istate]; + return farfield_conservative[istate]; } }; template class FreeStreamInitialConditions ; @@ -182,8 +184,8 @@ ::run_test () const // Generate grid and mapping dealii::parallel::distributed::Triangulation grid(this->mpi_communicator, typename dealii::Triangulation::MeshSmoothing( - dealii::Triangulation::smoothing_on_refinement | - dealii::Triangulation::smoothing_on_coarsening)); + dealii::Triangulation::MeshSmoothing::smoothing_on_refinement | + dealii::Triangulation::MeshSmoothing::smoothing_on_coarsening)); const unsigned int n_cells_circle = n_1d_cells[0]; const unsigned int n_cells_radial = 2*n_cells_circle; @@ -200,6 +202,7 @@ ::run_test () const std::shared_ptr> ode_solver = ODE::ODESolverFactory::create_ODESolver(dg); ode_solver->initialize_steady_polynomial_ramping(poly_degree); + dealii::Vector estimated_error_per_cell(grid.n_active_cells()); for (unsigned int igrid=0; igrid, dealii::hp::DoFHandler> solution_transfer(dg->dof_handler); solution_transfer.prepare_for_coarsening_and_refinement(old_solution); dg->high_order_grid.prepare_for_coarsening_and_refinement(); + grid.refine_global (1); + //dealii::GridRefinement::refine_and_coarsen_fixed_number(grid, + // estimated_error_per_cell, + // 0.3, + // 0.03); + //grid.execute_coarsening_and_refinement(); dg->high_order_grid.execute_coarsening_and_refinement(); dg->allocate_system (); dg->solution.zero_out_ghosts(); @@ -236,43 +245,56 @@ ::run_test () const // Solve the steady state problem ode_solver->steady_state(); + const auto mapping = (*(dg->high_order_grid.mapping_fe_field)); + dealii::hp::MappingCollection mapping_collection(mapping); + dealii::hp::FEValues fe_values_collection_volume (mapping_collection, dg->fe_collection, dg->volume_quadrature_collection, dealii::update_values | dealii::update_JxW_values); ///< FEValues of volume. // Overintegrate the error to make sure there is not integration error in the error estimate - int overintegrate = 10; - dealii::QGauss quad_extra(dg->max_degree+1+overintegrate); - dealii::FEValues fe_values_extra(*(dg->high_order_grid.mapping_fe_field), dg->fe_collection[poly_degree], quad_extra, - dealii::update_values | dealii::update_JxW_values | dealii::update_quadrature_points); - const unsigned int n_quad_pts = fe_values_extra.n_quadrature_points; + //int overintegrate = 0; + //dealii::QGauss quad_extra(dg->max_degree+1+overintegrate); + //dealii::FEValues fe_values_extra(*(dg->high_order_grid.mapping_fe_field), dg->fe_collection[poly_degree], quad_extra, + // dealii::update_values | dealii::update_JxW_values | dealii::update_quadrature_points); std::array soln_at_q; double l2error = 0; double area = 0; const double exact_area = (std::pow(outer_radius+inner_radius, 2.0) - std::pow(inner_radius,2.0))*dealii::numbers::PI / 2.0; - std::vector dofs_indices (fe_values_extra.dofs_per_cell); const double entropy_inf = euler_physics_double.entropy_inf; + estimated_error_per_cell.reinit(grid.n_active_cells()); // Integrate solution error and output error for (auto cell = dg->dof_handler.begin_active(); cell!=dg->dof_handler.end(); ++cell) { if (!cell->is_locally_owned()) continue; - fe_values_extra.reinit (cell); + const int i_fele = cell->active_fe_index(); + const int i_quad = i_fele; + const int i_mapp = 0; + + fe_values_collection_volume.reinit (cell, i_quad, i_mapp, i_fele); + const dealii::FEValues &fe_values_volume = fe_values_collection_volume.get_present_fe_values(); + + //fe_values_volume.reinit (cell); + std::vector dofs_indices (fe_values_volume.dofs_per_cell); cell->get_dof_indices (dofs_indices); - for (unsigned int iquad=0; iquadsolution[dofs_indices[idof]] * fe_values_extra.shape_value_component(idof, iquad, istate); + for (unsigned int idof=0; idofsolution[dofs_indices[idof]] * fe_values_volume.shape_value_component(idof, iquad, istate); } const double entropy = euler_physics_double.compute_entropy_measure(soln_at_q); const double uexact = entropy_inf; - l2error += pow(entropy - uexact, 2) * fe_values_extra.JxW(iquad); + cell_l2error += pow(entropy - uexact, 2) * fe_values_volume.JxW(iquad); + estimated_error_per_cell[cell->active_cell_index()] = cell_l2error; + l2error += cell_l2error; - area += fe_values_extra.JxW(iquad); + area += fe_values_volume.JxW(iquad); } } const double l2error_mpi_sum = std::sqrt(dealii::Utilities::MPI::sum(l2error, mpi_communicator)); diff --git a/src/testing/euler_entropy_waves.h b/src/testing/euler_entropy_waves.h index 2aefd507d..91c363e44 100644 --- a/src/testing/euler_entropy_waves.h +++ b/src/testing/euler_entropy_waves.h @@ -9,7 +9,9 @@ namespace PHiLiP { namespace Tests { -// Masatsuka2018 section 6.3 & section 7.13.3 +/// Exact entropy waves solution. +/** Masatsuka2018 section 6.3 & section 7.13.3 + */ template class EulerEntropyWavesFunction : public dealii::Function { @@ -21,10 +23,10 @@ class EulerEntropyWavesFunction : public dealii::Function */ EulerEntropyWavesFunction (const Physics::Euler euler_physics, const real dimensional_density_inf); - const Physics::Euler euler_physics; + const Physics::Euler euler_physics; ///< Euler physics. - const real dimensional_density_inf; - real Q_inf; + const real dimensional_density_inf; ///< Dimensional density at infinity. + real Q_inf; ///< Velocity at infinity. /// Destructor ~EulerEntropyWavesFunction() {}; @@ -54,6 +56,7 @@ class EulerEntropyWavesFunction : public dealii::Function real value (const dealii::Point &point, const unsigned int istate = 0) const; private: + /// Exact solution using the current \p time provided by the dealii::Function class dealii::Point<2> advected_location(const dealii::Point<2> old_location) const; }; @@ -70,8 +73,6 @@ class EulerEntropyWaves: public TestsBase */ EulerEntropyWaves(const Parameters::AllParameters *const parameters_input); - ~EulerEntropyWaves() {}; ///< Destructor. - /// Manufactured grid convergence /** Will run the a grid convergence test for various p * on multiple grids to determine the order of convergence. diff --git a/src/testing/euler_gaussian_bump.cpp b/src/testing/euler_gaussian_bump.cpp index bd36c9c68..963798b84 100644 --- a/src/testing/euler_gaussian_bump.cpp +++ b/src/testing/euler_gaussian_bump.cpp @@ -25,18 +25,21 @@ #include "dg/dg.h" #include "ode_solver/ode_solver.h" -//#include "template_instantiator.h" - namespace PHiLiP { namespace Tests { +/// Function used to evaluate farfield conservative solution template class FreeStreamInitialConditions : public dealii::Function { public: - std::array far_field_conservative; + /// Farfield conservative solution + std::array farfield_conservative; + /// Constructor. + /** Evaluates the primary farfield solution and converts it into the store farfield_conservative solution + */ FreeStreamInitialConditions (const Physics::Euler euler_physics) : dealii::Function(nstate) { @@ -46,14 +49,13 @@ class FreeStreamInitialConditions : public dealii::Function primitive_boundary_values[0] = density_bc; for (int d=0;d &/*point*/, const unsigned int istate) const { - return far_field_conservative[istate]; + return farfield_conservative[istate]; } }; template class FreeStreamInitialConditions ; @@ -167,7 +169,7 @@ ::run_test () const FreeStreamInitialConditions initial_conditions(euler_physics_double); pcout << "Farfield conditions: "<< std::endl; for (int s=0;s fail_conv_poly; diff --git a/src/testing/euler_gaussian_bump.h b/src/testing/euler_gaussian_bump.h index 8e5e6295b..6cb82c889 100644 --- a/src/testing/euler_gaussian_bump.h +++ b/src/testing/euler_gaussian_bump.h @@ -11,13 +11,14 @@ namespace PHiLiP { namespace Tests { +/// Gaussian bump manifold. class BumpManifold: public dealii::ChartManifold<2,2,2> { public: - virtual dealii::Point<2> pull_back(const dealii::Point<2> &space_point) const override; - virtual dealii::Point<2> push_forward(const dealii::Point<2> &chart_point) const override; - virtual dealii::DerivativeForm<1,2,2> push_forward_gradient(const dealii::Point<2> &chart_point) const override; + virtual dealii::Point<2> pull_back(const dealii::Point<2> &space_point) const override; ///< See dealii::Manifold. + virtual dealii::Point<2> push_forward(const dealii::Point<2> &chart_point) const override; ///< See dealii::Manifold. + virtual dealii::DerivativeForm<1,2,2> push_forward_gradient(const dealii::Point<2> &chart_point) const override; ///< See dealii::Manifold. - virtual std::unique_ptr > clone() const override; + virtual std::unique_ptr > clone() const override; ///< See dealii::Manifold. }; /// Performs grid convergence for various polynomial degrees. @@ -32,9 +33,7 @@ class EulerGaussianBump: public TestsBase */ EulerGaussianBump(const Parameters::AllParameters *const parameters_input); - ~EulerGaussianBump() {}; ///< Destructor. - - // Warp grid into Gaussian bump + /// Warp grid into Gaussian bump static dealii::Point warp (const dealii::Point &p); /// Grid convergence on Euler Gaussian Bump @@ -52,7 +51,8 @@ class EulerGaussianBump: public TestsBase protected: - double integrate_entropy_over_domain(DGBase &dg) const; + // // Integrate entropy over the entire domain to use as a functional. + // double integrate_entropy_over_domain(DGBase &dg) const; }; diff --git a/src/testing/euler_split_inviscid_taylor_green_vortex.cpp b/src/testing/euler_split_inviscid_taylor_green_vortex.cpp index 9c057b985..5055bfebc 100644 --- a/src/testing/euler_split_inviscid_taylor_green_vortex.cpp +++ b/src/testing/euler_split_inviscid_taylor_green_vortex.cpp @@ -7,8 +7,6 @@ template EulerTaylorGreen::EulerTaylorGreen(const Parameters::AllParameters *const parameters_input) : TestsBase::TestsBase(parameters_input) -, mpi_communicator(MPI_COMM_WORLD) -, pcout(std::cout, dealii::Utilities::MPI::this_mpi_process(mpi_communicator)==0) {} //template @@ -85,7 +83,7 @@ template int EulerTaylorGreen::run_test() const { //dealii::Triangulation grid; - dealii::parallel::distributed::Triangulation grid(this->mpi_communicator); + dealii::parallel::distributed::Triangulation grid(mpi_communicator); double left = 0.0; double right = 2 * dealii::numbers::PI; diff --git a/src/testing/euler_split_inviscid_taylor_green_vortex.h b/src/testing/euler_split_inviscid_taylor_green_vortex.h index 58a3d9314..465292574 100644 --- a/src/testing/euler_split_inviscid_taylor_green_vortex.h +++ b/src/testing/euler_split_inviscid_taylor_green_vortex.h @@ -38,20 +38,33 @@ namespace PHiLiP { namespace Tests { +/// Euler Taylor Green Vortex +/** Ensure that the kinetic energy is bounded. + * Gassner 2016. + */ template class EulerTaylorGreen : public TestsBase { public: - EulerTaylorGreen() = delete; + /// Constructor. + /** Simply calls the TestsBase constructor to set its parameters = parameters_input + */ EulerTaylorGreen(const Parameters::AllParameters *const parameters_input); + +/// Ensure that the kinetic energy is bounded. +/* If the kinetic energy increases about its initial value, then the test should fail. + * CURRENTLY PASSES NO MATTER WHAT. TO BE FIXED. + * Gassner 2016. + */ int run_test() const override; private: + /// Computes an integral of the kinetic energy (density * velocity squared) in the entire domain. + /** Overintegration of kinetic energy. + */ double compute_kinetic_energy(std::shared_ptr < DGBase > &dg, unsigned int poly_degree) const; //double compute_quadrature_kinetic_energy(std::array soln_at_q) const ; //const Parameters::AllParameters *const all_parameters; ///< Pointer to all parameters - const MPI_Comm mpi_communicator; - dealii::ConditionalOStream pcout; }; diff --git a/src/testing/euler_vortex.cpp b/src/testing/euler_vortex.cpp index 29624e2e2..60128b942 100644 --- a/src/testing/euler_vortex.cpp +++ b/src/testing/euler_vortex.cpp @@ -30,8 +30,6 @@ #include "dg/dg.h" #include "ode_solver/ode_solver.h" -//#include "template_instantiator.h" - namespace PHiLiP { namespace Tests { diff --git a/src/testing/euler_vortex.h b/src/testing/euler_vortex.h index 44d6e9f03..0b12155e7 100644 --- a/src/testing/euler_vortex.h +++ b/src/testing/euler_vortex.h @@ -29,15 +29,12 @@ class EulerVortexFunction : public dealii::Function const real vortex_strength, const real vortex_stddev_decay); - const Physics::Euler euler_physics; - const real vortex_characteristic_length; // R - const dealii::Point initial_vortex_center; // x_c, y_c - const real vortex_strength; // beta - const real vortex_stddev_decay; // sigma + const Physics::Euler euler_physics; ///< Euler physics. + const real vortex_characteristic_length; ///< R + const dealii::Point initial_vortex_center; ///< x_c, y_c + const real vortex_strength; ///< beta + const real vortex_stddev_decay; ///< sigma - /// Destructor - ~EulerVortexFunction() {}; - /// Manufactured solution exact value /** \code * u[s] = A[s]*sin(freq[s][0]*x)*sin(freq[s][1]*y)*sin(freq[s][2]*z); @@ -46,6 +43,7 @@ class EulerVortexFunction : public dealii::Function real value (const dealii::Point &point, const unsigned int istate = 0) const; private: + /// Exact solution using the current \p time provided by the dealii::Function class dealii::Point advected_location(const dealii::Point old_location) const; }; @@ -62,8 +60,6 @@ class EulerVortex: public TestsBase */ EulerVortex(const Parameters::AllParameters *const parameters_input); - ~EulerVortex() {}; ///< Destructor. - /// Manufactured grid convergence /** Will run the a grid convergence test for various p * on multiple grids to determine the order of convergence. diff --git a/src/testing/grid_study.cpp b/src/testing/grid_study.cpp index 66dc2da32..e3fd75243 100644 --- a/src/testing/grid_study.cpp +++ b/src/testing/grid_study.cpp @@ -9,6 +9,7 @@ #include #include +#include #include #include #include @@ -27,8 +28,6 @@ #include "dg/dg.h" #include "ode_solver/ode_solver.h" -//#include "template_instantiator.h" - namespace PHiLiP { namespace Tests { @@ -167,38 +166,75 @@ ::run_test () const dealii::ConvergenceTable convergence_table; - for (unsigned int igrid=0; igrid does not work for 1D - dealii::Triangulation grid( - typename dealii::Triangulation::MeshSmoothing( - dealii::Triangulation::smoothing_on_refinement | - dealii::Triangulation::smoothing_on_coarsening)); + dealii::Triangulation grid( + typename dealii::Triangulation::MeshSmoothing( + dealii::Triangulation::smoothing_on_refinement | + dealii::Triangulation::smoothing_on_coarsening)); #else - dealii::parallel::distributed::Triangulation grid( - this->mpi_communicator, - typename dealii::Triangulation::MeshSmoothing( - dealii::Triangulation::smoothing_on_refinement | - dealii::Triangulation::smoothing_on_coarsening)); + dealii::parallel::distributed::Triangulation grid( + this->mpi_communicator, + typename dealii::Triangulation::MeshSmoothing( + dealii::Triangulation::MeshSmoothing::smoothing_on_refinement)); + //dealii::Triangulation::smoothing_on_refinement | + //dealii::Triangulation::smoothing_on_coarsening)); #endif - - // Generate hypercube - if ( manu_grid_conv_param.grid_type == GridEnum::hypercube || manu_grid_conv_param.grid_type == GridEnum::sinehypercube ) { - - dealii::GridGenerator::subdivided_hyper_cube(grid, n_1d_cells[igrid]); - for (auto cell = grid.begin_active(); cell != grid.end(); ++cell) { - // Set a dummy boundary ID - cell->set_material_id(9002); - for (unsigned int face=0; face::faces_per_cell; ++face) { - if (cell->face(face)->at_boundary()) cell->face(face)->set_boundary_id (1000); - } + dealii::Vector estimated_error_per_cell; + for (unsigned int igrid=0; igridset_material_id(9002); + for (unsigned int face=0; face::faces_per_cell; ++face) { + if (cell->face(face)->at_boundary()) cell->face(face)->set_boundary_id (1000); } - // Warp grid if requested in input file - if (manu_grid_conv_param.grid_type == GridEnum::sinehypercube) dealii::GridTools::transform (&warp, grid); } + // Warp grid if requested in input file + if (manu_grid_conv_param.grid_type == GridEnum::sinehypercube) dealii::GridTools::transform (&warp, grid); + + // // Generate hypercube + // if ( igrid==0 && (manu_grid_conv_param.grid_type == GridEnum::hypercube || manu_grid_conv_param.grid_type == GridEnum::sinehypercube ) ) { + + // grid.clear(); + // dealii::GridGenerator::subdivided_hyper_cube(grid, n_1d_cells[igrid]); + // for (auto cell = grid.begin_active(); cell != grid.end(); ++cell) { + // // Set a dummy boundary ID + // cell->set_material_id(9002); + // for (unsigned int face=0; face::faces_per_cell; ++face) { + // if (cell->face(face)->at_boundary()) cell->face(face)->set_boundary_id (1000); + // } + // } + // // Warp grid if requested in input file + // if (manu_grid_conv_param.grid_type == GridEnum::sinehypercube) dealii::GridTools::transform (&warp, grid); + // } else { + // dealii::GridRefinement::refine_and_coarsen_fixed_number(grid, + // estimated_error_per_cell, + // 0.3, + // 0.03); + // grid.execute_coarsening_and_refinement(); + // } + + //for (int i=0; i<5;i++) { + // int icell = 0; + // for (auto cell = grid.begin_active(grid.n_levels()-1); cell!=grid.end(); ++cell) { + // if (!cell->is_locally_owned()) continue; + // icell++; + // if (icell < 2) { + // cell->set_refine_flag(); + // } + // //else if (icell%2 == 0) { + // // cell->set_refine_flag(); + // //} else if (icell%3 == 0) { + // // //cell->set_coarsen_flag(); + // //} + // } + // grid.execute_coarsening_and_refinement(); + //} // Distort grid by random amount if requested const double random_factor = manu_grid_conv_param.random_distortion; @@ -258,7 +294,7 @@ ::run_test () const // Overintegrate the error to make sure there is not integration error in the error estimate int overintegrate = 10; - dealii::QGauss quad_extra(dg->max_degree+1+overintegrate); + dealii::QGauss quad_extra(dg->max_degree+overintegrate); //dealii::MappingQ mappingq(dg->max_degree+1); dealii::FEValues fe_values_extra(*(dg->high_order_grid.mapping_fe_field), dg->fe_collection[poly_degree], quad_extra, dealii::update_values | dealii::update_JxW_values | dealii::update_quadrature_points); @@ -270,6 +306,7 @@ ::run_test () const // Integrate solution error and output error std::vector dofs_indices (fe_values_extra.dofs_per_cell); + estimated_error_per_cell.reinit(grid.n_active_cells()); for (auto cell = dg->dof_handler.begin_active(); cell!=dg->dof_handler.end(); ++cell) { if (!cell->is_locally_owned()) continue; @@ -277,6 +314,7 @@ ::run_test () const fe_values_extra.reinit (cell); cell->get_dof_indices (dofs_indices); + double cell_l2error = 0; for (unsigned int iquad=0; iquadmanufactured_solution_function.value(qpoint, istate); - l2error += pow(soln_at_q[istate] - uexact, 2) * fe_values_extra.JxW(iquad); + //l2error += pow(soln_at_q[istate] - uexact, 2) * fe_values_extra.JxW(iquad); + + cell_l2error += pow(soln_at_q[istate] - uexact, 2) * fe_values_extra.JxW(iquad); } } + estimated_error_per_cell[cell->active_cell_index()] = cell_l2error; + l2error += cell_l2error; } const double l2error_mpi_sum = std::sqrt(dealii::Utilities::MPI::sum(l2error, mpi_communicator)); diff --git a/src/testing/grid_study.h b/src/testing/grid_study.h index 1e6e8e3f0..68c6ba679 100644 --- a/src/testing/grid_study.h +++ b/src/testing/grid_study.h @@ -45,7 +45,14 @@ class GridStudy: public TestsBase */ static dealii::Point warp (const dealii::Point &p); + /// Initialize the solution with the exact solution + /** This is not an issue since the discretized solution will not be the exact solution. + * Therefore, the residual does not start at 0. + */ void initialize_perturbed_solution(DGBase &dg, const Physics::PhysicsBase &physics) const; + /// L2-Integral of the solution over the entire domain. + /** Used to evaluate error of a functional. + */ double integrate_solution_over_domain(DGBase &dg) const; }; diff --git a/src/testing/tests.h b/src/testing/tests.h index ded08b57f..f216fd985 100644 --- a/src/testing/tests.h +++ b/src/testing/tests.h @@ -22,6 +22,8 @@ class TestsBase /// Constructor. Deleted the default constructor since it should not be used TestsBase () = delete; /// Constructor. + /** @param[in] parameters_input Input parameters. + */ TestsBase(const Parameters::AllParameters *const parameters_input); /// Destructor. @@ -33,10 +35,16 @@ class TestsBase */ virtual int run_test() const = 0; protected: - const MPI_Comm mpi_communicator; + const MPI_Comm mpi_communicator; ///< MPI communicator. + /// ConditionalOStream. + /** Used as std::cout, but only prints if mpi_rank == 0 + */ dealii::ConditionalOStream pcout; /// Evaluates the number of cells to generate the grids for 1D grid based on input file. + /** @param[in] ngrids Number of grid sequences to generate. + * \return Vector of 1D grid sizes + */ std::vector get_number_1d_cells(const int ngrids) const; // /// Evaluates the number of cells to generate the grids for 1D grid based on input file. @@ -54,9 +62,15 @@ class TestsFactory * * TestBase test = TestFactory::create_test<3,5>(parameters_input) */ + /** @param[in] parameters_input Input parameters. + * \return Smart pointer to the test + */ static std::unique_ptr< TestsBase > create_test(const Parameters::AllParameters *const parameters_input); /// Selects the actual test such as grid convergence, numerical flux conversation, etc. + /** @param[in] parameters_input Input parameters. + * \return Smart pointer to the test + */ static std::unique_ptr< TestsBase > select_test(const Parameters::AllParameters *const parameters_input); }; diff --git a/tests/unit_tests/euler_unit_test/CMakeLists.txt b/tests/unit_tests/euler_unit_test/CMakeLists.txt index 59c89fbda..3029d4bd2 100644 --- a/tests/unit_tests/euler_unit_test/CMakeLists.txt +++ b/tests/unit_tests/euler_unit_test/CMakeLists.txt @@ -2,6 +2,7 @@ set(TEST_SRC euler_convert_primitive_conservative.cpp ) +if(NOT DOC_ONLY) foreach(dim RANGE 1 3) # Output executable @@ -99,3 +100,4 @@ foreach(dim RANGE 1 3) unset(PhysicsLib) endforeach() +endif() diff --git a/tests/unit_tests/grid/CMakeLists.txt b/tests/unit_tests/grid/CMakeLists.txt index c9b599e6b..fa2ab8331 100644 --- a/tests/unit_tests/grid/CMakeLists.txt +++ b/tests/unit_tests/grid/CMakeLists.txt @@ -21,7 +21,9 @@ foreach(dim RANGE 2 3) target_link_libraries(${TEST_TARGET} ${ParametersLib}) target_link_libraries(${TEST_TARGET} ${DiscontinuousGalerkinLib}) # Setup target with deal.II - DEAL_II_SETUP_TARGET(${TEST_TARGET}) + if (NOT DOC_ONLY) + DEAL_II_SETUP_TARGET(${TEST_TARGET}) + endif() add_test( NAME ${TEST_TARGET} @@ -34,3 +36,79 @@ foreach(dim RANGE 2 3) unset(ParametersLib) endforeach() + +set(TEST_SRC + RBF_mesh_movement.cpp + ) + +foreach(dim RANGE 2 3) + + # Output executable + string(CONCAT TEST_TARGET ${dim}D_RBF_mesh_movement) + message("Adding executable " ${TEST_TARGET} " with files " ${TEST_SRC} "\n") + add_executable(${TEST_TARGET} ${TEST_SRC}) + # Replace occurences of PHILIP_DIM with 1, 2, or 3 in the code + target_compile_definitions(${TEST_TARGET} PRIVATE PHILIP_DIM=${dim}) + + # Compile this executable when 'make unit_tests' + add_dependencies(unit_tests ${TEST_TARGET}) + add_dependencies(${dim}D ${TEST_TARGET}) + + # Library dependency + set(ParametersLib ParametersLibrary) + string(CONCAT DiscontinuousGalerkinLib DiscontinuousGalerkin_${dim}D) + target_link_libraries(${TEST_TARGET} ${ParametersLib}) + target_link_libraries(${TEST_TARGET} ${DiscontinuousGalerkinLib}) + # Setup target with deal.II + if (NOT DOC_ONLY) + DEAL_II_SETUP_TARGET(${TEST_TARGET}) + endif() + + add_test( + NAME ${TEST_TARGET} + COMMAND mpirun -n ${MPIMAX} ${EXECUTABLE_OUTPUT_PATH}/${TEST_TARGET} + WORKING_DIRECTORY ${TEST_OUTPUT_DIR} + ) + + unset(TEST_TARGET) + unset(DiscontinuousGalerkinLib) + unset(ParametersLib) + +endforeach() + +set(TEST_SRC + make_cells_valid.cpp + ) + +set(dim 2) +# Output executable +string(CONCAT TEST_TARGET ${dim}D_make_cells_valid) +message("Adding executable " ${TEST_TARGET} " with files " ${TEST_SRC} "\n") +add_executable(${TEST_TARGET} ${TEST_SRC}) +# Replace occurences of PHILIP_DIM with 1, 2, or 3 in the code +target_compile_definitions(${TEST_TARGET} PRIVATE PHILIP_DIM=${dim}) + +# Compile this executable when 'make unit_tests' +add_dependencies(unit_tests ${TEST_TARGET}) +add_dependencies(${dim}D ${TEST_TARGET}) + +# Library dependency +set(ParametersLib ParametersLibrary) +string(CONCAT DiscontinuousGalerkinLib DiscontinuousGalerkin_${dim}D) +target_link_libraries(${TEST_TARGET} ${ParametersLib}) +target_link_libraries(${TEST_TARGET} ${DiscontinuousGalerkinLib}) +# Setup target with deal.II +if (NOT DOC_ONLY) + DEAL_II_SETUP_TARGET(${TEST_TARGET}) +endif() + +add_test( + NAME ${TEST_TARGET} + COMMAND mpirun -n ${MPIMAX} ${EXECUTABLE_OUTPUT_PATH}/${TEST_TARGET} + WORKING_DIRECTORY ${TEST_OUTPUT_DIR} +) + +unset(TEST_TARGET) +unset(DiscontinuousGalerkinLib) +unset(ParametersLib) + diff --git a/tests/unit_tests/grid/RBF_mesh_movement.cpp b/tests/unit_tests/grid/RBF_mesh_movement.cpp new file mode 100644 index 000000000..1e32d2698 --- /dev/null +++ b/tests/unit_tests/grid/RBF_mesh_movement.cpp @@ -0,0 +1,175 @@ +#include +#include +#include + +#include + +#include +#include +#include +#include +#include +#include + +#include +#include + +#include +#include +#include // ExcTransformationFailed + +#include +#include + +#include "dg/high_order_grid.h" +#include "parameters/all_parameters.h" + +// template +// void output_grid(const std::string prefix, const unsigned int poly_degree, PHiLiP::HighOrderGrid &high_order_grid) +// { +// const int mpi_rank = dealii::Utilities::MPI::this_mpi_process(MPI_COMM_WORLD); +// // Output for Paraview visualization +// dealii::DataOut data_out; +// data_out.attach_dof_handler(high_order_grid.dof_handler_grid); +// std::vector solution_names; +// for (int d=0;d::CurvedCellRegion::curved_inner_cells); +// std::string filename = prefix + dealii::Utilities::int_to_string(dim, 1) +"D-"; +// filename += "Degree" + dealii::Utilities::int_to_string(poly_degree, 2) + "."; +// filename += dealii::Utilities::int_to_string(mpi_rank, 4); +// filename += ".vtu"; +// std::ofstream output(filename); +// data_out.write_vtu(output); +// +// if (mpi_rank == 0) { +// std::vector filenames; +// for (unsigned int mpi_rank = 0; mpi_rank < dealii::Utilities::MPI::n_mpi_processes(MPI_COMM_WORLD); ++mpi_rank) { +// std::string fn = prefix + dealii::Utilities::int_to_string(dim, 1) +"D-"; +// fn += "Degree" + dealii::Utilities::int_to_string(poly_degree, 2) + "."; +// fn += dealii::Utilities::int_to_string(mpi_rank, 4); +// fn += ".vtu"; +// filenames.push_back(fn); +// } +// std::string master_fn = prefix + dealii::Utilities::int_to_string(dim, 1) +"D-"; +// master_fn += "Degree" + dealii::Utilities::int_to_string(poly_degree, 2) + "."; +// master_fn += ".pvtu"; +// std::ofstream master_output(master_fn); +// data_out.write_pvtu_record(master_output, filenames); +// } +// } + +int main (int argc, char * argv[]) +{ + const int dim = PHILIP_DIM; + int fail_bool = false; + + dealii::Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, 1); + const int mpi_rank = dealii::Utilities::MPI::this_mpi_process(MPI_COMM_WORLD); + dealii::ConditionalOStream pcout(std::cout, mpi_rank==0); + + using namespace PHiLiP; + + dealii::ParameterHandler parameter_handler; + Parameters::AllParameters::declare_parameters (parameter_handler); + Parameters::AllParameters all_parameters; + all_parameters.parse_parameters (parameter_handler); + + std::vector convergence_table_vector; + + const int poly_degree = 4; + const unsigned int n_grids = 3; + //const std::vector n_1d_cells = {2,4,8,16}; + + //const unsigned int n_cells_circle = n_1d_cells[0]; + //const unsigned int n_cells_radial = 3*n_cells_circle; + + dealii::Point center; // Constructor initializes Point at the origin + const double inner_radius = 1, outer_radius = inner_radius*10; + + // Generate the original grid and assign a manifold to it + dealii::parallel::distributed::Triangulation grid( + MPI_COMM_WORLD, + typename dealii::Triangulation::MeshSmoothing( + dealii::Triangulation::smoothing_on_refinement | + dealii::Triangulation::smoothing_on_coarsening)); + + const int n_cells = 0; + const bool colorize = true; + dealii::GridGenerator::quarter_hyper_shell(grid, center, inner_radius, outer_radius, n_cells, colorize); + // Set a spherical manifold + // Works but we will usually not have a volume manifold that we can provide. + //grid.set_all_manifold_ids(0); + //grid.set_manifold(0, dealii::SphericalManifold(center)); + + // Set a spherical manifold on the boundary and a TransfiniteInterpolationManifold in the domain + // This is more realistic with what we will be doing since we will usually provide and boundary parametrization + // but have no idea of the volume parametrization. The initial curving of the TransfiniteInterpolationManifold + // will ensure that no cells initially have a negative Jacobians. This curvature will be transfered to a polynomial + // representation using MappingFEField + grid.set_all_manifold_ids(1); + grid.set_all_manifold_ids_on_boundary(0); + grid.set_manifold(0, dealii::SphericalManifold(center)); + dealii::TransfiniteInterpolationManifold transfinite_interpolation; + transfinite_interpolation.initialize(grid); + grid.set_manifold(1, transfinite_interpolation); + + + HighOrderGrid high_order_grid(&all_parameters, poly_degree, &grid); + //std::shared_ptr < dealii::MappingFEField, dealii::DoFHandler> > mapping = (high_order_grid.mapping_fe_field); + grid.reset_all_manifolds(); + + dealii::ConvergenceTable convergence_table; + std::vector grid_size(n_grids); + std::vector volume_error(n_grids); + + for (unsigned int igrid=0; igridis_locally_owned()) continue; + // icell++; + // if (icell < 3) { + // cell->set_refine_flag(); + // } else if (icell%5 == 0) { + // cell->set_refine_flag(); + // } else if (icell%3 == 0) { + // cell->set_coarsen_flag(); + // } + //} + //grid.execute_coarsening_and_refinement(); + + high_order_grid.execute_coarsening_and_refinement(); + + high_order_grid.prepare_for_coarsening_and_refinement(); + grid.repartition(); + high_order_grid.execute_coarsening_and_refinement(true); + } + + //output_grid("before", poly_degree, high_order_grid); + // Deform the y = 0 face + for (auto indices = high_order_grid.locally_owned_surface_nodes_indices.begin(); indices!=high_order_grid.locally_owned_surface_nodes_indices.end(); ++indices) { + } + + + //const unsigned int n_dofs = high_order_grid.dof_handler_grid.n_dofs(); + //const unsigned int n_global_active_cells = grid.n_global_active_cells(); + + + if (fail_bool) { + pcout << "Test failed. The estimated error should be the same for a given p, even after refinement and translation." << std::endl; + } else { + pcout << "Test successful." << std::endl; + } + return fail_bool; +} + diff --git a/tests/unit_tests/grid/high_order_grid_test.cpp b/tests/unit_tests/grid/high_order_grid_test.cpp index 53323ff8a..eefa7ef39 100644 --- a/tests/unit_tests/grid/high_order_grid_test.cpp +++ b/tests/unit_tests/grid/high_order_grid_test.cpp @@ -123,6 +123,10 @@ int main (int argc, char * argv[]) high_order_grid.execute_coarsening_and_refinement(); + high_order_grid.prepare_for_coarsening_and_refinement(); + grid.repartition(); + high_order_grid.execute_coarsening_and_refinement(true); + const unsigned int n_dofs = high_order_grid.dof_handler_grid.n_dofs(); const unsigned int n_global_active_cells = grid.n_global_active_cells(); diff --git a/tests/unit_tests/grid/make_cells_valid.cpp b/tests/unit_tests/grid/make_cells_valid.cpp new file mode 100644 index 000000000..ec23e36cb --- /dev/null +++ b/tests/unit_tests/grid/make_cells_valid.cpp @@ -0,0 +1,162 @@ +#include +#include + +#include + +#include +#include +#include +#include +#include + +#include + +#include + +#include "dg/high_order_grid.h" +#include "parameters/all_parameters.h" + +dealii::Point<2> center(0.0,0.0); +const double inner_radius = 1, outer_radius = inner_radius*20; + +dealii::Point<2> warp_cylinder (const dealii::Point<2> &p)//, const double inner_radius, const double outer_radius) +{ + const double rectangle_height = 1.0; + //const double original_radius = std::abs(p[0]); + const double angle = p[1]/rectangle_height * dealii::numbers::PI; + + //const double radius = std::abs(p[0]); + + const double power = 2.25; + const double radius = outer_radius*(inner_radius/outer_radius + pow(std::abs(p[0]), power)); + + dealii::Point<2> q = p; + q[0] = -radius*cos(angle); + q[1] = radius*sin(angle); + return q; +} + +void half_cylinder(dealii::parallel::distributed::Triangulation<2> & tria, + const unsigned int n_cells_circle, + const unsigned int n_cells_radial) +{ + dealii::Point<2> p1(-1,0.0), p2(-0.0,1.0); + + const bool colorize = true; + + std::vector n_subdivisions(2); + n_subdivisions[0] = n_cells_radial; + n_subdivisions[1] = n_cells_circle; + dealii::GridGenerator::subdivided_hyper_rectangle (tria, n_subdivisions, p1, p2, colorize); + + dealii::GridTools::transform (&warp_cylinder, tria); + + tria.set_all_manifold_ids(1); + tria.set_all_manifold_ids_on_boundary(0); + tria.set_all_manifold_ids(0); + dealii::TransfiniteInterpolationManifold<2> inner_manifold; + inner_manifold.initialize(tria); + tria.set_manifold(0, dealii::SphericalManifold<2>(center)); + tria.set_manifold(1, inner_manifold); + + + // Assign BC + for (auto cell = tria.begin_active(); cell != tria.end(); ++cell) { + //if (!cell->is_locally_owned()) continue; + for (unsigned int face=0; face::faces_per_cell; ++face) { + if (cell->face(face)->at_boundary()) { + unsigned int current_id = cell->face(face)->boundary_id(); + if (current_id == 0) { + cell->face(face)->set_boundary_id (1004); // x_left, Farfield + } else if (current_id == 1) { + cell->face(face)->set_boundary_id (1001); // x_right, Symmetry/Wall + } else if (current_id == 2) { + cell->face(face)->set_boundary_id (1001); // y_bottom, Symmetry/Wall + } else if (current_id == 3) { + cell->face(face)->set_boundary_id (1001); // y_top, Wall + } else { + std::abort(); + } + } + } + } +} + +int main (int argc, char * argv[]) +{ + const int dim = PHILIP_DIM; + + dealii::Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, 1); + const int mpi_rank = dealii::Utilities::MPI::this_mpi_process(MPI_COMM_WORLD); + dealii::ConditionalOStream pcout(std::cout, mpi_rank==0); + + const unsigned int p_start = 1; + const unsigned int p_end = 4; + + const unsigned int n_grids = 3; + + bool has_invalid_poly = false; + for (unsigned int poly_degree = p_start; poly_degree <= p_end; ++poly_degree) { + + // Generate grid and mapping + dealii::parallel::distributed::Triangulation grid(MPI_COMM_WORLD, + typename dealii::Triangulation::MeshSmoothing( + dealii::Triangulation::MeshSmoothing::smoothing_on_refinement | + dealii::Triangulation::MeshSmoothing::smoothing_on_coarsening)); + + // This grid is known to generated a negative Jacobian as seen + // in the EulerCylinder test case described in + // https://github.com/dougshidong/PHiLiP/issues/27 + const unsigned int n_cells_circle = 2; + const unsigned int n_cells_radial = 2*n_cells_circle; + half_cylinder(grid, n_cells_circle, n_cells_radial); + + dealii::ParameterHandler parameter_handler; + PHiLiP::Parameters::AllParameters::declare_parameters (parameter_handler); + PHiLiP::Parameters::AllParameters all_parameters; + all_parameters.parse_parameters (parameter_handler); + + PHiLiP::HighOrderGrid high_order_grid(&all_parameters, poly_degree, &grid); + + //bool has_invalid_grid = false; + for (unsigned int igrid=0; igrid0) { + high_order_grid.prepare_for_coarsening_and_refinement(); + grid.refine_global (1); + high_order_grid.execute_coarsening_and_refinement(); + + high_order_grid.prepare_for_coarsening_and_refinement(); + grid.repartition(); + high_order_grid.execute_coarsening_and_refinement(true); + } + + auto cell = high_order_grid.dof_handler_grid.begin_active(); + auto endcell = high_order_grid.dof_handler_grid.end(); + + for (; cell!=endcell; ++cell) { + if (!cell->is_locally_owned()) continue; + const bool is_invalid_cell = high_order_grid.check_valid_cell(cell); + + if ( !is_invalid_cell ) { + std::cout << " Poly: " << poly_degree + << " Grid: " << igrid + << " Cell: " << cell->active_cell_index() << " has an invalid Jacobian." << std::endl; + //has_invalid_grid = true; + bool fixed_invalid_cell = high_order_grid.fix_invalid_cell(cell); + if (fixed_invalid_cell) std::cout << "Fixed it." << std::endl; + else has_invalid_poly = true; + } + } + high_order_grid.nodes.update_ghost_values(); + if (has_invalid_poly) std::abort(); + } + } + + bool mpi_has_invalid_poly; + MPI_Allreduce(&has_invalid_poly, &mpi_has_invalid_poly, 1, MPI::BOOL, MPI_LOR, MPI_COMM_WORLD); + + return mpi_has_invalid_poly; +} + diff --git a/tests/unit_tests/numerical_flux/CMakeLists.txt b/tests/unit_tests/numerical_flux/CMakeLists.txt index b5b4e26b9..ffff4077f 100644 --- a/tests/unit_tests/numerical_flux/CMakeLists.txt +++ b/tests/unit_tests/numerical_flux/CMakeLists.txt @@ -24,7 +24,9 @@ foreach(dim RANGE 1 3) target_link_libraries(${TEST_TARGET} ${PhysicsLib}) target_link_libraries(${TEST_TARGET} ${NumericalFluxLib}) # Setup target with deal.II - DEAL_II_SETUP_TARGET(${TEST_TARGET}) + if(NOT DOC_ONLY) + DEAL_II_SETUP_TARGET(${TEST_TARGET}) + endif() add_test( NAME ${TEST_TARGET} diff --git a/tests/unit_tests/regression/CMakeLists.txt b/tests/unit_tests/regression/CMakeLists.txt index ed000e473..3cf9ea4fc 100644 --- a/tests/unit_tests/regression/CMakeLists.txt +++ b/tests/unit_tests/regression/CMakeLists.txt @@ -23,7 +23,9 @@ foreach(dim RANGE 1 3) target_link_libraries(${TEST_TARGET} ${ParametersLib}) target_link_libraries(${TEST_TARGET} ${DiscontinuousGalerkinLib}) # Setup target with deal.II - DEAL_II_SETUP_TARGET(${TEST_TARGET}) + if(NOT DOC_ONLY) + DEAL_II_SETUP_TARGET(${TEST_TARGET}) + endif() add_test( NAME ${TEST_TARGET}