diff --git a/.clang-tidy b/.clang-tidy index 6919bf1b..03b6ebb1 100644 --- a/.clang-tidy +++ b/.clang-tidy @@ -9,6 +9,7 @@ Checks: > cppcoreguidelines-*, -cppcoreguidelines-avoid-const-or-ref-data-members, -cppcoreguidelines-non-private-member-variables-in-classes, + -cppcoreguidelines-avoid-magic-numbers, hicpp-*, misc-*, -misc-non-private-member-variables-in-classes, @@ -19,5 +20,6 @@ Checks: > readability-*, -readability-else-after-return, -readability-identifier-length, + -readability-magic-numbers, WarningsAsErrors: '*' diff --git a/.github/workflows/clang-tidy.yml b/.github/workflows/clang-tidy.yml index 243e51b9..38500019 100644 --- a/.github/workflows/clang-tidy.yml +++ b/.github/workflows/clang-tidy.yml @@ -7,7 +7,6 @@ on: pull_request: branches: - master - - devel types: - opened - reopened diff --git a/.github/workflows/linux.yml b/.github/workflows/linux.yml index 48b3b329..cf7bf9d7 100644 --- a/.github/workflows/linux.yml +++ b/.github/workflows/linux.yml @@ -7,7 +7,6 @@ on: pull_request: branches: - master - - devel types: - opened - reopened diff --git a/.gitignore b/.gitignore index 33d26c1c..231a8aca 100644 --- a/.gitignore +++ b/.gitignore @@ -59,6 +59,7 @@ main #Output files applications/*/*.vtk *.vtu +*.pvtu *.vts *.frt integratedFields.txt diff --git a/applications/main.cc b/applications/main.cc index f5526507..c9fd00bd 100644 --- a/applications/main.cc +++ b/applications/main.cc @@ -5,8 +5,8 @@ #include "core/variableAttributes.h" #include "equations.cc" -#include #include +#include #include // Header file for postprocessing that may or may not exist @@ -38,7 +38,7 @@ main(int argc, char **argv) std::string parameters_filename; try { - ParseCommandLineOpts cli_options(argc, argv); + parseCMDOptions cli_options(argc, argv); parameters_filename = cli_options.getParametersFilename(); } catch (const char *msg) diff --git a/applications/multigrid_test/CMakeLists.txt b/applications/multigrid_test/CMakeLists.txt new file mode 100644 index 00000000..7780059a --- /dev/null +++ b/applications/multigrid_test/CMakeLists.txt @@ -0,0 +1,157 @@ +## +# CMake script for the PRISMS-PF applications +# Adapted from the ASPECT CMake file +## + +cmake_minimum_required(VERSION 3.3.0) + +project(myapp) + +message(STATUS "") +message(STATUS "=========================================================") +message(STATUS "Configuring PRISMS-PF application") +message(STATUS "=========================================================") +message(STATUS "") + +# Setup the build type (debug, release, debugrelease) +if("${CMAKE_BUILD_TYPE}" STREQUAL "") + set(CMAKE_BUILD_TYPE "DebugRelease" + CACHE STRING + "Choose the type of build, options are: Debug, Release and DebugRelease." + FORCE) +endif() + +# Convert build type into the debug and release builds, which may or may +# not be built. +if("${CMAKE_BUILD_TYPE}" STREQUAL "Release" OR + "${CMAKE_BUILD_TYPE}" STREQUAL "Debug" OR + "${CMAKE_BUILD_TYPE}" STREQUAL "DebugRelease" ) + message(STATUS "Setting up PRISMS-PF application for ${CMAKE_BUILD_TYPE} mode.") +else() + message(FATAL_ERROR + "CMAKE_BUILD_TYPE must either be 'Release', 'Debug', or 'DebugRelease', but is set to '${CMAKE_BUILD_TYPE}'.") +endif() + +if("${CMAKE_BUILD_TYPE}" STREQUAL "Debug" OR + "${CMAKE_BUILD_TYPE}" STREQUAL "DebugRelease") + set(PRISMS_PF_BUILD_DEBUG "ON") +else() + set(PRISMS_PF_BUILD_DEBUG "OFF") +endif() + +if("${CMAKE_BUILD_TYPE}" STREQUAL "Release" OR + "${CMAKE_BUILD_TYPE}" STREQUAL "DebugRelease") + set(PRISMS_PF_BUILD_RELEASE "ON") +else() + set(PRISMS_PF_BUILD_RELEASE "OFF") +endif() + +# Find deal.II installation +find_package(deal.II 9.6.0 QUIET + HINTS ${DEAL_II_DIR} $ENV{DEAL_II_DIR} + ) +if(NOT ${deal.II_FOUND}) + message(FATAL_ERROR "\n*** Could not find a recent version of deal.II. ***\n" + "You may want to either pass a flag -DDEAL_II_DIR=/path/to/deal.II to cmake " + "or set an environment variable \"DEAL_II_DIR\" that contains a path to a " + "recent version of deal.II." + ) +endif() + +message(STATUS "Found deal.II version ${DEAL_II_PACKAGE_VERSION} at '${deal.II_DIR}'") + +set(DEALII_INSTALL_VALID ON) + +if(NOT DEAL_II_WITH_P4EST) + message(SEND_ERROR + "\n**deal.II was built without support for p4est!\n" + ) + set(DEALII_INSTALL_VALID OFF) +endif() + +if(NOT DEALII_INSTALL_VALID) + message(FATAL_ERROR + "\nPRISMS-PF requires a deal.II installation with certain features enabled!\n" + ) +endif() + +deal_ii_initialize_cached_variables() + +# Make and ninja build options +if(CMAKE_GENERATOR MATCHES "Ninja") + set(_make_command "$ ninja") +else() + set(_make_command " $ make") +endif() + +# Debug and release targets +if(${DEAL_II_BUILD_TYPE} MATCHES "DebugRelease") + add_custom_target(release + COMMAND ${CMAKE_COMMAND} -D CMAKE_BUILD_TYPE=Release . + COMMAND ${CMAKE_COMMAND} -E echo "***" + COMMAND ${CMAKE_COMMAND} -E echo "*** Switched to Release mode. Now recompile with: ${_make_command}" + COMMAND ${CMAKE_COMMAND} -E echo "***" + WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR} + VERBATIM + COMMENT "switching to RELEASE mode..." + ) + + add_custom_target(debug + COMMAND ${CMAKE_COMMAND} -D CMAKE_BUILD_TYPE=Debug . + COMMAND ${CMAKE_COMMAND} -E echo "***" + COMMAND ${CMAKE_COMMAND} -E echo "*** Switched to Debug mode. Now recompile with: ${_make_command}" + COMMAND ${CMAKE_COMMAND} -E echo "***" + WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR} + VERBATIM + COMMENT "switching to DEBUG mode..." + ) + + add_custom_target(debugrelease + COMMAND ${CMAKE_COMMAND} -D CMAKE_BUILD_TYPE=DebugRelease . + COMMAND ${CMAKE_COMMAND} -E echo "***" + COMMAND ${CMAKE_COMMAND} -E echo "*** Switched to Debug and Release mode. Now recompile with: ${_make_command}" + COMMAND ${CMAKE_COMMAND} -E echo "***" + WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR} + VERBATIM + COMMENT "switching to DEBUG/RELEASE mode..." + ) +endif() + +# Add distclean target to clean build +add_custom_target(distclean + COMMAND ${CMAKE_COMMAND} --build ${CMAKE_BINARY_DIR} --target clean + COMMAND ${CMAKE_COMMAND} -E remove_directory CMakeFiles + COMMAND ${CMAKE_COMMAND} -E remove + CMakeCache.txt cmake_install.cmake Makefile + build.ninja rules.ninja .ninja_deps .ninja_log + COMMENT "distclean invoked" + ) + +# Set location of files +include_directories(${CMAKE_SOURCE_DIR}/../../include) +include_directories(${CMAKE_SOURCE_DIR}/../../src) +include_directories(${CMAKE_SOURCE_DIR}) + +# Set the location of the main.cc file +set(TARGET_SRC "${CMAKE_SOURCE_DIR}/main.cc" "${CMAKE_SOURCE_DIR}/equations.cc") + +# Check if there has been updates to main library +set(dir ${PROJECT_SOURCE_DIR}/../..) +EXECUTE_PROCESS(COMMAND "rm" "CMakeCache.txt" WORKING_DIRECTORY ${dir}) +EXECUTE_PROCESS(COMMAND "cmake" "CMakeLists.txt" WORKING_DIRECTORY ${dir}) +EXECUTE_PROCESS(COMMAND "make" WORKING_DIRECTORY ${dir}) + +# Set targets & link libraries for the build type +if(${PRISMS_PF_BUILD_DEBUG} STREQUAL "ON") + add_executable(main_debug ${TARGET_SRC}) + set_property(TARGET main_debug PROPERTY OUTPUT_NAME main-debug) + deal_ii_setup_target(main_debug DEBUG) + target_link_libraries(main_debug ${CMAKE_SOURCE_DIR}/../../libprisms-pf-debug.a) +endif() + +if(${PRISMS_PF_BUILD_RELEASE} STREQUAL "ON") + add_executable(main_release ${TARGET_SRC}) + set_property(TARGET main_release PROPERTY OUTPUT_NAME main) + deal_ii_setup_target(main_release RELEASE) + target_link_libraries(main_release ${CMAKE_SOURCE_DIR}/../../libprisms-pf-release.a) +endif() \ No newline at end of file diff --git a/applications/multigrid_test/custom_pde.h b/applications/multigrid_test/custom_pde.h new file mode 100644 index 00000000..ee568561 --- /dev/null +++ b/applications/multigrid_test/custom_pde.h @@ -0,0 +1,57 @@ +#ifndef CUSTOM_PDE_H_ +#define CUSTOM_PDE_H_ + +#include + +using namespace dealii; + +/** + * \brief This is a derived class of `matrixFreeOperator` where the user implements their + * PDEs. + * + * \tparam dim The number of dimensions in the problem. + * \tparam degree The polynomial degree of the shape functions. + * \tparam number Datatype to use. Either double or float. + */ +template +class customPDE : public matrixFreeOperator +{ +public: + using scalarValue = VectorizedArray; + using scalarGrad = Tensor<1, dim, VectorizedArray>; + using scalarHess = Tensor<2, dim, VectorizedArray>; + using vectorValue = Tensor<1, dim, VectorizedArray>; + using vectorGrad = Tensor<2, dim, VectorizedArray>; + using vectorHess = Tensor<3, dim, VectorizedArray>; + + /** + * \brief Constructor. + */ + customPDE() = default; + + /** + * \brief User-implemented class for the RHS of explicit equations. + */ + void + compute_explicit_RHS(variableContainer &variable_list, + const dealii::Point> + &q_point_loc) const override; + + /** + * \brief User-implemented class for the RHS of nonexplicit equations. + */ + void + compute_nonexplicit_RHS(variableContainer &variable_list, + const dealii::Point> + &q_point_loc) const override; + + /** + * \brief User-implemented class for the LHS of nonexplicit equations. + */ + void + compute_nonexplicit_LHS(variableContainer &variable_list, + const dealii::Point> + &q_point_loc) const override; +}; + +#endif \ No newline at end of file diff --git a/applications/multigrid_test/equations.cc b/applications/multigrid_test/equations.cc new file mode 100644 index 00000000..ef428c76 --- /dev/null +++ b/applications/multigrid_test/equations.cc @@ -0,0 +1,48 @@ +#include "custom_pde.h" + +#include +#include + +void +customAttributeLoader::loadVariableAttributes() +{ + set_variable_name(0, "phi"); + set_variable_type(0, SCALAR); + set_variable_equation_type(0, EXPLICIT_TIME_DEPENDENT); + + set_dependencies_value_term_RHS(0, ""); + set_dependencies_gradient_term_RHS(0, ""); +} + +void +customAttributeLoader::loadPostProcessorVariableAttributes() +{} + +template +void +customPDE::compute_explicit_RHS( + variableContainer &variable_list, + const Point> &q_point_loc) const +{} + +template +void +customPDE::compute_nonexplicit_RHS( + variableContainer &variable_list, + const Point> &q_point_loc) const +{ + scalarGrad phi = variable_list.get_scalar_gradient(0); + + scalarValue coefficient = 1.0 / (0.05 + 2.0 * q_point_loc.square()); + + variable_list.set_scalar_gradient_term(0, phi * coefficient); +} + +template +void +customPDE::compute_nonexplicit_LHS( + variableContainer &variable_list, + const Point> &q_point_loc) const +{} + +INSTANTIATE_TRI_TEMPLATE(customPDE) \ No newline at end of file diff --git a/applications/multigrid_test/main.cc b/applications/multigrid_test/main.cc new file mode 100644 index 00000000..f3c561f0 --- /dev/null +++ b/applications/multigrid_test/main.cc @@ -0,0 +1,38 @@ +#include "custom_pde.h" + +#include + +int +main(int argc, char *argv[]) +{ + try + { + Utilities::MPI::MPI_InitFinalize mpi_init(argc, argv, 1); + + LaplaceProblem laplace_problem; + laplace_problem.run(); + } + catch (std::exception &exc) + { + std::cerr << '\n' + << '\n' + << "----------------------------------------------------" << '\n'; + std::cerr << "Exception on processing: " << '\n' + << exc.what() << '\n' + << "Aborting!" << '\n' + << "----------------------------------------------------" << '\n'; + return 1; + } + catch (...) + { + std::cerr << '\n' + << '\n' + << "----------------------------------------------------" << '\n'; + std::cerr << "Unknown exception!" << '\n' + << "Aborting!" << '\n' + << "----------------------------------------------------" << '\n'; + return 1; + } + + return 0; +} \ No newline at end of file diff --git a/automatic_tests/main.cc b/automatic_tests/main.cc index f5526507..c9fd00bd 100644 --- a/automatic_tests/main.cc +++ b/automatic_tests/main.cc @@ -5,8 +5,8 @@ #include "core/variableAttributes.h" #include "equations.cc" -#include #include +#include #include // Header file for postprocessing that may or may not exist @@ -38,7 +38,7 @@ main(int argc, char **argv) std::string parameters_filename; try { - ParseCommandLineOpts cli_options(argc, argv); + parseCMDOptions cli_options(argc, argv); parameters_filename = cli_options.getParametersFilename(); } catch (const char *msg) diff --git a/include/config.h.in b/include/config.h.in index 622a2f0f..bf643496 100644 --- a/include/config.h.in +++ b/include/config.h.in @@ -10,4 +10,46 @@ #cmakedefine PRISMS_PF_WITH_CALIPER +// Macro for template instantations with +#define INSTANTIATE_TRI_TEMPLATE(class_name) \ + template class class_name<2, 1, float>; \ + template class class_name<2, 2, float>; \ + template class class_name<2, 3, float>; \ + template class class_name<2, 4, float>; \ + template class class_name<2, 5, float>; \ + template class class_name<2, 6, float>; \ + template class class_name<2, 1, double>; \ + template class class_name<2, 2, double>; \ + template class class_name<2, 3, double>; \ + template class class_name<2, 4, double>; \ + template class class_name<2, 5, double>; \ + template class class_name<2, 6, double>; \ + template class class_name<3, 1, float>; \ + template class class_name<3, 2, float>; \ + template class class_name<3, 3, float>; \ + template class class_name<3, 4, float>; \ + template class class_name<3, 5, float>; \ + template class class_name<3, 6, float>; \ + template class class_name<3, 1, double>; \ + template class class_name<3, 2, double>; \ + template class class_name<3, 3, double>; \ + template class class_name<3, 4, double>; \ + template class class_name<3, 5, double>; \ + template class class_name<3, 6, double>; + +// Macro for template instantations with +#define INSTANTIATE_BI_TEMPLATE(class_name) \ + template class class_name<2, 1>; \ + template class class_name<2, 2>; \ + template class class_name<2, 3>; \ + template class class_name<2, 4>; \ + template class class_name<2, 5>; \ + template class class_name<2, 6>; \ + template class class_name<3, 1>; \ + template class class_name<3, 2>; \ + template class class_name<3, 3>; \ + template class class_name<3, 4>; \ + template class class_name<3, 5>; \ + template class class_name<3, 6>; + #endif diff --git a/include/core/conditional_ostreams.h b/include/core/conditional_ostreams.h new file mode 100644 index 00000000..ed27efbc --- /dev/null +++ b/include/core/conditional_ostreams.h @@ -0,0 +1,47 @@ +#ifndef CONDITIONAL_OSTREAMS_H_ +#define CONDITIONAL_OSTREAMS_H_ + +#include +#include + +#include + +/** + * \brief A class that allows printing to different output streams that are classified + * based on their verbosity. For now, this consists of two stream the release and debug. + * The debug stream provides more information that may be useful when debugging. + */ +class conditionalOStreams +{ +public: + /** + * \brief Constructor. + */ + conditionalOStreams(); + + /** + * \brief Generic parallel output stream. Used for essential information in release and + * debug mode. + */ + static const dealii::ConditionalOStream pout_base; + + /** + * \brief Verbose parallel output stream. Used for additional information in debug mode. + */ + static const dealii::ConditionalOStream pout_verbose; +}; + +// NOLINTBEGIN +const dealii::ConditionalOStream conditionalOStreams::pout_base( + std::cout, + dealii::Utilities::MPI::this_mpi_process(MPI_COMM_WORLD) == 0); + +const dealii::ConditionalOStream conditionalOStreams::pout_verbose( + std::cout, +#ifndef DEBUG + false && +#endif + dealii::Utilities::MPI::this_mpi_process(MPI_COMM_WORLD) == 0); +// NOLINTEND + +#endif \ No newline at end of file diff --git a/include/core/fields.h b/include/core/fields.h index 9a2cf8a0..b2335b0a 100644 --- a/include/core/fields.h +++ b/include/core/fields.h @@ -3,7 +3,7 @@ #include -#include +#include #include #include #include diff --git a/include/core/matrixFreePDE.h b/include/core/matrixFreePDE.h index bdc070cb..617087c8 100644 --- a/include/core/matrixFreePDE.h +++ b/include/core/matrixFreePDE.h @@ -4,7 +4,6 @@ // dealii headers #include #include -#include #include #include #include @@ -20,19 +19,20 @@ #include #include #include +#include #include #include #include #include // PRISMS headers -#include #include #include #include #include #include #include +#include #include #include @@ -78,18 +78,6 @@ class MatrixFreePDE : public Subscriptor virtual void create_triangulation(parallel::distributed::Triangulation &tria) const; - /** - * \brief Initializes the data structures for enabling unit tests. - * - * \details This method initializes the MatrixFreePDE object with a fixed geometry, - * discretization and other custom selected options specifically to help with unit - * tests, and should not be called in any of the physical models. - * - * \param _fields Vector of PDE descriptions (e.g., scalar/vector) for each field. - */ - void - initForTests(std::vector> _fields); - /** * \brief This method implements the time stepping algorithm and invokes the * solveIncrement() method. @@ -97,37 +85,6 @@ class MatrixFreePDE : public Subscriptor void solve(); - /** - * \brief This method essentially converts the MatrixFreePDE object into a matrix object - * which can be used with matrix free iterative solvers. Provides the A*x functionality - * for solving the system of equations Ax=b. - * - * \param dst The destination vector. - * \param src The source vector. - */ - void - vmult(dealii::LinearAlgebra::distributed::Vector &dst, - const dealii::LinearAlgebra::distributed::Vector &src) const; - - /** - * \brief Vector of all the physical fields in the problem. Fields are identified by - * dimentionality (SCALAR/VECTOR), the kind of PDE (ELLIPTIC/PARABOLIC) used to compute - * them and a character identifier (e.g.: "c" for composition) which is used to write - * the fields to the output files. - */ - std::vector> fields; - - /** - * \brief Create the vector of all physical fields in the problem. - */ - void - buildFields(); - - /** - * \brief Parallel message stream. - */ - ConditionalOStream pcout; - /** * \brief Set the initial condition for all fields. This function is overriden in each * application. @@ -189,6 +146,10 @@ class MatrixFreePDE : public Subscriptor std::vector> simplified_grain_representations; + // =========================================================================== + // METHODS FOR SOLVING + // =========================================================================== + /** * Method to solve each time increment of a time-dependent problem. For * time-independent problems this method is called only once. This method @@ -197,7 +158,51 @@ class MatrixFreePDE : public Subscriptor * problems, Implicit (matrix-free) solver for Elliptic problems. */ virtual void - solveIncrement(bool skip_time_dependent); + solveIncrement(uint current_increment, bool skip_time_dependent); + + void + linearSolveOnce(unsigned int var_index, SolverControl &solver_control); + void + auxiliaryOnce(unsigned int var_index, uint current_increment); + bool + nonlinearIncrement(unsigned int var_index, + SolverControl &solver_control, + uint current_increment); + bool + auxiliaryIncrement(unsigned int var_index, uint current_increment); + void + print_explicit_update(uint var_index, uint current_increment); + void + print_linear_update(uint var_index, + uint current_increment, + SolverControl &solver_control); + void + print_nonlinear_update(uint var_index, uint current_increment); + void + print_nonlinear_status(uint var_index, + uint current_increment, + uint nonlinear_iteration_index, + double diff); + + /*Method to compute an explicit timestep*/ + void + updateExplicitSolution(unsigned int var_index); + + void + updateLinearSolution(unsigned int var_index); + + /*Method to compute an implicit timestep*/ + bool + updateImplicitSolution(unsigned int fieldIndex, unsigned int nonlinear_iteration_index); + + /*Method to apply boundary conditions*/ + void + applyBCs(unsigned int var_index); + + // =========================================================================== + // + // =========================================================================== + /* Method to write solution fields to vtu and pvtu (parallel) files. * * This method can be enabled/disabled by setting the flag writeOutput to @@ -281,18 +286,6 @@ class MatrixFreePDE : public Subscriptor void computeInvM(); - /*Method to compute an explicit timestep*/ - void - updateExplicitSolution(unsigned int fieldIndex); - - /*Method to compute an implicit timestep*/ - bool - updateImplicitSolution(unsigned int fieldIndex, unsigned int nonlinear_iteration_index); - - /*Method to apply boundary conditions*/ - void - applyBCs(unsigned int fieldIndex); - /** * \brief Compute element volume for the triangulation */ @@ -304,93 +297,6 @@ class MatrixFreePDE : public Subscriptor */ dealii::AlignedVector> element_volume; - /*Method to compute the right hand side (RHS) residual vectors*/ - void - computeExplicitRHS(); - void - computeNonexplicitRHS(); - - // virtual methods to be implemented in the derived class - /*Method to calculate LHS(implicit solve)*/ - void - getLHS(const MatrixFree &data, - dealii::LinearAlgebra::distributed::Vector &dst, - const dealii::LinearAlgebra::distributed::Vector &src, - const std::pair &cell_range) const; - - bool generatingInitialGuess; - void - getLaplaceLHS(const MatrixFree &data, - dealii::LinearAlgebra::distributed::Vector &dst, - const dealii::LinearAlgebra::distributed::Vector &src, - const std::pair &cell_range) const; - - void - setNonlinearEqInitialGuess(); - void - computeLaplaceRHS(unsigned int fieldIndex); - void - getLaplaceRHS(const MatrixFree &data, - dealii::LinearAlgebra::distributed::Vector &dst, - const dealii::LinearAlgebra::distributed::Vector &src, - const std::pair &cell_range) const; - - /*Method to calculate RHS (implicit/explicit). This is an abstract method, so - * every model which inherits MatrixFreePDE has to implement this - * method.*/ - void - getExplicitRHS( - const MatrixFree &data, - std::vector *> &dst, - const std::vector *> &src, - const std::pair &cell_range) const; - - void - getNonexplicitRHS( - const MatrixFree &data, - std::vector *> &dst, - const std::vector *> &src, - const std::pair &cell_range) const; - - virtual void - explicitEquationRHS( - [[maybe_unused]] variableContainer> - &variable_list, - [[maybe_unused]] const Point> q_point_loc, - [[maybe_unused]] const VectorizedArray element_volume) const = 0; - - virtual void - nonExplicitEquationRHS( - [[maybe_unused]] variableContainer> - &variable_list, - [[maybe_unused]] const Point> q_point_loc, - [[maybe_unused]] const VectorizedArray element_volume) const = 0; - - virtual void - equationLHS([[maybe_unused]] variableContainer> - &variable_list, - [[maybe_unused]] const Point> q_point_loc, - [[maybe_unused]] const VectorizedArray element_volume) const = 0; - - virtual void - postProcessedFields( - [[maybe_unused]] const variableContainer> - &variable_list, - [[maybe_unused]] variableContainer> - &pp_variable_list, - [[maybe_unused]] const Point> q_point_loc, - [[maybe_unused]] const VectorizedArray element_volume) const {}; - void - computePostProcessedFields( - std::vector *> &postProcessedSet); - - void - getPostProcessedFields( - const MatrixFree &data, - std::vector *> &dst, - const std::vector *> &src, - const std::pair &cell_range); - // methods to apply dirichlet BC's /*Map of degrees of freedom to the corresponding Dirichlet boundary * conditions, if any.*/ @@ -495,11 +401,6 @@ class MatrixFreePDE : public Subscriptor return 0.0; }; - // utility functions - /*Returns index of given field name if exists, else throw error.*/ - unsigned int - getFieldIndex(std::string _name); - /*Method to compute the integral of a field.*/ void computeIntegral( @@ -520,9 +421,6 @@ class MatrixFreePDE : public Subscriptor unsigned int currentIncrement, currentOutput, currentCheckpoint, current_grain_reassignment; - /*Timer and logging object*/ - mutable TimerOutput computing_timer; - bool first_integrated_var_output_complete; // Methods and variables for integration diff --git a/include/core/matrix_free_operator.h b/include/core/matrix_free_operator.h new file mode 100644 index 00000000..ad605062 --- /dev/null +++ b/include/core/matrix_free_operator.h @@ -0,0 +1,192 @@ +#ifndef MATRIX_FREE_OPERATOR_H_ +#define MATRIX_FREE_OPERATOR_H_ + +#include +#include +#include + +#include + +/** + * \brief This is the abstract base class for the matrix free implementation of some PDE. + * + * \tparam dim The number of dimensions in the problem. + * \tparam degree The polynomial degree of the shape functions. + * \tparam number Datatype to use for `LinearAlgebra::distributed::Vector`. Either + * double or float. + */ +template +class matrixFreeOperator + : public dealii::MatrixFreeOperators:: + Base> +{ +public: + /** + * \brief Constructor. + */ + matrixFreeOperator(); + + /** + * \brief Release all memory and return to state like having called the default + * constructor. + */ + void + clear() override; + + /** + * \brief Compute the diagonal of this operator. + */ + void + compute_diagonal() override; + +protected: + /** + * \brief User-implemented class for the RHS of explicit equations. + */ + virtual void + compute_explicit_RHS( + variableContainer &variable_list, + const dealii::Point> &q_point_loc) const = 0; + + /** + * \brief User-implemented class for the RHS of nonexplicit equations. + */ + virtual void + compute_nonexplicit_RHS( + variableContainer &variable_list, + const dealii::Point> &q_point_loc) const = 0; + + /** + * \brief User-implemented class for the LHS of nonexplicit equations. + */ + virtual void + compute_nonexplicit_LHS( + variableContainer &variable_list, + const dealii::Point> &q_point_loc) const = 0; + +private: + /** + * \brief Apply operator to src and add result in dst. + */ + void + apply_add(dealii::LinearAlgebra::distributed::Vector &dst, + const dealii::LinearAlgebra::distributed::Vector &src) const override; + + /** + * \brief Local application of the operator. + */ + void + local_apply(const dealii::MatrixFree &data, + dealii::LinearAlgebra::distributed::Vector &dst, + const dealii::LinearAlgebra::distributed::Vector &src, + const std::pair &cell_range) const; + + /** + * \brief Local computation of the diagonal of the operator. + */ + void + local_compute_diagonal(const dealii::MatrixFree &data, + dealii::LinearAlgebra::distributed::Vector &dst, + const unsigned int &dummy, + const std::pair &cell_range) const; +}; + +template +matrixFreeOperator::matrixFreeOperator() + : dealii::MatrixFreeOperators:: + Base>() +{} + +template +void +matrixFreeOperator::clear() +{ + dealii::MatrixFreeOperators:: + Base>::clear(); +} + +template +void +matrixFreeOperator::local_apply( + const dealii::MatrixFree &data, + dealii::LinearAlgebra::distributed::Vector &dst, + const dealii::LinearAlgebra::distributed::Vector &src, + const std::pair &cell_range) const +{ + // Constructor for FEEvaluation objects + variableContainer variable_list(data); + + // Initialize, evaluate, and submit based on user function. + variable_list.eval_local_operator( + [this](variableContainer &var_list, + const dealii::Point> &q_point_loc) + { + this->compute_nonexplicit_RHS(var_list, q_point_loc); + }, + dst, + src, + cell_range); +} + +template +void +matrixFreeOperator::apply_add( + dealii::LinearAlgebra::distributed::Vector &dst, + const dealii::LinearAlgebra::distributed::Vector &src) const +{ + this->data->cell_loop(&matrixFreeOperator::local_apply, this, dst, src); +} + +template +void +matrixFreeOperator::compute_diagonal() +{ + this->inverse_diagonal_entries.reset( + new DiagonalMatrix>()); + dealii::LinearAlgebra::distributed::Vector &inverse_diagonal = + this->inverse_diagonal_entries->get_vector(); + this->data->initialize_dof_vector(inverse_diagonal); + unsigned int dummy = 0; + this->data->cell_loop(&matrixFreeOperator::local_compute_diagonal, + this, + inverse_diagonal, + dummy); + + this->set_constrained_entries_to_one(inverse_diagonal); + + for (unsigned int i = 0; i < inverse_diagonal.locally_owned_size(); ++i) + { + Assert(inverse_diagonal.local_element(i) > 0.0, + dealii::ExcMessage( + "No diagonal entry in a positive definite operator should be zero")); + inverse_diagonal.local_element(i) = 1.0 / inverse_diagonal.local_element(i); + } +} + +template +void +matrixFreeOperator::local_compute_diagonal( + const dealii::MatrixFree &data, + dealii::LinearAlgebra::distributed::Vector &dst, + [[maybe_unused]] const unsigned int &dummy, + const std::pair &cell_range) const +{ + // Field index + unsigned int field_index = 0; + + // Constructor for FEEvaluation objects + variableContainer variable_list(data); + + // Initialize, evaluate, and submit diagonal based on user function. + variable_list.eval_local_diagonal( + [this](variableContainer &var_list, + const dealii::Point> &q_point_loc) + { + this->compute_nonexplicit_RHS(var_list, q_point_loc); + }, + dst, + cell_range, + field_index); +} + +#endif \ No newline at end of file diff --git a/include/core/model_variables.h b/include/core/model_variables.h index 24aed755..983ca225 100644 --- a/include/core/model_variables.h +++ b/include/core/model_variables.h @@ -1,22 +1,33 @@ -// Model Variables Class - #ifndef INCLUDE_MODELVARIABLE_H_ #define INCLUDE_MODELVARIABLE_H_ -#include -#include #include -struct variable_info +#include + +struct variableInfo { + /** + * \brief Global field index. + */ unsigned int global_var_index = 0; - bool is_scalar = true; - bool var_needed = false; + /** + * \brief Variable field type (SCALAR/VECTOR). + */ + fieldType var_type = UNDEFINED_FIELD; + + /** + * \brief Evaluation flags. + */ dealii::EvaluationFlags::EvaluationFlags evaluation_flags = dealii::EvaluationFlags::nothing; + + /** + * \brief Residual evaluation flags. + */ dealii::EvaluationFlags::EvaluationFlags residual_flags = dealii::EvaluationFlags::nothing; }; -#endif /* INCLUDE_MODELVARIABLE_H_ */ +#endif diff --git a/include/core/ParseCommandLineOpts.h b/include/core/parse_cmd_options.h similarity index 61% rename from include/core/ParseCommandLineOpts.h rename to include/core/parse_cmd_options.h index bd17defe..2cc9bfa4 100644 --- a/include/core/ParseCommandLineOpts.h +++ b/include/core/parse_cmd_options.h @@ -1,22 +1,25 @@ -#ifndef INCLUDE_PARSECOMMANDLINE_OPTS_H_ -#define INCLUDE_PARSECOMMANDLINE_OPTS_H_ +#ifndef PARSE_CMD_OPTIONS_H_ +#define PARSE_CMD_OPTIONS_H_ #include #include +#include #include #include #include #include -class ParseCommandLineOpts +class parseCMDOptions { public: - ParseCommandLineOpts(int &_argc, char **argv) + parseCMDOptions(int &_argc, char **argv) + : argc(_argc) { - argc = _argc; for (int i = 1; i < argc; ++i) - tokens.push_back(std::string(argv[i])); + { + tokens.emplace_back(argv[i]); + } } std::string @@ -29,11 +32,8 @@ class ParseCommandLineOpts if (cmdOptionExists("-i")) { parameters_filename = getCmdOption("-i"); - if (dealii::Utilities::MPI::this_mpi_process(MPI_COMM_WORLD) == 0) - { - std::cout << "Using the input parameter file: " << parameters_filename - << std::endl; - } + conditionalOStreams::pout_base + << "Using the input parameter file: " << parameters_filename << '\n'; } else { @@ -49,14 +49,13 @@ class ParseCommandLineOpts std::ifstream ifs_prm(parameters_filename); std::ifstream ifs_in("parameters.in"); if (!ifs_prm && ifs_in) - throw("The previous extension .in for the parameters file is no " - "longer accepted. Please rename parameters.in as " - "parameters.prm"); - if (dealii::Utilities::MPI::this_mpi_process(MPI_COMM_WORLD) == 0) { - std::cout << "Using the input parameter file: " << parameters_filename - << std::endl; + throw("The previous extension .in for the parameters file is no " + "longer accepted. Please rename parameters.in as " + "parameters.prm"); } + conditionalOStreams::pout_base + << "Using the input parameter file: " << parameters_filename << '\n'; } else { @@ -70,13 +69,13 @@ class ParseCommandLineOpts std::ifstream ifs_prm(parameters_filename); std::ifstream ifs_in("parameters.in"); if (!ifs_prm && ifs_in) - throw("The previous extension .in for the parameters file is no longer " - "accepted. Please rename parameters.in as parameters.prm"); - if (dealii::Utilities::MPI::this_mpi_process(MPI_COMM_WORLD) == 0) { - std::cout << "Using the input parameter file: " << parameters_filename - << std::endl; + throw("The previous extension .in for the parameters file is no longer " + "accepted. Please rename parameters.in as parameters.prm"); } + + conditionalOStreams::pout_base + << "Using the input parameter file: " << parameters_filename << '\n'; } else { @@ -91,7 +90,7 @@ class ParseCommandLineOpts int argc; std::vector tokens; - const std::string & + [[nodiscard]] const std::string & getCmdOption(const std::string &option) const { std::vector::const_iterator itr; @@ -100,11 +99,11 @@ class ParseCommandLineOpts { return *itr; } - static const std::string empty_string(""); + static const std::string empty_string; return empty_string; } - bool + [[nodiscard]] bool cmdOptionExists(const std::string &option) const { return std::find(tokens.begin(), tokens.end(), option) != tokens.end(); diff --git a/include/core/refinement/AdaptiveRefinement.h b/include/core/refinement/AdaptiveRefinement.h index fad85f84..0ab00343 100644 --- a/include/core/refinement/AdaptiveRefinement.h +++ b/include/core/refinement/AdaptiveRefinement.h @@ -10,7 +10,6 @@ #include #include -#include #include using namespace dealii; @@ -28,7 +27,6 @@ class AdaptiveRefinement AdaptiveRefinement( const userInputParameters &_userInputs, parallel::distributed::Triangulation &_triangulation, - std::vector> &_fields, std::vector *> &_solutionSet, std::vector &triangulation; - std::vector> &fields; - std::vector *> &solutionSet; std::vector +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include + +/** + * Forward declaration for user-implemented PDE class. + */ +template +class customPDE; + +const unsigned int degree = 2; +const unsigned int dimension = 3; + +template +class LaplaceProblem +{ +public: + using SystemMatrixType = customPDE; + using LevelMatrixType = customPDE; + + LaplaceProblem(); + void + run(); + +private: + void + setup_system(); + void + assemble_rhs(); + void + solve(); + void + output_results(unsigned int cycle) const; + + dealii::parallel::distributed::Triangulation triangulation; + + const dealii::FE_Q fe; + dealii::DoFHandler dof_handler; + + dealii::MappingQ1 mapping; + + dealii::AffineConstraints constraints; + + SystemMatrixType system_matrix; + + dealii::MGConstrainedDoFs mg_constrained_dofs; + + dealii::MGLevelObject mg_matrices; + + dealii::LinearAlgebra::distributed::Vector solution; + dealii::LinearAlgebra::distributed::Vector system_rhs; + + double setup_time {}; +}; + +template +LaplaceProblem::LaplaceProblem() + : triangulation( + MPI_COMM_WORLD, + dealii::Triangulation::limit_level_difference_at_vertices, + dealii::parallel::distributed::Triangulation::construct_multigrid_hierarchy) + , fe(degree) + , dof_handler(triangulation) +{} + +template +void +LaplaceProblem::setup_system() +{ + Timer time; + setup_time = 0; + + system_matrix.clear(); + mg_matrices.clear_elements(); + + dof_handler.distribute_dofs(fe); + dof_handler.distribute_mg_dofs(); + + conditionalOStreams::pout_base + << "Number of degrees of freedom: " << dof_handler.n_dofs() << std::endl; + + constraints.clear(); + constraints.reinit(dof_handler.locally_owned_dofs(), + dealii::DoFTools::extract_locally_relevant_dofs(dof_handler)); + dealii::DoFTools::make_hanging_node_constraints(dof_handler, constraints); + dealii::VectorTools::interpolate_boundary_values(mapping, + dof_handler, + 0, + dealii::Functions::ZeroFunction(), + constraints); + constraints.close(); + + setup_time += time.wall_time(); + conditionalOStreams::pout_base << "Distribute DoFs & B.C. (CPU/wall) " + << time.cpu_time() << "s/" << time.wall_time() << 's' + << '\n'; + time.restart(); + + typename dealii::MatrixFree::AdditionalData additional_data; + additional_data.tasks_parallel_scheme = + dealii::MatrixFree::AdditionalData::none; + additional_data.mapping_update_flags = + (update_gradients | update_JxW_values | update_quadrature_points); + std::shared_ptr> system_mf_storage( + new dealii::MatrixFree()); + system_mf_storage->reinit(mapping, + dof_handler, + constraints, + dealii::QGauss<1>(fe.degree + 1), + additional_data); + system_matrix.initialize(system_mf_storage); + + system_matrix.initialize_dof_vector(solution); + system_matrix.initialize_dof_vector(system_rhs); + + setup_time += time.wall_time(); + conditionalOStreams::pout_base << "Setup matrix-free system (CPU/wall) " + << time.cpu_time() << "s/" << time.wall_time() << 's' + << '\n'; + time.restart(); + + const unsigned int nlevels = triangulation.n_global_levels(); + mg_matrices.resize(0, nlevels - 1); + + const std::set dirichlet_boundary_ids = {0}; + mg_constrained_dofs.initialize(dof_handler); + mg_constrained_dofs.make_zero_boundary_constraints(dof_handler, dirichlet_boundary_ids); + + for (unsigned int level = 0; level < nlevels; ++level) + { + dealii::AffineConstraints level_constraints( + dof_handler.locally_owned_mg_dofs(level), + dealii::DoFTools::extract_locally_relevant_level_dofs(dof_handler, level)); + for (const dealii::types::global_dof_index dof_index : + mg_constrained_dofs.get_boundary_indices(level)) + { + level_constraints.constrain_dof_to_zero(dof_index); + } + level_constraints.close(); + + typename dealii::MatrixFree::AdditionalData additional_data; + additional_data.tasks_parallel_scheme = + dealii::MatrixFree::AdditionalData::none; + additional_data.mapping_update_flags = + (update_gradients | update_JxW_values | update_quadrature_points); + additional_data.mg_level = level; + std::shared_ptr> mg_mf_storage_level = + std::make_shared>(); + mg_mf_storage_level->reinit(mapping, + dof_handler, + level_constraints, + dealii::QGauss<1>(fe.degree + 1), + additional_data); + + mg_matrices[level].initialize(mg_mf_storage_level, mg_constrained_dofs, level); + } + + setup_time += time.wall_time(); + conditionalOStreams::pout_base << "Setup matrix-free levels (CPU/wall) " + << time.cpu_time() << "s/" << time.wall_time() << 's' + << '\n'; +} + +template +void +LaplaceProblem::assemble_rhs() +{ + Timer time; + + system_rhs = 0; + dealii::FEEvaluation phi(*system_matrix.get_matrix_free()); + for (unsigned int cell = 0; cell < system_matrix.get_matrix_free()->n_cell_batches(); + ++cell) + { + phi.reinit(cell); + for (const unsigned int q : phi.quadrature_point_indices()) + { + phi.submit_value(dealii::make_vectorized_array(1.0), q); + } + phi.integrate(dealii::EvaluationFlags::EvaluationFlags::values); + phi.distribute_local_to_global(system_rhs); + } + system_rhs.compress(dealii::VectorOperation::add); + + setup_time += time.wall_time(); + conditionalOStreams::pout_base << "Assemble right hand side (CPU/wall) " + << time.cpu_time() << "s/" << time.wall_time() << 's' + << '\n'; +} + +template +void +LaplaceProblem::solve() +{ + Timer time; + dealii::MGTransferMatrixFree mg_transfer(mg_constrained_dofs); + mg_transfer.build(dof_handler); + setup_time += time.wall_time(); + conditionalOStreams::pout_base << "MG build transfer time (CPU/wall) " + << time.cpu_time() << "s/" << time.wall_time() << "s\n"; + time.restart(); + + using SmootherType = + dealii::PreconditionChebyshev>; + dealii::mg::SmootherRelaxation> + mg_smoother; + dealii::MGLevelObject smoother_data; + smoother_data.resize(0, triangulation.n_global_levels() - 1); + for (unsigned int level = 0; level < triangulation.n_global_levels(); ++level) + { + if (level > 0) + { + smoother_data[level].smoothing_range = 15.; + smoother_data[level].degree = 5; + smoother_data[level].eig_cg_n_iterations = 10; + } + else + { + smoother_data[0].smoothing_range = 1e-3; + smoother_data[0].degree = dealii::numbers::invalid_unsigned_int; + smoother_data[0].eig_cg_n_iterations = mg_matrices[0].m(); + } + mg_matrices[level].compute_diagonal(); + smoother_data[level].preconditioner = + mg_matrices[level].get_matrix_diagonal_inverse(); + } + mg_smoother.initialize(mg_matrices, smoother_data); + + dealii::MGCoarseGridApplySmoother> + mg_coarse; + mg_coarse.initialize(mg_smoother); + + dealii::mg::Matrix> mg_matrix( + mg_matrices); + + dealii::MGLevelObject> + mg_interface_matrices; + mg_interface_matrices.resize(0, triangulation.n_global_levels() - 1); + for (unsigned int level = 0; level < triangulation.n_global_levels(); ++level) + { + mg_interface_matrices[level].initialize(mg_matrices[level]); + } + dealii::mg::Matrix> mg_interface( + mg_interface_matrices); + + Multigrid> mg(mg_matrix, + mg_coarse, + mg_transfer, + mg_smoother, + mg_smoother); + mg.set_edge_matrices(mg_interface, mg_interface); + + dealii::PreconditionMG, + dealii::MGTransferMatrixFree> + preconditioner(dof_handler, mg, mg_transfer); + + dealii::SolverControl solver_control(100, 1e-12 * system_rhs.l2_norm()); + dealii::SolverCG> cg(solver_control); + setup_time += time.wall_time(); + conditionalOStreams::pout_base << "MG build smoother time (CPU/wall) " + << time.cpu_time() << "s/" << time.wall_time() << "s\n"; + conditionalOStreams::pout_base << "Total setup time (wall) " << setup_time + << "s\n"; + + time.reset(); + time.start(); + constraints.set_zero(solution); + cg.solve(system_matrix, solution, system_rhs, preconditioner); + + constraints.distribute(solution); + + conditionalOStreams::pout_base + << "Time solve (" << solver_control.last_step() << " iterations)" + << (solver_control.last_step() < 10 ? " " : " ") << "(CPU/wall) " << time.cpu_time() + << "s/" << time.wall_time() << "s\n"; +} + +template +void +LaplaceProblem::output_results(unsigned int cycle) const +{ + Timer time; + if (triangulation.n_global_active_cells() > 1000000) + { + return; + } + + dealii::DataOut data_out; + + solution.update_ghost_values(); + data_out.attach_dof_handler(dof_handler); + data_out.add_data_vector(solution, "solution"); + data_out.build_patches(mapping); + + dealii::DataOutBase::VtkFlags flags; + flags.compression_level = dealii::DataOutBase::CompressionLevel::best_speed; + data_out.set_flags(flags); + data_out.write_vtu_with_pvtu_record("./", "solution", cycle, MPI_COMM_WORLD, 3); + + conditionalOStreams::pout_base << "Time write output (CPU/wall) " + << time.cpu_time() << "s/" << time.wall_time() << "s\n"; +} + +template +void +LaplaceProblem::run() +{ + const unsigned int n_vect_doubles = dealii::VectorizedArray::size(); + const unsigned int n_vect_bits = 8 * sizeof(double) * n_vect_doubles; + + conditionalOStreams::pout_base + << "Vectorization over " << n_vect_doubles << " doubles = " << n_vect_bits + << " bits (" << Utilities::System::get_current_vectorization_level() << ')' << '\n'; + + dealii::GridGenerator::hyper_cube(triangulation, 0.0, 1.0); + triangulation.refine_global(6); + + setup_system(); + assemble_rhs(); + solve(); + output_results(0); +} + +#endif \ No newline at end of file diff --git a/include/core/typeDefs.h b/include/core/typeDefs.h deleted file mode 100644 index 73211048..00000000 --- a/include/core/typeDefs.h +++ /dev/null @@ -1,50 +0,0 @@ -/* - * typeDefs.h - * - * Created on: Feb 24, 2017 - * Author: stephendewitt - */ - -// #ifndef INCLUDE_TYPEDEFS_H_ -// #define INCLUDE_TYPEDEFS_H_ - -// define FE system types -#ifndef typeScalar -using typeScalar = dealii::FEEvaluation; -#endif -#ifndef typeVector -using typeVector = dealii::FEEvaluation; -#endif -// define data value types -#ifndef scalarvalueType -using scalarvalueType = dealii::VectorizedArray; -#endif -#ifndef vectorvalueType -using vectorvalueType = dealii::Tensor<1, dim, dealii::VectorizedArray>; -#endif -#if problemDIM == 1 -# ifndef scalargradType -using scalargradType = dealii::VectorizedArray; -# endif -# ifndef vectorgradType -using vectorgradType = dealii::VectorizedArray; -# endif -# ifndef vectorhessType -using vectorhessType = dealii::VectorizedArray; -# endif -#else -# ifndef scalargradType -using scalargradType = dealii::Tensor<1, dim, dealii::VectorizedArray>; -# endif -# ifndef scalarhessType -using scalarhessType = dealii::Tensor<2, dim, dealii::VectorizedArray>; -# endif -# ifndef vectorgradType -using vectorgradType = dealii::Tensor<2, dim, dealii::VectorizedArray>; -# endif -# ifndef vectorhessType -using vectorhessType = dealii::Tensor<3, dim, dealii::VectorizedArray>; -# endif -#endif - -// #endif /* INCLUDE_TYPEDEFS_H_ */ diff --git a/include/core/varTypeEnums.h b/include/core/typeEnums.h similarity index 56% rename from include/core/varTypeEnums.h rename to include/core/typeEnums.h index 6bf238cc..8452b832 100644 --- a/include/core/varTypeEnums.h +++ b/include/core/typeEnums.h @@ -1,5 +1,5 @@ -#ifndef INCLUDE_VARTYPEENUMS_H_ -#define INCLUDE_VARTYPEENUMS_H_ +#ifndef INCLUDE_TYPEENUMS_H_ +#define INCLUDE_TYPEENUMS_H_ enum fieldType { @@ -17,4 +17,12 @@ enum PDEType AUXILIARY }; +enum solveType +{ + EXPLICIT_RHS, + NONEXPLICIT_RHS, + NONEXPLICIT_LHS, + POSTPROCESS +}; + #endif diff --git a/include/core/userInputParameters.h b/include/core/userInputParameters.h index 5a988658..fe39424d 100644 --- a/include/core/userInputParameters.h +++ b/include/core/userInputParameters.h @@ -31,7 +31,7 @@ using InputVariant = boost::variant, dealii::Tensor<2, dim>, - dealii::Tensor<2, 2 * dim - 1 + dim / 3>>; + dealii::Tensor<2, (2 * dim) - 1 + (dim / 3)>>; enum elasticityModel { @@ -163,7 +163,7 @@ class userInputParameters * * \param constant_name Name of the constant to retrieve. */ - [[nodiscard]] dealii::Tensor<2, 2 * dim - 1 + dim / 3> + [[nodiscard]] dealii::Tensor<2, (2 * dim) - 1 + (dim / 3)> get_model_constant_elasticity_tensor(const std::string &constant_name) const { Assert(model_constants.find(constant_name) != model_constants.end(), @@ -172,7 +172,7 @@ class userInputParameters "customPDE.h. The constant that you attempted to access was " + constant_name + ".")); - return boost::get>( + return boost::get>( model_constants.at(constant_name)); }; @@ -264,13 +264,13 @@ class userInputParameters const AttributesList &pp_attributes; // Variables needed to calculate the RHS - unsigned int num_var_explicit_RHS, num_var_nonexplicit_RHS; - std::vector varInfoListExplicitRHS, varInfoListNonexplicitRHS; + unsigned int num_var_explicit_RHS, num_var_nonexplicit_RHS; + std::vector varInfoListExplicitRHS, varInfoListNonexplicitRHS; // Variables needed to calculate the LHS - unsigned int num_var_LHS; - std::vector varInfoListLHS; - std::vector varChangeInfoListLHS; + unsigned int num_var_LHS; + std::vector varInfoListLHS; + std::vector varChangeInfoListLHS; // Variables for loading in initial conditions std::vector load_ICs; @@ -288,8 +288,8 @@ class userInputParameters std::vector integrated_field_indices; // Variable and residual info - std::vector pp_varInfoList; - std::vector pp_baseVarInfoList; + std::vector pp_varInfoList; + std::vector pp_baseVarInfoList; // List of boundary conditions std::vector> BC_list; @@ -304,8 +304,7 @@ class userInputParameters bool evolution_before_nucleation; // Declare later // bool multiple_nuclei_per_order_parameter; - double min_distance_between_nuclei; // Only enforced for nuclei placed during - // the same time step + double min_distance_between_nuclei; double nucleation_order_parameter_cutoff; unsigned int steps_between_nucleation_attempts; double nucleation_start_time; @@ -313,15 +312,13 @@ class userInputParameters // Grain remapping parameters bool grain_remapping_activated; - std::vector variables_for_remapping; // Note: this should be a sorted list + std::vector variables_for_remapping; unsigned int skip_grain_reassignment_steps; double order_parameter_threshold; double buffer_between_grains; - bool load_grain_structure; - std::string load_vtk_file_type; // adding this string to know what type of vtk file you - // want to read, it will be passed to - // initialconditions.cc + bool load_grain_structure; + std::string load_vtk_file_type; double min_radius_for_loading_grains; std::string grain_structure_filename; std::string grain_structure_variable_name; @@ -440,14 +437,12 @@ class userInputParameters InputVariant primitive_model_constant(std::vector &model_constants_strings); - [[nodiscard]] dealii::Tensor<2, 2 * dim - 1 + dim / 3> + [[nodiscard]] dealii::Tensor<2, (2 * dim) - 1 + (dim / 3)> get_Cij_tensor(std::vector elastic_constants, const std::string &elastic_const_symmetry) const; - dealii::Tensor<2, 2 * dim - 1 + dim / 3> - getCIJMatrix(const elasticityModel model, - const std::vector &constants, - dealii::ConditionalOStream &pcout) const; + dealii::Tensor<2, (2 * dim) - 1 + (dim / 3)> + getCIJMatrix(const elasticityModel model, const std::vector &constants) const; // Private nucleation variables std::vector> nucleation_parameters_list; diff --git a/include/core/variableAttributeLoader.h b/include/core/variableAttributeLoader.h index 926ca033..8dc291d3 100644 --- a/include/core/variableAttributeLoader.h +++ b/include/core/variableAttributeLoader.h @@ -1,19 +1,19 @@ #ifndef VARIABLEATTRIBUTELOADER_H #define VARIABLEATTRIBUTELOADER_H -#include +#include #include #include #include -using EvalFlags = dealii::EvaluationFlags::EvaluationFlags; - /** * \brief Class to manage the variable attributes that the user specifies. */ class variableAttributeLoader { public: + using EvalFlags = dealii::EvaluationFlags::EvaluationFlags; + /** * \brief Constructor. */ diff --git a/include/core/variableAttributes.h b/include/core/variableAttributes.h index 5ebd93c3..8eab0889 100644 --- a/include/core/variableAttributes.h +++ b/include/core/variableAttributes.h @@ -3,7 +3,7 @@ #include -#include +#include #include #include #include diff --git a/include/core/variableContainer.h b/include/core/variableContainer.h index f7f311d6..10d402ce 100644 --- a/include/core/variableContainer.h +++ b/include/core/variableContainer.h @@ -1,5 +1,3 @@ -// This class permits the access of a subset of indexed fields and gives an -// error if any non-allowed fields are requested #ifndef VARIBLECONTAINER_H #define VARIBLECONTAINER_H @@ -8,140 +6,213 @@ #include #include -#include - +#include #include - -template +#include + +/** + * \brief This class permits the access of a subset of indexed fields and gives an error + * if any non-allowed fields are requested. + * + * \tparam dim The number of dimensions in the problem. + * \tparam degree The polynomial degree of the shape functions. + * \tparam number Datatype to use for `dealii::VectorizedArray`. Either + * double or float. + */ +template class variableContainer { public: -#include - // Constructors - - // Standard contructor, used for most situations - variableContainer(const dealii::MatrixFree &data, - const std::vector &_varInfoList, - const std::vector &_varChangeInfoList); - - variableContainer(const dealii::MatrixFree &data, - const std::vector &_varInfoList); - // Nonstandard constructor, used when only one index of "data" should be used, - // use with care! - variableContainer(const dealii::MatrixFree &data, - const std::vector &_varInfoList, - const unsigned int &fixed_index); - - // Methods to get the value/grad/hess in the residual method (this is how the - // user gets these values in equations.h) - T + using EvalFlags = dealii::EvaluationFlags::EvaluationFlags; + + /** + * \brief Constructor. + */ + variableContainer(const dealii::MatrixFree &data); + + // Methods to get the value/grad/hess in the residual method + dealii::VectorizedArray get_scalar_value(unsigned int global_variable_index) const; - dealii::Tensor<1, dim, T> + dealii::Tensor<1, dim, dealii::VectorizedArray> get_scalar_gradient(unsigned int global_variable_index) const; - dealii::Tensor<2, dim, T> + dealii::Tensor<2, dim, dealii::VectorizedArray> get_scalar_hessian(unsigned int global_variable_index) const; - dealii::Tensor<1, dim, T> + dealii::Tensor<1, dim, dealii::VectorizedArray> get_vector_value(unsigned int global_variable_index) const; - dealii::Tensor<2, dim, T> + dealii::Tensor<2, dim, dealii::VectorizedArray> get_vector_gradient(unsigned int global_variable_index) const; - dealii::Tensor<3, dim, T> + dealii::Tensor<3, dim, dealii::VectorizedArray> get_vector_hessian(unsigned int global_variable_index) const; - T - get_change_in_scalar_value(unsigned int global_variable_index) const; - dealii::Tensor<1, dim, T> - get_change_in_scalar_gradient(unsigned int global_variable_index) const; - dealii::Tensor<2, dim, T> - get_change_in_scalar_hessian(unsigned int global_variable_index) const; - dealii::Tensor<1, dim, T> - get_change_in_vector_value(unsigned int global_variable_index) const; - dealii::Tensor<2, dim, T> - get_change_in_vector_gradient(unsigned int global_variable_index) const; - dealii::Tensor<3, dim, T> - get_change_in_vector_hessian(unsigned int global_variable_index) const; - - // Methods to set the value residual and the gradient residual (this is how - // the user sets these values in equations.h) + // Methods to set the value residual and the gradient residual void - set_scalar_value_term_RHS(unsigned int global_variable_index, T val); + set_scalar_value_term(unsigned int global_variable_index, + dealii::VectorizedArray val); void - set_scalar_gradient_term_RHS(unsigned int global_variable_index, - dealii::Tensor<1, dim, T> grad); + set_scalar_gradient_term(unsigned int global_variable_index, + dealii::Tensor<1, dim, dealii::VectorizedArray> grad); void - set_vector_value_term_RHS(unsigned int global_variable_index, - dealii::Tensor<1, dim, T> val); + set_vector_value_term(unsigned int global_variable_index, + dealii::Tensor<1, dim, dealii::VectorizedArray> val); void - set_vector_gradient_term_RHS(unsigned int global_variable_index, - dealii::Tensor<2, dim, T> grad); + set_vector_gradient_term(unsigned int global_variable_index, + dealii::Tensor<2, dim, dealii::VectorizedArray> grad); + /** + * \brief Initialize, read DOFs, and set evaulation flags for each variable. + */ void - set_scalar_value_term_LHS(unsigned int global_variable_index, T val); - void - set_scalar_gradient_term_LHS(unsigned int global_variable_index, - dealii::Tensor<1, dim, T> grad); - void - set_vector_value_term_LHS(unsigned int global_variable_index, - dealii::Tensor<1, dim, T> val); - void - set_vector_gradient_term_LHS(unsigned int global_variable_index, - dealii::Tensor<2, dim, T> grad); + reinit_and_eval(const dealii::LinearAlgebra::distributed::Vector &src, + unsigned int cell); - // Initialize, read DOFs, and set evaulation flags for each variable - void - reinit_and_eval( - const std::vector *> &src, - unsigned int cell); + /** + * \brief Integrate the residuals and distribute from local to global. + */ void - reinit_and_eval_change_in_solution( - const dealii::LinearAlgebra::distributed::Vector &src, - unsigned int cell, - unsigned int var_being_solved); + integrate_and_distribute(dealii::LinearAlgebra::distributed::Vector &dst); - // Only initialize the FEEvaluation object for each variable (used for - // post-processing) + /** + * \brief TODO: Add comments + */ void - reinit(unsigned int cell); - - // Integrate the residuals and distribute from local to global - void - integrate_and_distribute( - std::vector *> &dst); + eval_local_operator( + const std::function> &)> + &func, + dealii::LinearAlgebra::distributed::Vector &dst, + const LinearAlgebra::distributed::Vector &src, + const std::pair &cell_range) + { + for (unsigned int cell = cell_range.first; cell < cell_range.second; ++cell) + { + // Initialize, read DOFs, and set evaulation flags for each variable + reinit_and_eval(src, cell); + + // Number of quadrature points + unsigned int n_q_points = get_num_q_points(); + + for (unsigned int q = 0; q < n_q_points; ++q) + { + // Set the quadrature point + q_point = q; + + // Grab the quadrature point location + Point> q_point_loc = get_q_point_location(); + + // Calculate the residuals + func(*this, q_point_loc); + } + + // Integrate and add to global vector dst + integrate_and_distribute(dst); + } + }; + + /** + * \brief TODO: Add comments + */ void - integrate_and_distribute_change_in_solution_LHS( - dealii::LinearAlgebra::distributed::Vector &dst, - const unsigned int var_being_solved); + eval_local_diagonal( + const std::function> &)> + &func, + dealii::LinearAlgebra::distributed::Vector &dst, + const std::pair &cell_range, + const unsigned int &global_var_index) + { + const auto &var_info = varInfoList[global_var_index]; + + AssertThrow(var_info.var_type != fieldType::VECTOR, + FeatureNotImplemented("Vector multigrid")); + + auto *scalar_FEEval_ptr = scalar_vars_map.at(global_var_index).get(); + + n_dofs_per_cell = scalar_FEEval_ptr->dofs_per_cell; + diagonal = std::make_unique>>( + n_dofs_per_cell); + + for (unsigned int cell = cell_range.first; cell < cell_range.second; ++cell) + { + scalar_FEEval_ptr->reinit(cell); + + for (unsigned int i = 0; i < n_dofs_per_cell; ++i) + { + for (unsigned int j = 0; j < n_dofs_per_cell; ++j) + { + scalar_FEEval_ptr->submit_dof_value(dealii::VectorizedArray(), j); + } + scalar_FEEval_ptr->submit_dof_value(dealii::make_vectorized_array( + 1.0), + i); + + scalar_FEEval_ptr->evaluate(var_info.evaluation_flags); + + // Number of quadrature points + unsigned int n_q_points = get_num_q_points(); + + for (unsigned int q = 0; q < n_q_points; ++q) + { + // Set the quadrature point + q_point = q; + + // Grab the quadrature point location + dealii::Point> q_point_loc = + get_q_point_location(); + + // Calculate the residuals + func(*this, q_point_loc); + } + + scalar_FEEval_ptr->integrate(var_info.residual_flags); + (*diagonal)[i] = scalar_FEEval_ptr->get_dof_value(i); + } + + for (unsigned int i = 0; i < n_dofs_per_cell; ++i) + { + scalar_FEEval_ptr->submit_dof_value((*diagonal)[i], i); + } + scalar_FEEval_ptr->distribute_local_to_global(dst); + } + }; // The quadrature point index, a method to get the number of quadrature points // per cell, and a method to get the xyz coordinates for the quadrature point - unsigned int q_point; + unsigned int q_point = 0; - unsigned int + // Number of DoFs per cell + unsigned int n_dofs_per_cell; + + [[nodiscard]] unsigned int get_num_q_points() const; - dealii::Point + [[nodiscard]] dealii::Point> get_q_point_location() const; private: + /** + * \brief Check that a variable value/gradient/hessians was marked as needed and thus + * properly initialized. + */ + void + access_valid(const unsigned int &global_variable_index, const EvalFlags &flag) const; + // Vectors of the actual FEEvaluation objects for each active variable, split // into scalar variables and vector variables for type reasons - using scalar_FEEval = dealii::FEEvaluation; - using vector_FEEval = dealii::FEEvaluation; - - boost::unordered_map> scalar_vars_map; - boost::unordered_map> vector_vars_map; + using scalar_FEEval = dealii::FEEvaluation; + using vector_FEEval = dealii::FEEvaluation; - boost::unordered_map> - scalar_change_in_vars_map; - boost::unordered_map> - vector_change_in_vars_map; + std::unordered_map> scalar_vars_map; + std::unordered_map> vector_vars_map; // Object containing some information about each variable (indices, whether // the val/grad/hess is needed, etc) - std::vector varInfoList; - std::vector varChangeInfoList; + std::vector varInfoList; // The number of variables unsigned int num_var; + + // Diagonal matrix that is used for preconditioning + std::unique_ptr>> diagonal; }; #endif diff --git a/src/core/CMakeLists.txt b/src/core/CMakeLists.txt index 2c8d5b1b..50ba40a3 100644 --- a/src/core/CMakeLists.txt +++ b/src/core/CMakeLists.txt @@ -5,15 +5,10 @@ list(APPEND PRISMS_PF_SOURCE_FILES ${CMAKE_CURRENT_SOURCE_DIR}/boundary_conditions/vectorBCFunction.cc ${CMAKE_CURRENT_SOURCE_DIR}/initial_conditions/initialConditions.cc ${CMAKE_CURRENT_SOURCE_DIR}/postprocessing/computeIntegral.cc - ${CMAKE_CURRENT_SOURCE_DIR}/postprocessing/postprocessor.cc ${CMAKE_CURRENT_SOURCE_DIR}/refinement/AdaptiveRefinement.cc - ${CMAKE_CURRENT_SOURCE_DIR}/solvers/computeLHS.cc - ${CMAKE_CURRENT_SOURCE_DIR}/solvers/computeRHS.cc - ${CMAKE_CURRENT_SOURCE_DIR}/solvers/setNonlinearEqInitialGuess.cc ${CMAKE_CURRENT_SOURCE_DIR}/solvers/solve.cc ${CMAKE_CURRENT_SOURCE_DIR}/solvers/solveIncrement.cc ${CMAKE_CURRENT_SOURCE_DIR}/solvers/SolverParameters.cc - ${CMAKE_CURRENT_SOURCE_DIR}/buildFields.cc ${CMAKE_CURRENT_SOURCE_DIR}/checkpoint.cc ${CMAKE_CURRENT_SOURCE_DIR}/init.cc ${CMAKE_CURRENT_SOURCE_DIR}/inputFileReader.cc diff --git a/src/core/boundary_conditions/boundaryConditions.cc b/src/core/boundary_conditions/boundaryConditions.cc index 90b55028..9fe08c7c 100644 --- a/src/core/boundary_conditions/boundaryConditions.cc +++ b/src/core/boundary_conditions/boundaryConditions.cc @@ -1,8 +1,8 @@ -// methods to apply boundary conditons - +#include #include #include #include +#include #include // ================================================================================= @@ -236,7 +236,8 @@ MatrixFreePDE::setPeriodicity() } triangulation.add_periodicity(periodicity_vector); - pcout << "periodic facepairs: " << periodicity_vector.size() << std::endl; + conditionalOStreams::pout_base << "periodic facepairs: " << periodicity_vector.size() + << std::endl; } // Set constraints to enforce periodic boundary conditions @@ -313,20 +314,4 @@ MatrixFreePDE::set_rigid_body_mode_constraints( } } -template class MatrixFreePDE<2, 1>; -template class MatrixFreePDE<3, 1>; - -template class MatrixFreePDE<2, 2>; -template class MatrixFreePDE<3, 2>; - -template class MatrixFreePDE<3, 3>; -template class MatrixFreePDE<2, 3>; - -template class MatrixFreePDE<3, 4>; -template class MatrixFreePDE<2, 4>; - -template class MatrixFreePDE<3, 5>; -template class MatrixFreePDE<2, 5>; - -template class MatrixFreePDE<3, 6>; -template class MatrixFreePDE<2, 6>; \ No newline at end of file +INSTANTIATE_BI_TEMPLATE(MatrixFreePDE) \ No newline at end of file diff --git a/src/core/boundary_conditions/markBoundaries.cc b/src/core/boundary_conditions/markBoundaries.cc index 4723df5f..d4ec0618 100644 --- a/src/core/boundary_conditions/markBoundaries.cc +++ b/src/core/boundary_conditions/markBoundaries.cc @@ -1,5 +1,4 @@ -// methods to mark boundaries - +#include #include #include @@ -32,20 +31,4 @@ MatrixFreePDE::markBoundaries( } } -template class MatrixFreePDE<2, 1>; -template class MatrixFreePDE<3, 1>; - -template class MatrixFreePDE<2, 2>; -template class MatrixFreePDE<3, 2>; - -template class MatrixFreePDE<3, 3>; -template class MatrixFreePDE<2, 3>; - -template class MatrixFreePDE<3, 4>; -template class MatrixFreePDE<2, 4>; - -template class MatrixFreePDE<3, 5>; -template class MatrixFreePDE<2, 5>; - -template class MatrixFreePDE<3, 6>; -template class MatrixFreePDE<2, 6>; +INSTANTIATE_BI_TEMPLATE(MatrixFreePDE) \ No newline at end of file diff --git a/src/core/buildFields.cc b/src/core/buildFields.cc deleted file mode 100644 index 497b7aed..00000000 --- a/src/core/buildFields.cc +++ /dev/null @@ -1,41 +0,0 @@ -/* - * buildFields.cc - * - * Created on: Feb 22, 2017 - * Author: stephendewitt - */ - -// ===================================================================== -// FUNCTION TO BUILD THE VECTOR OF FIELDS -// ===================================================================== - -#include - -template -void -MatrixFreePDE::buildFields() -{ - // Build each of the fields in the system - for (const auto &[index, variable] : var_attributes) - { - fields.push_back(Field(variable.var_type, variable.eq_type, variable.name)); - } -} - -template class MatrixFreePDE<2, 1>; -template class MatrixFreePDE<3, 1>; - -template class MatrixFreePDE<2, 2>; -template class MatrixFreePDE<3, 2>; - -template class MatrixFreePDE<3, 3>; -template class MatrixFreePDE<2, 3>; - -template class MatrixFreePDE<3, 4>; -template class MatrixFreePDE<2, 4>; - -template class MatrixFreePDE<3, 5>; -template class MatrixFreePDE<2, 5>; - -template class MatrixFreePDE<3, 6>; -template class MatrixFreePDE<2, 6>; \ No newline at end of file diff --git a/src/core/checkpoint.cc b/src/core/checkpoint.cc index 096b5b45..ed9d88ea 100644 --- a/src/core/checkpoint.cc +++ b/src/core/checkpoint.cc @@ -1,6 +1,8 @@ #include #include +#include +#include #include #include #include @@ -15,7 +17,6 @@ template void MatrixFreePDE::save_checkpoint() { - computing_timer.enter_subsection("matrixFreePDE: save_checkpoint"); unsigned int my_id = Utilities::MPI::this_mpi_process(MPI_COMM_WORLD); if (my_id == 0) @@ -118,8 +119,7 @@ MatrixFreePDE::save_checkpoint() time_info_file.close(); } - pcout << "*** Checkpoint created! ***\n\n"; - computing_timer.leave_subsection("matrixFreePDE: save_checkpoint"); + conditionalOStreams::pout_base << "*** Checkpoint created! ***\n\n"; } // Load from a previously created checkpoint @@ -132,7 +132,7 @@ MatrixFreePDE::load_checkpoint_triangulation() verify_checkpoint_file_exists("restart.mesh"); verify_checkpoint_file_exists("restart.mesh.info"); - pcout << "\n*** Resuming from a checkpoint! ***\n\n"; + conditionalOStreams::pout_base << "\n*** Resuming from a checkpoint! ***\n\n"; try { @@ -274,20 +274,4 @@ MatrixFreePDE::verify_checkpoint_file_exists(const std::string &fil } } -template class MatrixFreePDE<2, 1>; -template class MatrixFreePDE<3, 1>; - -template class MatrixFreePDE<2, 2>; -template class MatrixFreePDE<3, 2>; - -template class MatrixFreePDE<3, 3>; -template class MatrixFreePDE<2, 3>; - -template class MatrixFreePDE<3, 4>; -template class MatrixFreePDE<2, 4>; - -template class MatrixFreePDE<3, 5>; -template class MatrixFreePDE<2, 5>; - -template class MatrixFreePDE<3, 6>; -template class MatrixFreePDE<2, 6>; \ No newline at end of file +INSTANTIATE_BI_TEMPLATE(MatrixFreePDE) \ No newline at end of file diff --git a/src/core/init.cc b/src/core/init.cc index f4fe66c8..054da32a 100644 --- a/src/core/init.cc +++ b/src/core/init.cc @@ -1,8 +1,8 @@ -// init() method for MatrixFreePDE class - #include +#include #include +#include #include // populate with fields and setup matrix free system @@ -10,11 +10,9 @@ template void MatrixFreePDE::init() { - computing_timer.enter_subsection("matrixFreePDE: initialization"); - // creating mesh - pcout << "creating problem mesh...\n"; + conditionalOStreams::pout_base << "creating problem mesh...\n"; // Create the coarse mesh and mark the boundaries create_triangulation(triangulation); @@ -37,197 +35,200 @@ MatrixFreePDE::init() // elements if (dim < 3) { - pcout << "problem dimensions: " << userInputs.domain_size[0] << "x" - << userInputs.domain_size[1] << "\n"; + conditionalOStreams::pout_base + << "problem dimensions: " << userInputs.domain_size[0] << "x" + << userInputs.domain_size[1] << "\n"; } else { - pcout << "problem dimensions: " << userInputs.domain_size[0] << "x" - << userInputs.domain_size[1] << "x" << userInputs.domain_size[2] << "\n"; + conditionalOStreams::pout_base + << "problem dimensions: " << userInputs.domain_size[0] << "x" + << userInputs.domain_size[1] << "x" << userInputs.domain_size[2] << "\n"; } - pcout << "number of elements: " << triangulation.n_global_active_cells() << "\n\n"; + conditionalOStreams::pout_base + << "number of elements: " << triangulation.n_global_active_cells() << "\n\n"; // Setup system - pcout << "initializing matrix free object\n"; + conditionalOStreams::pout_base << "initializing matrix free object\n"; totalDOFs = 0; - for (auto &field : fields) - { - currentFieldIndex = field.index; - - char buffer[100]; - - // print to std::out - std::string var_type; - if (field.pdetype == EXPLICIT_TIME_DEPENDENT) - { - var_type = "EXPLICIT_TIME_DEPENDENT"; - } - else if (field.pdetype == IMPLICIT_TIME_DEPENDENT) - { - var_type = "IMPLICIT_TIME_DEPENDENT"; - } - else if (field.pdetype == TIME_INDEPENDENT) - { - var_type = "TIME_INDEPENDENT"; - } - else if (field.pdetype == AUXILIARY) - { - var_type = "AUXILIARY"; - } - - snprintf(buffer, - sizeof(buffer), - "initializing finite element space P^%u for %9s:%6s field '%s'\n", - degree, - var_type.c_str(), - (field.type == SCALAR ? "SCALAR" : "VECTOR"), - field.name.c_str()); - pcout << buffer; - - // Check if any time dependent fields present - if (field.pdetype == EXPLICIT_TIME_DEPENDENT) - { - isTimeDependentBVP = true; - hasExplicitEquation = true; - } - else if (field.pdetype == IMPLICIT_TIME_DEPENDENT) - { - isTimeDependentBVP = true; - hasNonExplicitEquation = true; - std::cerr << "PRISMS-PF Error: IMPLICIT_TIME_DEPENDENT equation " - "types are not currently supported\n"; - abort(); - } - else if (field.pdetype == AUXILIARY) - { - hasNonExplicitEquation = true; - } - else if (field.pdetype == TIME_INDEPENDENT) - { - isEllipticBVP = true; - hasNonExplicitEquation = true; - } - - // create FESystem - FESystem *fe = nullptr; - - if (field.type == SCALAR) - { - fe = new FESystem(FE_Q(QGaussLobatto<1>(degree + 1)), 1); - } - else if (field.type == VECTOR) - { - fe = new FESystem(FE_Q(QGaussLobatto<1>(degree + 1)), dim); - } - else - { - pcout << "\nmatrixFreePDE.h: unknown field type\n"; - exit(-1); - } - FESet.push_back(fe); - - // distribute DOFs - DoFHandler *dof_handler = nullptr; - - dof_handler = new DoFHandler(triangulation); - dofHandlersSet.push_back(dof_handler); - dofHandlersSet_nonconst.push_back(dof_handler); - - dof_handler->distribute_dofs(*fe); - totalDOFs += dof_handler->n_dofs(); - - // Extract locally_relevant_dofs - IndexSet *locally_relevant_dofs = nullptr; - - locally_relevant_dofs = new IndexSet; - locally_relevant_dofsSet.push_back(locally_relevant_dofs); - locally_relevant_dofsSet_nonconst.push_back(locally_relevant_dofs); - - locally_relevant_dofs->clear(); - DoFTools::extract_locally_relevant_dofs(*dof_handler, *locally_relevant_dofs); - - // Create constraints - AffineConstraints *constraintsDirichlet = nullptr; - AffineConstraints *constraintsOther = nullptr; - - constraintsDirichlet = new AffineConstraints; - constraintsDirichletSet.push_back(constraintsDirichlet); - constraintsDirichletSet_nonconst.push_back(constraintsDirichlet); - constraintsOther = new AffineConstraints; - constraintsOtherSet.push_back(constraintsOther); - constraintsOtherSet_nonconst.push_back(constraintsOther); - valuesDirichletSet.push_back(new std::map); - - constraintsDirichlet->clear(); - constraintsDirichlet->reinit(*locally_relevant_dofs); - constraintsOther->clear(); - constraintsOther->reinit(*locally_relevant_dofs); - - // Get hanging node constraints - DoFTools::make_hanging_node_constraints(*dof_handler, *constraintsOther); - - // Pin solution - if (userInputs.pinned_point.find(currentFieldIndex) != - userInputs.pinned_point.end()) - { - set_rigid_body_mode_constraints(constraintsOther, - dof_handler, - userInputs.pinned_point[currentFieldIndex]); - } - - // Get constraints for periodic BCs - setPeriodicityConstraints(constraintsOther, dof_handler); - - // Check if Dirichlet BCs are used - for (unsigned int i = 0; i < userInputs.BC_list.size(); i++) - { - for (unsigned int direction = 0; direction < 2 * dim; direction++) - { - if (userInputs.BC_list[i].var_BC_type[direction] == DIRICHLET) - { - field.hasDirichletBCs = true; - } - else if (userInputs.BC_list[i].var_BC_type[direction] == - NON_UNIFORM_DIRICHLET) - { - field.hasnonuniformDirichletBCs = true; - } - else if (userInputs.BC_list[i].var_BC_type[direction] == NEUMANN) - { - field.hasNeumannBCs = true; - } - } - } - - // Get constraints for Dirichlet BCs - applyDirichletBCs(); - - constraintsDirichlet->close(); - constraintsOther->close(); - - // Store Dirichlet BC DOF's - valuesDirichletSet[field.index]->clear(); - for (types::global_dof_index i = 0; i < dof_handler->n_dofs(); i++) - { - if (locally_relevant_dofs->is_element(i)) - { - if (constraintsDirichlet->is_constrained(i)) - { - (*valuesDirichletSet[field.index])[i] = - constraintsDirichlet->get_inhomogeneity(i); - } - } - } - - snprintf(buffer, - sizeof(buffer), - "field '%2s' DOF : %u (Constraint DOF : %u)\n", - field.name.c_str(), - dof_handler->n_dofs(), - constraintsDirichlet->n_constraints()); - pcout << buffer; - } - pcout << "total DOF : " << totalDOFs << "\n"; + // for (auto &field : fields) + // { + // currentFieldIndex = field.index; + + // char buffer[100]; + + // // print to std::out + // std::string var_type; + // if (field.pdetype == EXPLICIT_TIME_DEPENDENT) + // { + // var_type = "EXPLICIT_TIME_DEPENDENT"; + // } + // else if (field.pdetype == IMPLICIT_TIME_DEPENDENT) + // { + // var_type = "IMPLICIT_TIME_DEPENDENT"; + // } + // else if (field.pdetype == TIME_INDEPENDENT) + // { + // var_type = "TIME_INDEPENDENT"; + // } + // else if (field.pdetype == AUXILIARY) + // { + // var_type = "AUXILIARY"; + // } + + // snprintf(buffer, + // sizeof(buffer), + // "initializing finite element space P^%u for %9s:%6s field '%s'\n", + // degree, + // var_type.c_str(), + // (field.type == SCALAR ? "SCALAR" : "VECTOR"), + // field.name.c_str()); + // conditionalOStreams::pout_base << buffer; + + // // Check if any time dependent fields present + // if (field.pdetype == EXPLICIT_TIME_DEPENDENT) + // { + // isTimeDependentBVP = true; + // hasExplicitEquation = true; + // } + // else if (field.pdetype == IMPLICIT_TIME_DEPENDENT) + // { + // isTimeDependentBVP = true; + // hasNonExplicitEquation = true; + // std::cerr << "PRISMS-PF Error: IMPLICIT_TIME_DEPENDENT equation " + // "types are not currently supported\n"; + // abort(); + // } + // else if (field.pdetype == AUXILIARY) + // { + // hasNonExplicitEquation = true; + // } + // else if (field.pdetype == TIME_INDEPENDENT) + // { + // isEllipticBVP = true; + // hasNonExplicitEquation = true; + // } + + // // create FESystem + // FESystem *fe = nullptr; + + // if (field.type == SCALAR) + // { + // fe = new FESystem(FE_Q(QGaussLobatto<1>(degree + 1)), 1); + // } + // else if (field.type == VECTOR) + // { + // fe = new FESystem(FE_Q(QGaussLobatto<1>(degree + 1)), dim); + // } + // else + // { + // conditionalOStreams::pout_base << "\nmatrixFreePDE.h: unknown field type\n"; + // exit(-1); + // } + // FESet.push_back(fe); + + // // distribute DOFs + // DoFHandler *dof_handler = nullptr; + + // dof_handler = new DoFHandler(triangulation); + // dofHandlersSet.push_back(dof_handler); + // dofHandlersSet_nonconst.push_back(dof_handler); + + // dof_handler->distribute_dofs(*fe); + // totalDOFs += dof_handler->n_dofs(); + + // // Extract locally_relevant_dofs + // IndexSet *locally_relevant_dofs = nullptr; + + // locally_relevant_dofs = new IndexSet; + // locally_relevant_dofsSet.push_back(locally_relevant_dofs); + // locally_relevant_dofsSet_nonconst.push_back(locally_relevant_dofs); + + // locally_relevant_dofs->clear(); + // DoFTools::extract_locally_relevant_dofs(*dof_handler, *locally_relevant_dofs); + + // // Create constraints + // AffineConstraints *constraintsDirichlet = nullptr; + // AffineConstraints *constraintsOther = nullptr; + + // constraintsDirichlet = new AffineConstraints; + // constraintsDirichletSet.push_back(constraintsDirichlet); + // constraintsDirichletSet_nonconst.push_back(constraintsDirichlet); + // constraintsOther = new AffineConstraints; + // constraintsOtherSet.push_back(constraintsOther); + // constraintsOtherSet_nonconst.push_back(constraintsOther); + // valuesDirichletSet.push_back(new std::map); + + // constraintsDirichlet->clear(); + // constraintsDirichlet->reinit(*locally_relevant_dofs); + // constraintsOther->clear(); + // constraintsOther->reinit(*locally_relevant_dofs); + + // // Get hanging node constraints + // DoFTools::make_hanging_node_constraints(*dof_handler, *constraintsOther); + + // // Pin solution + // if (userInputs.pinned_point.find(currentFieldIndex) != + // userInputs.pinned_point.end()) + // { + // set_rigid_body_mode_constraints(constraintsOther, + // dof_handler, + // userInputs.pinned_point[currentFieldIndex]); + // } + + // // Get constraints for periodic BCs + // setPeriodicityConstraints(constraintsOther, dof_handler); + + // // Check if Dirichlet BCs are used + // for (unsigned int i = 0; i < userInputs.BC_list.size(); i++) + // { + // for (unsigned int direction = 0; direction < 2 * dim; direction++) + // { + // if (userInputs.BC_list[i].var_BC_type[direction] == DIRICHLET) + // { + // field.hasDirichletBCs = true; + // } + // else if (userInputs.BC_list[i].var_BC_type[direction] == + // NON_UNIFORM_DIRICHLET) + // { + // field.hasnonuniformDirichletBCs = true; + // } + // else if (userInputs.BC_list[i].var_BC_type[direction] == NEUMANN) + // { + // field.hasNeumannBCs = true; + // } + // } + // } + + // // Get constraints for Dirichlet BCs + // applyDirichletBCs(); + + // constraintsDirichlet->close(); + // constraintsOther->close(); + + // // Store Dirichlet BC DOF's + // valuesDirichletSet[field.index]->clear(); + // for (types::global_dof_index i = 0; i < dof_handler->n_dofs(); i++) + // { + // if (locally_relevant_dofs->is_element(i)) + // { + // if (constraintsDirichlet->is_constrained(i)) + // { + // (*valuesDirichletSet[field.index])[i] = + // constraintsDirichlet->get_inhomogeneity(i); + // } + // } + // } + + // snprintf(buffer, + // sizeof(buffer), + // "field '%2s' DOF : %u (Constraint DOF : %u)\n", + // field.name.c_str(), + // dof_handler->n_dofs(), + // constraintsDirichlet->n_constraints()); + // conditionalOStreams::pout_base << buffer; + // } + conditionalOStreams::pout_base << "total DOF : " << totalDOFs << "\n"; // Setup the matrix free object typename MatrixFree::AdditionalData additional_data; @@ -247,47 +248,48 @@ MatrixFreePDE::init() bool dU_vector_init = false; // Setup solution vectors - pcout << "initializing parallel::distributed residual and solution vectors\n"; - for (unsigned int fieldIndex = 0; fieldIndex < fields.size(); fieldIndex++) - { - dealii::LinearAlgebra::distributed::Vector *U = nullptr; - dealii::LinearAlgebra::distributed::Vector *R = nullptr; - - U = new dealii::LinearAlgebra::distributed::Vector; - R = new dealii::LinearAlgebra::distributed::Vector; - solutionSet.push_back(U); - residualSet.push_back(R); - matrixFreeObject.initialize_dof_vector(*R, fieldIndex); - *R = 0; - - matrixFreeObject.initialize_dof_vector(*U, fieldIndex); - *U = 0; - - // Initializing temporary dU vector required for implicit solves of the - // elliptic equation. - if (fields[fieldIndex].pdetype == TIME_INDEPENDENT || - fields[fieldIndex].pdetype == IMPLICIT_TIME_DEPENDENT || - (fields[fieldIndex].pdetype == AUXILIARY && - var_attributes.at(fieldIndex).is_nonlinear)) - { - if (fields[fieldIndex].type == SCALAR) - { - if (!dU_scalar_init) - { - matrixFreeObject.initialize_dof_vector(dU_scalar, fieldIndex); - dU_scalar_init = true; - } - } - else - { - if (!dU_vector_init) - { - matrixFreeObject.initialize_dof_vector(dU_vector, fieldIndex); - dU_vector_init = true; - } - } - } - } + conditionalOStreams::pout_base + << "initializing parallel::distributed residual and solution vectors\n"; + // for (unsigned int fieldIndex = 0; fieldIndex < fields.size(); fieldIndex++) + // { + // dealii::LinearAlgebra::distributed::Vector *U = nullptr; + // dealii::LinearAlgebra::distributed::Vector *R = nullptr; + + // U = new dealii::LinearAlgebra::distributed::Vector; + // R = new dealii::LinearAlgebra::distributed::Vector; + // solutionSet.push_back(U); + // residualSet.push_back(R); + // matrixFreeObject.initialize_dof_vector(*R, fieldIndex); + // *R = 0; + + // matrixFreeObject.initialize_dof_vector(*U, fieldIndex); + // *U = 0; + + // // Initializing temporary dU vector required for implicit solves of the + // // elliptic equation. + // if (fields[fieldIndex].pdetype == TIME_INDEPENDENT || + // fields[fieldIndex].pdetype == IMPLICIT_TIME_DEPENDENT || + // (fields[fieldIndex].pdetype == AUXILIARY && + // var_attributes.at(fieldIndex).is_nonlinear)) + // { + // if (fields[fieldIndex].type == SCALAR) + // { + // if (!dU_scalar_init) + // { + // matrixFreeObject.initialize_dof_vector(dU_scalar, fieldIndex); + // dU_scalar_init = true; + // } + // } + // else + // { + // if (!dU_vector_init) + // { + // matrixFreeObject.initialize_dof_vector(dU_vector, fieldIndex); + // dU_vector_init = true; + // } + // } + // } + // } // check if time dependent BVP and compute invM if (isTimeDependentBVP) @@ -310,29 +312,27 @@ MatrixFreePDE::init() // Create new solution transfer sets (needed for the "refine_grid" call, might // be able to move this elsewhere) soltransSet.clear(); - for (unsigned int fieldIndex = 0; fieldIndex < fields.size(); fieldIndex++) - { - soltransSet.push_back( - new parallel::distributed:: - SolutionTransfer>( - *dofHandlersSet_nonconst[fieldIndex])); - } - - // Ghost the solution vectors. Also apply the constraints (if any) on the - // solution vectors - for (unsigned int fieldIndex = 0; fieldIndex < fields.size(); fieldIndex++) - { - constraintsDirichletSet[fieldIndex]->distribute(*solutionSet[fieldIndex]); - constraintsOtherSet[fieldIndex]->distribute(*solutionSet[fieldIndex]); - solutionSet[fieldIndex]->update_ghost_values(); - } + // for (unsigned int fieldIndex = 0; fieldIndex < fields.size(); fieldIndex++) + // { + // soltransSet.push_back( + // new parallel::distributed:: + // SolutionTransfer>( + // *dofHandlersSet_nonconst[fieldIndex])); + // } + + // // Ghost the solution vectors. Also apply the constraints (if any) on the + // // solution vectors + // for (unsigned int fieldIndex = 0; fieldIndex < fields.size(); fieldIndex++) + // { + // constraintsDirichletSet[fieldIndex]->distribute(*solutionSet[fieldIndex]); + // constraintsOtherSet[fieldIndex]->distribute(*solutionSet[fieldIndex]); + // solutionSet[fieldIndex]->update_ghost_values(); + // } // If not resuming from a checkpoint, check and perform adaptive mesh refinement, which // reinitializes the system with the new mesh if (!userInputs.resume_from_checkpoint && userInputs.h_adaptivity == true) { - computing_timer.enter_subsection("matrixFreePDE: AMR"); - unsigned int numDoF_preremesh = totalDOFs; for (unsigned int remesh_index = 0; remesh_index < @@ -347,8 +347,6 @@ MatrixFreePDE::init() } numDoF_preremesh = totalDOFs; } - - computing_timer.leave_subsection("matrixFreePDE: AMR"); } // If resuming from a checkpoint, load the proper starting increment and time @@ -359,8 +357,6 @@ MatrixFreePDE::init() // Once the initial triangulation has been set, compute element volume compute_element_volume(); - - computing_timer.leave_subsection("matrixFreePDE: initialization"); } template @@ -444,20 +440,4 @@ MatrixFreePDE::compute_element_volume() } } -template class MatrixFreePDE<2, 1>; -template class MatrixFreePDE<3, 1>; - -template class MatrixFreePDE<2, 2>; -template class MatrixFreePDE<3, 2>; - -template class MatrixFreePDE<3, 3>; -template class MatrixFreePDE<2, 3>; - -template class MatrixFreePDE<3, 4>; -template class MatrixFreePDE<2, 4>; - -template class MatrixFreePDE<3, 5>; -template class MatrixFreePDE<2, 5>; - -template class MatrixFreePDE<3, 6>; -template class MatrixFreePDE<2, 6>; +INSTANTIATE_BI_TEMPLATE(MatrixFreePDE) \ No newline at end of file diff --git a/src/core/initial_conditions/initialConditions.cc b/src/core/initial_conditions/initialConditions.cc index 923e72b3..19d55496 100644 --- a/src/core/initial_conditions/initialConditions.cc +++ b/src/core/initial_conditions/initialConditions.cc @@ -1,4 +1,6 @@ #include +#include +#include #include #include #include @@ -100,14 +102,15 @@ MatrixFreePDE::applyInitialConditions() } else { - pcout << "Error in vtk file type: Use either UNSTRUCTURED OR RECTILINEAR\n"; + conditionalOStreams::pout_base + << "Error in vtk file type: Use either UNSTRUCTURED OR RECTILINEAR\n"; abort(); } // new section ends ScalarField &id_field = body.find_scalar_field(userInputs.grain_structure_variable_name); - pcout << "Applying PField initial condition...\n"; + conditionalOStreams::pout_base << "Applying PField initial condition...\n"; VectorTools::interpolate(*dofHandlersSet[scalar_field_index], InitialConditionPField(0, id_field), @@ -124,11 +127,11 @@ MatrixFreePDE::applyInitialConditions() QGaussLobatto quadrature2(degree + 1); FloodFiller flood_filler(*FESet.at(scalar_field_index), quadrature2); - pcout << "Locating the grains...\n"; + conditionalOStreams::pout_base << "Locating the grains...\n"; std::vector> grain_sets; for (unsigned int id = min_id; id < max_id + 1; id++) { - pcout << "Locating grain " << id << "...\n"; + conditionalOStreams::pout_base << "Locating grain " << id << "...\n"; std::vector> grain_sets_single_id; @@ -151,7 +154,8 @@ MatrixFreePDE::applyInitialConditions() grain_sets_single_id.end()); } - pcout << "Generating simplified representations of the grains...\n"; + conditionalOStreams::pout_base + << "Generating simplified representations of the grains...\n"; for (unsigned int g = 0; g < grain_sets.size(); g++) { SimplifiedGrainRepresentation simplified_grain_representation( @@ -159,22 +163,24 @@ MatrixFreePDE::applyInitialConditions() if (dim == 2) { - pcout << "Grain: " << simplified_grain_representation.getGrainId() << " " - << simplified_grain_representation.getOrderParameterId() - << " Center: " << simplified_grain_representation.getCenter()(0) - << " " << simplified_grain_representation.getCenter()(1) - << " Radius: " << simplified_grain_representation.getRadius() - << std::endl; + conditionalOStreams::pout_base + << "Grain: " << simplified_grain_representation.getGrainId() << " " + << simplified_grain_representation.getOrderParameterId() + << " Center: " << simplified_grain_representation.getCenter()(0) << " " + << simplified_grain_representation.getCenter()(1) + << " Radius: " << simplified_grain_representation.getRadius() + << std::endl; } else { - pcout << "Grain: " << simplified_grain_representation.getGrainId() << " " - << simplified_grain_representation.getOrderParameterId() - << " Center: " << simplified_grain_representation.getCenter()(0) - << " " << simplified_grain_representation.getCenter()(1) << " " - << simplified_grain_representation.getCenter()(2) - << " Radius: " << simplified_grain_representation.getRadius() - << std::endl; + conditionalOStreams::pout_base + << "Grain: " << simplified_grain_representation.getGrainId() << " " + << simplified_grain_representation.getOrderParameterId() + << " Center: " << simplified_grain_representation.getCenter()(0) << " " + << simplified_grain_representation.getCenter()(1) << " " + << simplified_grain_representation.getCenter()(2) + << " Radius: " << simplified_grain_representation.getRadius() + << std::endl; } simplified_grain_representations.push_back(simplified_grain_representation); @@ -193,35 +199,38 @@ MatrixFreePDE::applyInitialConditions() } } - pcout << "Reassigning the grains to new order parameters...\n"; + conditionalOStreams::pout_base + << "Reassigning the grains to new order parameters...\n"; SimplifiedGrainManipulator simplified_grain_manipulator; simplified_grain_manipulator.reassignGrains(simplified_grain_representations, userInputs.buffer_between_grains, userInputs.variables_for_remapping); - pcout << "After reassignment: " << std::endl; + conditionalOStreams::pout_base << "After reassignment: " << std::endl; for (unsigned int g = 0; g < simplified_grain_representations.size(); g++) { if (dim == 2) { - pcout << "Grain: " << simplified_grain_representations.at(g).getGrainId() - << " " << simplified_grain_representations.at(g).getOrderParameterId() - << " Center: " - << simplified_grain_representations.at(g).getCenter()(0) << " " - << simplified_grain_representations.at(g).getCenter()(1) << std::endl; + conditionalOStreams::pout_base + << "Grain: " << simplified_grain_representations.at(g).getGrainId() << " " + << simplified_grain_representations.at(g).getOrderParameterId() + << " Center: " << simplified_grain_representations.at(g).getCenter()(0) + << " " << simplified_grain_representations.at(g).getCenter()(1) + << std::endl; } else { - pcout << "Grain: " << simplified_grain_representations.at(g).getGrainId() - << " " << simplified_grain_representations.at(g).getOrderParameterId() - << " Center: " - << simplified_grain_representations.at(g).getCenter()(0) << " " - << simplified_grain_representations.at(g).getCenter()(1) << " " - << simplified_grain_representations.at(g).getCenter()(2) << std::endl; + conditionalOStreams::pout_base + << "Grain: " << simplified_grain_representations.at(g).getGrainId() << " " + << simplified_grain_representations.at(g).getOrderParameterId() + << " Center: " << simplified_grain_representations.at(g).getCenter()(0) + << " " << simplified_grain_representations.at(g).getCenter()(1) << " " + << simplified_grain_representations.at(g).getCenter()(2) << std::endl; } } - pcout << "Placing the grains in their new order parameters...\n"; + conditionalOStreams::pout_base + << "Placing the grains in their new order parameters...\n"; OrderParameterRemapper order_parameter_remapper; order_parameter_remapper.remap_from_index_field( simplified_grain_representations, @@ -243,55 +252,55 @@ MatrixFreePDE::applyInitialConditions() double dt_for_smoothing = 0.25 * min_dx * min_dx / (1.0 * dim); op_list_index = 0; - for (unsigned int fieldIndex = 0; fieldIndex < fields.size(); fieldIndex++) - { - if (op_list_index < userInputs.variables_for_remapping.size()) - { - if (fieldIndex == userInputs.variables_for_remapping.at(op_list_index)) - { - for (unsigned int cycle = 0; - cycle < userInputs.num_grain_smoothing_cycles; - cycle++) - { - // Calculates the Laplace RHS and stores the information - // in residualSet - computeLaplaceRHS(fieldIndex); - if (fields[fieldIndex].type == SCALAR) - { - unsigned int invM_size = invMscalar.locally_owned_size(); - for (unsigned int dof = 0; - dof < solutionSet[fieldIndex]->locally_owned_size(); - ++dof) - { - solutionSet[fieldIndex]->local_element(dof) = - solutionSet[fieldIndex]->local_element(dof) - - invMscalar.local_element(dof % invM_size) * - residualSet[fieldIndex]->local_element(dof) * - dt_for_smoothing; - } - } - else if (fields[fieldIndex].type == VECTOR) - { - unsigned int invM_size = invMvector.locally_owned_size(); - for (unsigned int dof = 0; - dof < solutionSet[fieldIndex]->locally_owned_size(); - ++dof) - { - solutionSet[fieldIndex]->local_element(dof) = - solutionSet[fieldIndex]->local_element(dof) - - invMvector.local_element(dof % invM_size) * - residualSet[fieldIndex]->local_element(dof) * - dt_for_smoothing; - } - } - - solutionSet[fieldIndex]->update_ghost_values(); - } - - op_list_index++; - } - } - } + // for (unsigned int fieldIndex = 0; fieldIndex < fields.size(); fieldIndex++) + // { + // if (op_list_index < userInputs.variables_for_remapping.size()) + // { + // if (fieldIndex == userInputs.variables_for_remapping.at(op_list_index)) + // { + // for (unsigned int cycle = 0; + // cycle < userInputs.num_grain_smoothing_cycles; + // cycle++) + // { + // // Calculates the Laplace RHS and stores the information + // // in residualSet + // // computeLaplaceRHS(fieldIndex); + // // if (fields[fieldIndex].type == SCALAR) + // // { + // // unsigned int invM_size = invMscalar.locally_owned_size(); + // // for (unsigned int dof = 0; + // // dof < solutionSet[fieldIndex]->locally_owned_size(); + // // ++dof) + // // { + // // solutionSet[fieldIndex]->local_element(dof) = + // // solutionSet[fieldIndex]->local_element(dof) - + // // invMscalar.local_element(dof % invM_size) * + // // residualSet[fieldIndex]->local_element(dof) * + // // dt_for_smoothing; + // // } + // // } + // // else if (fields[fieldIndex].type == VECTOR) + // // { + // // unsigned int invM_size = invMvector.locally_owned_size(); + // // for (unsigned int dof = 0; + // // dof < solutionSet[fieldIndex]->locally_owned_size(); + // // ++dof) + // // { + // // solutionSet[fieldIndex]->local_element(dof) = + // // solutionSet[fieldIndex]->local_element(dof) - + // // invMvector.local_element(dof % invM_size) * + // // residualSet[fieldIndex]->local_element(dof) * + // // dt_for_smoothing; + // // } + // // } + + // solutionSet[fieldIndex]->update_ghost_values(); + // } + + // op_list_index++; + // } + // } + // } } // Create map of vtk files that are read in @@ -304,12 +313,12 @@ MatrixFreePDE::applyInitialConditions() } } // Read out unique filenames - if (Utilities::MPI::this_mpi_process(MPI_COMM_WORLD) == 0 && !file_field_map.empty()) + if (!file_field_map.empty()) { - std::cout << "Unique VTK input files: " << std::endl; + conditionalOStreams::pout_base << "Unique VTK input files: " << std::endl; for (const auto &pair : file_field_map) { - std::cout << pair.first << ", "; + conditionalOStreams::pout_base << pair.first << ", "; } } // Read in each vtk once and apply initial conditions @@ -352,7 +361,8 @@ MatrixFreePDE::applyInitialConditions() } else { - pcout << "Error in vtk file type: Use either UNSTRUCTURED OR RECTILINEAR\n"; + conditionalOStreams::pout_base + << "Error in vtk file type: Use either UNSTRUCTURED OR RECTILINEAR\n"; abort(); } // new section ends @@ -365,8 +375,9 @@ MatrixFreePDE::applyInitialConditions() if (var_attributes.at(index).var_type == SCALAR) { - pcout << "Applying PField initial condition for " - << userInputs.load_field_name[index] << "...\n"; + conditionalOStreams::pout_base << "Applying PField initial condition for " + << userInputs.load_field_name[index] + << "...\n"; VectorTools::interpolate(*dofHandlersSet[index], InitialConditionPField(index, field), *solutionSet[index]); @@ -397,7 +408,8 @@ MatrixFreePDE::applyInitialConditions() { if (userInputs.load_ICs[var_index] == false) { - pcout << "Applying non-PField initial condition...\n"; + conditionalOStreams::pout_base + << "Applying non-PField initial condition...\n"; if (variable.var_type == SCALAR) { @@ -416,26 +428,11 @@ MatrixFreePDE::applyInitialConditions() *solutionSet[var_index]); } } - pcout << "Application of initial conditions for field number " << var_index - << " complete \n"; + conditionalOStreams::pout_base + << "Application of initial conditions for field number " << var_index + << " complete \n"; } } } -template class MatrixFreePDE<2, 1>; -template class MatrixFreePDE<3, 1>; - -template class MatrixFreePDE<2, 2>; -template class MatrixFreePDE<3, 2>; - -template class MatrixFreePDE<3, 3>; -template class MatrixFreePDE<2, 3>; - -template class MatrixFreePDE<3, 4>; -template class MatrixFreePDE<2, 4>; - -template class MatrixFreePDE<3, 5>; -template class MatrixFreePDE<2, 5>; - -template class MatrixFreePDE<3, 6>; -template class MatrixFreePDE<2, 6>; \ No newline at end of file +INSTANTIATE_BI_TEMPLATE(MatrixFreePDE) \ No newline at end of file diff --git a/src/core/inputFileReader.cc b/src/core/inputFileReader.cc index 16129d6f..fe483fc4 100644 --- a/src/core/inputFileReader.cc +++ b/src/core/inputFileReader.cc @@ -1,6 +1,7 @@ #include #include +#include #include #include #include @@ -16,12 +17,9 @@ inputFileReader::inputFileReader(const std::string &input_file_name, model_constant_names = get_model_constant_names(); uint num_constants = model_constant_names.size(); - if (dealii::Utilities::MPI::this_mpi_process(MPI_COMM_WORLD) == 0) - { - std::cout << "Number of constants: " << num_constants << "\n"; - std::cout << "Number of post-processing variables: " << pp_attributes.size() - << "\n"; - } + conditionalOStreams::pout_base + << "Number of constants: " << num_constants << "\n" + << "Number of post-processing variables: " << pp_attributes.size() << "\n"; // Read in all of the parameters now declare_parameters(); diff --git a/src/core/invM.cc b/src/core/invM.cc index b98f0ae4..d89e00a2 100644 --- a/src/core/invM.cc +++ b/src/core/invM.cc @@ -1,6 +1,7 @@ -// computeInvM() method for MatrixFreePDE class #include +#include +#include #include // compute inverse of the diagonal mass matrix and store in vector invM @@ -22,25 +23,25 @@ MatrixFreePDE::computeInvM() unsigned int parabolicScalarFieldIndex = 0; unsigned int parabolicVectorFieldIndex = 0; - for (unsigned int fieldIndex = 0; fieldIndex < fields.size(); fieldIndex++) - { - if (fields[fieldIndex].pdetype == EXPLICIT_TIME_DEPENDENT || - fields[fieldIndex].pdetype == AUXILIARY) - { - if (fields[fieldIndex].type == SCALAR && !invMscalarFound) - { - parabolicScalarFieldIndex = fieldIndex; - invMscalarFound = true; - continue; - } - else if (fields[fieldIndex].type == VECTOR && !invMvectorFound) - { - parabolicVectorFieldIndex = fieldIndex; - invMvectorFound = true; - continue; - } - } - } + // for (unsigned int fieldIndex = 0; fieldIndex < fields.size(); fieldIndex++) + // { + // if (fields[fieldIndex].pdetype == EXPLICIT_TIME_DEPENDENT || + // fields[fieldIndex].pdetype == AUXILIARY) + // { + // if (fields[fieldIndex].type == SCALAR && !invMscalarFound) + // { + // parabolicScalarFieldIndex = fieldIndex; + // invMscalarFound = true; + // continue; + // } + // else if (fields[fieldIndex].type == VECTOR && !invMvectorFound) + // { + // parabolicVectorFieldIndex = fieldIndex; + // invMvectorFound = true; + // continue; + // } + // } + // } // Initialize invM and clear its values matrixFreeObject.initialize_dof_vector(invMscalar, parabolicScalarFieldIndex); @@ -53,46 +54,46 @@ MatrixFreePDE::computeInvM() // Compute mass matrix for the given type of quadrature. Selecting gauss // lobatto quadrature points which are suboptimal but give diagonal M - if (fields[parabolicScalarFieldIndex].type == SCALAR) - { - VectorizedArray one = make_vectorized_array(1.0); - FEEvaluation fe_eval(matrixFreeObject, parabolicScalarFieldIndex); - - const unsigned int n_q_points = fe_eval.n_q_points; - for (unsigned int cell = 0; cell < matrixFreeObject.n_cell_batches(); ++cell) - { - fe_eval.reinit(cell); - for (unsigned int q = 0; q < n_q_points; ++q) - { - fe_eval.submit_value(one, q); - } - fe_eval.integrate(invM_flags); - fe_eval.distribute_local_to_global(invMscalar); - } - } - if (fields[parabolicVectorFieldIndex].type == VECTOR) - { - dealii::Tensor<1, dim, dealii::VectorizedArray> oneV; - for (unsigned int i = 0; i < dim; i++) - { - oneV[i] = 1.0; - } - - FEEvaluation fe_eval(matrixFreeObject, - parabolicVectorFieldIndex); - - const unsigned int n_q_points = fe_eval.n_q_points; - for (unsigned int cell = 0; cell < matrixFreeObject.n_cell_batches(); ++cell) - { - fe_eval.reinit(cell); - for (unsigned int q = 0; q < n_q_points; ++q) - { - fe_eval.submit_value(oneV, q); - } - fe_eval.integrate(invM_flags); - fe_eval.distribute_local_to_global(invMvector); - } - } + // if (fields[parabolicScalarFieldIndex].type == SCALAR) + // { + // VectorizedArray one = make_vectorized_array(1.0); + // FEEvaluation fe_eval(matrixFreeObject, parabolicScalarFieldIndex); + + // const unsigned int n_q_points = fe_eval.n_q_points; + // for (unsigned int cell = 0; cell < matrixFreeObject.n_cell_batches(); ++cell) + // { + // fe_eval.reinit(cell); + // for (unsigned int q = 0; q < n_q_points; ++q) + // { + // fe_eval.submit_value(one, q); + // } + // fe_eval.integrate(invM_flags); + // fe_eval.distribute_local_to_global(invMscalar); + // } + // } + // if (fields[parabolicVectorFieldIndex].type == VECTOR) + // { + // dealii::Tensor<1, dim, dealii::VectorizedArray> oneV; + // for (unsigned int i = 0; i < dim; i++) + // { + // oneV[i] = 1.0; + // } + + // FEEvaluation fe_eval(matrixFreeObject, + // parabolicVectorFieldIndex); + + // const unsigned int n_q_points = fe_eval.n_q_points; + // for (unsigned int cell = 0; cell < matrixFreeObject.n_cell_batches(); ++cell) + // { + // fe_eval.reinit(cell); + // for (unsigned int q = 0; q < n_q_points; ++q) + // { + // fe_eval.submit_value(oneV, q); + // } + // fe_eval.integrate(invM_flags); + // fe_eval.distribute_local_to_global(invMvector); + // } + // } invMscalar.compress(VectorOperation::add); invMvector.compress(VectorOperation::add); @@ -125,8 +126,9 @@ MatrixFreePDE::computeInvM() invMscalar.local_element(k) = 0; } } - pcout << "computed scalar mass matrix (using FE space for field: " - << parabolicScalarFieldIndex << ")\n"; + conditionalOStreams::pout_base + << "computed scalar mass matrix (using FE space for field: " + << parabolicScalarFieldIndex << ")\n"; // Invert vector mass matrix diagonal elements for (unsigned int k = 0; k < invMvector.locally_owned_size(); ++k) @@ -140,24 +142,9 @@ MatrixFreePDE::computeInvM() invMvector.local_element(k) = 0; } } - pcout << "computed vector mass matrix (using FE space for field: " - << parabolicVectorFieldIndex << ")\n"; + conditionalOStreams::pout_base + << "computed vector mass matrix (using FE space for field: " + << parabolicVectorFieldIndex << ")\n"; } -template class MatrixFreePDE<2, 1>; -template class MatrixFreePDE<3, 1>; - -template class MatrixFreePDE<2, 2>; -template class MatrixFreePDE<3, 2>; - -template class MatrixFreePDE<3, 3>; -template class MatrixFreePDE<2, 3>; - -template class MatrixFreePDE<3, 4>; -template class MatrixFreePDE<2, 4>; - -template class MatrixFreePDE<3, 5>; -template class MatrixFreePDE<2, 5>; - -template class MatrixFreePDE<3, 6>; -template class MatrixFreePDE<2, 6>; \ No newline at end of file +INSTANTIATE_BI_TEMPLATE(MatrixFreePDE) \ No newline at end of file diff --git a/src/core/matrixFreePDE.cc b/src/core/matrixFreePDE.cc index c73e7cbe..67b836b8 100644 --- a/src/core/matrixFreePDE.cc +++ b/src/core/matrixFreePDE.cc @@ -1,12 +1,11 @@ -// constructor and destructor for matrixFreePDE class - +#include +#include #include // constructor template MatrixFreePDE::MatrixFreePDE(userInputParameters _userInputs) : Subscriptor() - , pcout(std::cout, Utilities::MPI::this_mpi_process(MPI_COMM_WORLD) == 0) , userInputs(_userInputs) , var_attributes(_userInputs.var_attributes) , pp_attributes(_userInputs.pp_attributes) @@ -21,11 +20,9 @@ MatrixFreePDE::MatrixFreePDE(userInputParameters _userInputs) , currentOutput(0) , currentCheckpoint(0) , current_grain_reassignment(0) - , computing_timer(pcout, TimerOutput::summary, TimerOutput::wall_times) , first_integrated_var_output_complete(false) , AMR(_userInputs, triangulation, - fields, solutionSet, soltransSet, FESet, @@ -73,20 +70,4 @@ MatrixFreePDE::~MatrixFreePDE() } } -template class MatrixFreePDE<2, 1>; -template class MatrixFreePDE<3, 1>; - -template class MatrixFreePDE<2, 2>; -template class MatrixFreePDE<3, 2>; - -template class MatrixFreePDE<3, 3>; -template class MatrixFreePDE<2, 3>; - -template class MatrixFreePDE<3, 4>; -template class MatrixFreePDE<2, 4>; - -template class MatrixFreePDE<3, 5>; -template class MatrixFreePDE<2, 5>; - -template class MatrixFreePDE<3, 6>; -template class MatrixFreePDE<2, 6>; \ No newline at end of file +INSTANTIATE_BI_TEMPLATE(MatrixFreePDE) \ No newline at end of file diff --git a/src/core/outputResults.cc b/src/core/outputResults.cc index 68d1d9ed..5f3d7e37 100644 --- a/src/core/outputResults.cc +++ b/src/core/outputResults.cc @@ -1,7 +1,10 @@ #include #include +#include +#include #include +#include #include #include @@ -9,30 +12,27 @@ template void MatrixFreePDE::outputResults() { - // log time - computing_timer.enter_subsection("matrixFreePDE: output"); - // create DataOut object DataOut data_out; // loop over fields - for (unsigned int fieldIndex = 0; fieldIndex < fields.size(); fieldIndex++) - { - // mark field as scalar/vector - std::vector dataType( - fields[fieldIndex].numComponents, - (fields[fieldIndex].type == SCALAR - ? DataComponentInterpretation::component_is_scalar - : DataComponentInterpretation::component_is_part_of_vector)); - // add field to data_out - std::vector solutionNames(fields[fieldIndex].numComponents, - fields[fieldIndex].name.c_str()); - data_out.add_data_vector(*dofHandlersSet[fieldIndex], - *solutionSet[fieldIndex], - solutionNames, - dataType); - } + // for (unsigned int fieldIndex = 0; fieldIndex < fields.size(); fieldIndex++) + // { + // // mark field as scalar/vector + // std::vector dataType( + // fields[fieldIndex].numComponents, + // (fields[fieldIndex].type == SCALAR + // ? DataComponentInterpretation::component_is_scalar + // : DataComponentInterpretation::component_is_part_of_vector)); + // // add field to data_out + // std::vector solutionNames(fields[fieldIndex].numComponents, + // fields[fieldIndex].name.c_str()); + // data_out.add_data_vector(*dofHandlersSet[fieldIndex], + // *solutionSet[fieldIndex], + // solutionNames, + // dataType); + // } // Test section for outputting postprocessed fields // Currently there are hacks in place, using the matrixFreeObject, invM, @@ -42,7 +42,7 @@ MatrixFreePDE::outputResults() if (userInputs.postProcessingRequired) { std::vector *> postProcessedSet; - computePostProcessedFields(postProcessedSet); + // computePostProcessedFields(postProcessedSet); unsigned int invM_size = invMscalar.locally_owned_size(); for (auto &field : postProcessedSet) { @@ -82,8 +82,9 @@ MatrixFreePDE::outputResults() { double integrated_field = NAN; computeIntegral(integrated_field, pp_index, postProcessedSet); - pcout << "Integrated value of " << pp_variable.name << ": " - << integrated_field << std::endl; + conditionalOStreams::pout_base << "Integrated value of " + << pp_variable.name << ": " + << integrated_field << std::endl; if (Utilities::MPI::this_mpi_process(MPI_COMM_WORLD) == 0) { output_file << "\t" << pp_variable.name << "\t" << integrated_field; @@ -103,7 +104,7 @@ MatrixFreePDE::outputResults() { // mark field as scalar/vector unsigned int components = 0; - if (userInputs.pp_varInfoList[fieldIndex].is_scalar) + if (userInputs.pp_varInfoList[fieldIndex].var_type == fieldType::SCALAR) { components = 1; std::vector @@ -115,7 +116,7 @@ MatrixFreePDE::outputResults() solutionNames, dataType); } - else + else if (userInputs.pp_varInfoList[fieldIndex].var_type == fieldType::VECTOR) { components = dim; std::vector @@ -190,7 +191,8 @@ MatrixFreePDE::outputResults() std::ofstream master_output(pvtuFileName); data_out.write_pvtu_record(master_output, filenames); - pcout << "Output written to:" << pvtuFileName << "\n\n"; + conditionalOStreams::pout_base << "Output written to:" << pvtuFileName + << "\n\n"; } } else @@ -201,7 +203,8 @@ MatrixFreePDE::outputResults() std::string svtuFileName = svtuFileNameStream.str(); data_out.write_vtu_in_parallel(svtuFileName, MPI_COMM_WORLD); - pcout << "Output written to:" << svtuFileName << "\n\n"; + conditionalOStreams::pout_base << "Output written to:" << svtuFileName + << "\n\n"; } } else if (userInputs.output_file_type == "vtk") @@ -209,7 +212,7 @@ MatrixFreePDE::outputResults() // Write the results to separate files for each process std::ofstream output(vtuFileName); data_out.write_vtk(output); - pcout << "Output written to:" << vtuFileName << "\n\n"; + conditionalOStreams::pout_base << "Output written to:" << vtuFileName << "\n\n"; } else { @@ -217,25 +220,6 @@ MatrixFreePDE::outputResults() "either \"vtu\" or \"vtk\"\n"; abort(); } - - // log time - computing_timer.leave_subsection("matrixFreePDE: output"); } -template class MatrixFreePDE<2, 1>; -template class MatrixFreePDE<3, 1>; - -template class MatrixFreePDE<2, 2>; -template class MatrixFreePDE<3, 2>; - -template class MatrixFreePDE<3, 3>; -template class MatrixFreePDE<2, 3>; - -template class MatrixFreePDE<3, 4>; -template class MatrixFreePDE<2, 4>; - -template class MatrixFreePDE<3, 5>; -template class MatrixFreePDE<2, 5>; - -template class MatrixFreePDE<3, 6>; -template class MatrixFreePDE<2, 6>; \ No newline at end of file +INSTANTIATE_BI_TEMPLATE(MatrixFreePDE) \ No newline at end of file diff --git a/src/core/postprocessing/computeIntegral.cc b/src/core/postprocessing/computeIntegral.cc index e8b2e5a0..cf2b6cba 100644 --- a/src/core/postprocessing/computeIntegral.cc +++ b/src/core/postprocessing/computeIntegral.cc @@ -1,4 +1,5 @@ +#include #include template @@ -39,27 +40,7 @@ MatrixFreePDE::computeIntegral( value = Utilities::MPI::sum(value, MPI_COMM_WORLD); - // if (Utilities::MPI::this_mpi_process(MPI_COMM_WORLD) == 0){ - // std::cout<<"Integrated field: "<; -template class MatrixFreePDE<3, 1>; - -template class MatrixFreePDE<2, 2>; -template class MatrixFreePDE<3, 2>; - -template class MatrixFreePDE<3, 3>; -template class MatrixFreePDE<2, 3>; - -template class MatrixFreePDE<3, 4>; -template class MatrixFreePDE<2, 4>; - -template class MatrixFreePDE<3, 5>; -template class MatrixFreePDE<2, 5>; - -template class MatrixFreePDE<3, 6>; -template class MatrixFreePDE<2, 6>; \ No newline at end of file +INSTANTIATE_BI_TEMPLATE(MatrixFreePDE) \ No newline at end of file diff --git a/src/core/postprocessing/postprocessor.cc b/src/core/postprocessing/postprocessor.cc deleted file mode 100644 index 2f477573..00000000 --- a/src/core/postprocessing/postprocessor.cc +++ /dev/null @@ -1,87 +0,0 @@ -#include - -template -void -MatrixFreePDE::computePostProcessedFields( - std::vector *> &postProcessedSet) -{ - // Initialize the postProcessedSet - for (unsigned int fieldIndex = 0; fieldIndex < pp_attributes.size(); fieldIndex++) - { - dealii::LinearAlgebra::distributed::Vector *U = nullptr; - U = new dealii::LinearAlgebra::distributed::Vector; - postProcessedSet.push_back(U); - matrixFreeObject.initialize_dof_vector(*U, 0); - } - - // call to integrate and assemble - matrixFreeObject.cell_loop(&MatrixFreePDE::getPostProcessedFields, - this, - postProcessedSet, - solutionSet, - true); -} - -template -void -MatrixFreePDE::getPostProcessedFields( - const dealii::MatrixFree &data, - std::vector *> &dst, - const std::vector *> &src, - const std::pair &cell_range) -{ - // initialize FEEvaulation objects - variableContainer> variable_list( - data, - userInputs.pp_baseVarInfoList); - variableContainer> - pp_variable_list(data, userInputs.pp_varInfoList, 0); - - // loop over cells - for (unsigned int cell = cell_range.first; cell < cell_range.second; ++cell) - { - // Initialize, read DOFs, and set evaulation flags for each variable - variable_list.reinit_and_eval(src, cell); - pp_variable_list.reinit(cell); - - unsigned int num_q_points = variable_list.get_num_q_points(); - - dealii::VectorizedArray local_element_volume = element_volume[cell]; - - // loop over quadrature points - for (unsigned int q = 0; q < num_q_points; ++q) - { - variable_list.q_point = q; - pp_variable_list.q_point = q; - - dealii::Point> q_point_loc = - variable_list.get_q_point_location(); - - // Calculate the residuals - postProcessedFields(variable_list, - pp_variable_list, - q_point_loc, - local_element_volume); - } - - pp_variable_list.integrate_and_distribute(dst); - } -} - -template class MatrixFreePDE<2, 1>; -template class MatrixFreePDE<3, 1>; - -template class MatrixFreePDE<2, 2>; -template class MatrixFreePDE<3, 2>; - -template class MatrixFreePDE<3, 3>; -template class MatrixFreePDE<2, 3>; - -template class MatrixFreePDE<3, 4>; -template class MatrixFreePDE<2, 4>; - -template class MatrixFreePDE<3, 5>; -template class MatrixFreePDE<2, 5>; - -template class MatrixFreePDE<3, 6>; -template class MatrixFreePDE<2, 6>; \ No newline at end of file diff --git a/src/core/refinement/AdaptiveRefinement.cc b/src/core/refinement/AdaptiveRefinement.cc index 47cf0b46..95b1a547 100644 --- a/src/core/refinement/AdaptiveRefinement.cc +++ b/src/core/refinement/AdaptiveRefinement.cc @@ -1,3 +1,4 @@ +#include #include using namespace dealii; @@ -6,7 +7,6 @@ template AdaptiveRefinement::AdaptiveRefinement( const userInputParameters &_userInputs, parallel::distributed::Triangulation &_triangulation, - std::vector> &_fields, std::vector *> &_solutionSet, std::vector::AdaptiveRefinement( std::vector *> &_constraintsOtherSet) : userInputs(_userInputs) , triangulation(_triangulation) - , fields(_fields) , solutionSet(_solutionSet) , soltransSet(_soltransSet) , FESet(_FESet) @@ -33,12 +32,12 @@ AdaptiveRefinement::do_adaptive_refinement(unsigned int currentIncr // Apply constraints for the initial condition so they are considered when remeshing if (currentIncrement != 0) { - for (unsigned int fieldIndex = 0; fieldIndex < fields.size(); fieldIndex++) - { - constraintsDirichletSet[fieldIndex]->distribute(*solutionSet[fieldIndex]); - constraintsOtherSet[fieldIndex]->distribute(*solutionSet[fieldIndex]); - solutionSet[fieldIndex]->update_ghost_values(); - } + // for (unsigned int fieldIndex = 0; fieldIndex < fields.size(); fieldIndex++) + // { + // constraintsDirichletSet[fieldIndex]->distribute(*solutionSet[fieldIndex]); + // constraintsOtherSet[fieldIndex]->distribute(*solutionSet[fieldIndex]); + // solutionSet[fieldIndex]->update_ghost_values(); + // } } adaptive_refinement_criterion(); @@ -158,31 +157,14 @@ AdaptiveRefinement::refine_grid() triangulation.prepare_coarsening_and_refinement(); // Transfer solution - for (unsigned int fieldIndex = 0; fieldIndex < fields.size(); fieldIndex++) - { - soltransSet[fieldIndex]->prepare_for_coarsening_and_refinement( - *solutionSet[fieldIndex]); - } + // for (unsigned int fieldIndex = 0; fieldIndex < fields.size(); fieldIndex++) + // { + // soltransSet[fieldIndex]->prepare_for_coarsening_and_refinement( + // *solutionSet[fieldIndex]); + // } // Execute refinement triangulation.execute_coarsening_and_refinement(); } -// Explicit instantiation -template class AdaptiveRefinement<2, 1>; -template class AdaptiveRefinement<3, 1>; - -template class AdaptiveRefinement<2, 2>; -template class AdaptiveRefinement<3, 2>; - -template class AdaptiveRefinement<2, 3>; -template class AdaptiveRefinement<3, 3>; - -template class AdaptiveRefinement<2, 4>; -template class AdaptiveRefinement<3, 4>; - -template class AdaptiveRefinement<2, 5>; -template class AdaptiveRefinement<3, 5>; - -template class AdaptiveRefinement<2, 6>; -template class AdaptiveRefinement<3, 6>; \ No newline at end of file +INSTANTIATE_BI_TEMPLATE(AdaptiveRefinement) \ No newline at end of file diff --git a/src/core/reinit.cc b/src/core/reinit.cc index 10085d93..f0e5f09e 100644 --- a/src/core/reinit.cc +++ b/src/core/reinit.cc @@ -1,5 +1,5 @@ -// reinit() method for MatrixFreePDE class - +#include +#include #include // populate with fields and setup matrix free system @@ -7,91 +7,89 @@ template void MatrixFreePDE::reinit() { - computing_timer.enter_subsection("matrixFreePDE: reinitialization"); - // setup system - pcout << "Reinitializing matrix free object\n"; + conditionalOStreams::pout_base << "Reinitializing matrix free object\n"; totalDOFs = 0; - for (const auto &field : fields) - { - currentFieldIndex = field.index; - - char buffer[100]; - - // create FESystem - FESystem *fe = nullptr; - fe = FESet.at(field.index); - - // distribute DOFs - DoFHandler *dof_handler = nullptr; - dof_handler = dofHandlersSet_nonconst.at(field.index); - - dof_handler->distribute_dofs(*fe); - totalDOFs += dof_handler->n_dofs(); - - // extract locally_relevant_dofs - IndexSet *locally_relevant_dofs = nullptr; - locally_relevant_dofs = locally_relevant_dofsSet_nonconst.at(field.index); - - locally_relevant_dofs->clear(); - DoFTools::extract_locally_relevant_dofs(*dof_handler, *locally_relevant_dofs); - - // create constraints - AffineConstraints *constraintsDirichlet = nullptr; - AffineConstraints *constraintsOther = nullptr; - - constraintsDirichlet = constraintsDirichletSet_nonconst.at(field.index); - constraintsOther = constraintsOtherSet_nonconst.at(field.index); - - constraintsDirichlet->clear(); - constraintsDirichlet->reinit(*locally_relevant_dofs); - constraintsOther->clear(); - constraintsOther->reinit(*locally_relevant_dofs); - - // Get hanging node constraints - DoFTools::make_hanging_node_constraints(*dof_handler, *constraintsOther); - - // Pin solution - if (userInputs.pinned_point.find(currentFieldIndex) != - userInputs.pinned_point.end()) - { - set_rigid_body_mode_constraints(constraintsOther, - dof_handler, - userInputs.pinned_point[currentFieldIndex]); - } - - // Get constraints for periodic BCs - setPeriodicityConstraints(constraintsOther, dof_handler); - - // Get constraints for Dirichlet BCs - applyDirichletBCs(); - - constraintsDirichlet->close(); - constraintsOther->close(); - - // Store Dirichlet BC DOF's - valuesDirichletSet[field.index]->clear(); - for (types::global_dof_index i = 0; i < dof_handler->n_dofs(); i++) - { - if (locally_relevant_dofs->is_element(i)) - { - if (constraintsDirichlet->is_constrained(i)) - { - (*valuesDirichletSet[field.index])[i] = - constraintsDirichlet->get_inhomogeneity(i); - } - } - } - - snprintf(buffer, - sizeof(buffer), - "field '%2s' DOF : %u (Constraint DOF : %u)\n", - field.name.c_str(), - dof_handler->n_dofs(), - constraintsDirichlet->n_constraints()); - pcout << buffer; - } - pcout << "total DOF : " << totalDOFs << "\n"; + // for (const auto &field : fields) + // { + // currentFieldIndex = field.index; + + // char buffer[100]; + + // // create FESystem + // FESystem *fe = nullptr; + // fe = FESet.at(field.index); + + // // distribute DOFs + // DoFHandler *dof_handler = nullptr; + // dof_handler = dofHandlersSet_nonconst.at(field.index); + + // dof_handler->distribute_dofs(*fe); + // totalDOFs += dof_handler->n_dofs(); + + // // extract locally_relevant_dofs + // IndexSet *locally_relevant_dofs = nullptr; + // locally_relevant_dofs = locally_relevant_dofsSet_nonconst.at(field.index); + + // locally_relevant_dofs->clear(); + // DoFTools::extract_locally_relevant_dofs(*dof_handler, *locally_relevant_dofs); + + // // create constraints + // AffineConstraints *constraintsDirichlet = nullptr; + // AffineConstraints *constraintsOther = nullptr; + + // constraintsDirichlet = constraintsDirichletSet_nonconst.at(field.index); + // constraintsOther = constraintsOtherSet_nonconst.at(field.index); + + // constraintsDirichlet->clear(); + // constraintsDirichlet->reinit(*locally_relevant_dofs); + // constraintsOther->clear(); + // constraintsOther->reinit(*locally_relevant_dofs); + + // // Get hanging node constraints + // DoFTools::make_hanging_node_constraints(*dof_handler, *constraintsOther); + + // // Pin solution + // if (userInputs.pinned_point.find(currentFieldIndex) != + // userInputs.pinned_point.end()) + // { + // set_rigid_body_mode_constraints(constraintsOther, + // dof_handler, + // userInputs.pinned_point[currentFieldIndex]); + // } + + // // Get constraints for periodic BCs + // setPeriodicityConstraints(constraintsOther, dof_handler); + + // // Get constraints for Dirichlet BCs + // applyDirichletBCs(); + + // constraintsDirichlet->close(); + // constraintsOther->close(); + + // // Store Dirichlet BC DOF's + // valuesDirichletSet[field.index]->clear(); + // for (types::global_dof_index i = 0; i < dof_handler->n_dofs(); i++) + // { + // if (locally_relevant_dofs->is_element(i)) + // { + // if (constraintsDirichlet->is_constrained(i)) + // { + // (*valuesDirichletSet[field.index])[i] = + // constraintsDirichlet->get_inhomogeneity(i); + // } + // } + // } + + // snprintf(buffer, + // sizeof(buffer), + // "field '%2s' DOF : %u (Constraint DOF : %u)\n", + // field.name.c_str(), + // dof_handler->n_dofs(), + // constraintsDirichlet->n_constraints()); + // conditionalOStreams::pout_base << buffer; + // } + conditionalOStreams::pout_base << "total DOF : " << totalDOFs << "\n"; // Setup the matrix free object typename MatrixFree::AdditionalData additional_data; @@ -117,41 +115,42 @@ MatrixFreePDE::reinit() bool dU_vector_init = false; // Setup solution vectors - pcout << "initializing parallel::distributed residual and solution vectors\n"; - for (unsigned int fieldIndex = 0; fieldIndex < fields.size(); fieldIndex++) - { - dealii::LinearAlgebra::distributed::Vector *U = nullptr; - - U = solutionSet.at(fieldIndex); - - matrixFreeObject.initialize_dof_vector(*U, fieldIndex); - *U = 0; - - // Initializing temporary dU vector required for implicit solves of the - // elliptic equation. - if (fields[fieldIndex].pdetype == TIME_INDEPENDENT || - fields[fieldIndex].pdetype == IMPLICIT_TIME_DEPENDENT || - (fields[fieldIndex].pdetype == AUXILIARY && - var_attributes.at(fieldIndex).is_nonlinear)) - { - if (fields[fieldIndex].type == SCALAR) - { - if (!dU_scalar_init) - { - matrixFreeObject.initialize_dof_vector(dU_scalar, fieldIndex); - dU_scalar_init = true; - } - } - else - { - if (!dU_vector_init) - { - matrixFreeObject.initialize_dof_vector(dU_vector, fieldIndex); - dU_vector_init = true; - } - } - } - } + conditionalOStreams::pout_base + << "initializing parallel::distributed residual and solution vectors\n"; + // for (unsigned int fieldIndex = 0; fieldIndex < fields.size(); fieldIndex++) + // { + // dealii::LinearAlgebra::distributed::Vector *U = nullptr; + + // U = solutionSet.at(fieldIndex); + + // matrixFreeObject.initialize_dof_vector(*U, fieldIndex); + // *U = 0; + + // // Initializing temporary dU vector required for implicit solves of the + // // elliptic equation. + // if (fields[fieldIndex].pdetype == TIME_INDEPENDENT || + // fields[fieldIndex].pdetype == IMPLICIT_TIME_DEPENDENT || + // (fields[fieldIndex].pdetype == AUXILIARY && + // var_attributes.at(fieldIndex).is_nonlinear)) + // { + // if (fields[fieldIndex].type == SCALAR) + // { + // if (!dU_scalar_init) + // { + // matrixFreeObject.initialize_dof_vector(dU_scalar, fieldIndex); + // dU_scalar_init = true; + // } + // } + // else + // { + // if (!dU_vector_init) + // { + // matrixFreeObject.initialize_dof_vector(dU_vector, fieldIndex); + // dU_vector_init = true; + // } + // } + // } + // } // Compute invM in PDE is a time-dependent BVP if (isTimeDependentBVP) @@ -160,27 +159,27 @@ MatrixFreePDE::reinit() } // Transfer solution from previous mesh - for (unsigned int fieldIndex = 0; fieldIndex < fields.size(); fieldIndex++) - { - // interpolate and clear used solution transfer sets - soltransSet[fieldIndex]->interpolate(*solutionSet[fieldIndex]); - delete soltransSet[fieldIndex]; - - // reset residual vector - dealii::LinearAlgebra::distributed::Vector *R = residualSet.at(fieldIndex); - matrixFreeObject.initialize_dof_vector(*R, fieldIndex); - *R = 0; - } + // for (unsigned int fieldIndex = 0; fieldIndex < fields.size(); fieldIndex++) + // { + // // interpolate and clear used solution transfer sets + // soltransSet[fieldIndex]->interpolate(*solutionSet[fieldIndex]); + // delete soltransSet[fieldIndex]; + + // // reset residual vector + // dealii::LinearAlgebra::distributed::Vector *R = residualSet.at(fieldIndex); + // matrixFreeObject.initialize_dof_vector(*R, fieldIndex); + // *R = 0; + // } // Create new solution transfer sets soltransSet.clear(); - for (unsigned int fieldIndex = 0; fieldIndex < fields.size(); fieldIndex++) - { - soltransSet.push_back( - new parallel::distributed:: - SolutionTransfer>( - *dofHandlersSet_nonconst[fieldIndex])); - } + // for (unsigned int fieldIndex = 0; fieldIndex < fields.size(); fieldIndex++) + // { + // soltransSet.push_back( + // new parallel::distributed:: + // SolutionTransfer>( + // *dofHandlersSet_nonconst[fieldIndex])); + // } // If remeshing at the zeroth time step, re-apply initial conditions so the // starting values are correct on the refined mesh @@ -191,33 +190,15 @@ MatrixFreePDE::reinit() // Ghost the solution vectors. Also apply the Dirichet BC's (if any) on the // solution vectors - for (unsigned int fieldIndex = 0; fieldIndex < fields.size(); fieldIndex++) - { - constraintsDirichletSet[fieldIndex]->distribute(*solutionSet[fieldIndex]); - constraintsOtherSet[fieldIndex]->distribute(*solutionSet[fieldIndex]); - solutionSet[fieldIndex]->update_ghost_values(); - } + // for (unsigned int fieldIndex = 0; fieldIndex < fields.size(); fieldIndex++) + // { + // constraintsDirichletSet[fieldIndex]->distribute(*solutionSet[fieldIndex]); + // constraintsOtherSet[fieldIndex]->distribute(*solutionSet[fieldIndex]); + // solutionSet[fieldIndex]->update_ghost_values(); + // } // Once the initial triangulation has been set, compute element volume compute_element_volume(); - - computing_timer.leave_subsection("matrixFreePDE: reinitialization"); } -template class MatrixFreePDE<2, 1>; -template class MatrixFreePDE<3, 1>; - -template class MatrixFreePDE<2, 2>; -template class MatrixFreePDE<3, 2>; - -template class MatrixFreePDE<3, 3>; -template class MatrixFreePDE<2, 3>; - -template class MatrixFreePDE<3, 4>; -template class MatrixFreePDE<2, 4>; - -template class MatrixFreePDE<3, 5>; -template class MatrixFreePDE<2, 5>; - -template class MatrixFreePDE<3, 6>; -template class MatrixFreePDE<2, 6>; \ No newline at end of file +INSTANTIATE_BI_TEMPLATE(MatrixFreePDE) \ No newline at end of file diff --git a/src/core/solvers/computeLHS.cc b/src/core/solvers/computeLHS.cc deleted file mode 100644 index b8b30ca1..00000000 --- a/src/core/solvers/computeLHS.cc +++ /dev/null @@ -1,110 +0,0 @@ -// vmult() and getLHS() method for MatrixFreePDE class - -#include - -// vmult operation for LHS -template - -void -MatrixFreePDE::vmult( - dealii::LinearAlgebra::distributed::Vector &dst, - const dealii::LinearAlgebra::distributed::Vector &src) const -{ - // log time - computing_timer.enter_subsection("matrixFreePDE: computeLHS"); - - // create temporary copy of src vector as src2, as vector src is marked const - // and cannot be changed - dealii::LinearAlgebra::distributed::Vector src2; - matrixFreeObject.initialize_dof_vector(src2, currentFieldIndex); - src2 = src; - - // call cell_loop - if (!generatingInitialGuess) - { - matrixFreeObject.cell_loop(&MatrixFreePDE::getLHS, - this, - dst, - src2, - true); - } - else - { - matrixFreeObject.cell_loop(&MatrixFreePDE::getLaplaceLHS, - this, - dst, - src2, - true); - } - - // Account for Dirichlet BC's (essentially copy dirichlet DOF values present in src to - // dst, although it is unclear why the constraints can't just be distributed here) - for (auto &it : *valuesDirichletSet[currentFieldIndex]) - { - if (dst.in_local_range(it.first)) - { - dst(it.first) = src(it.first); //*jacobianDiagonal(it->first); - } - } - - // end log - computing_timer.leave_subsection("matrixFreePDE: computeLHS"); -} - -template -void -MatrixFreePDE::getLHS( - const MatrixFree &data, - dealii::LinearAlgebra::distributed::Vector &dst, - const dealii::LinearAlgebra::distributed::Vector &src, - const std::pair &cell_range) const -{ - variableContainer> - variable_list(data, userInputs.varInfoListLHS, userInputs.varChangeInfoListLHS); - - // loop over cells - for (unsigned int cell = cell_range.first; cell < cell_range.second; ++cell) - { - // Initialize, read DOFs, and set evaulation flags for each variable - variable_list.reinit_and_eval(solutionSet, cell); - variable_list.reinit_and_eval_change_in_solution(src, cell, currentFieldIndex); - - unsigned int num_q_points = variable_list.get_num_q_points(); - - dealii::VectorizedArray local_element_volume = element_volume[cell]; - - // loop over quadrature points - for (unsigned int q = 0; q < num_q_points; ++q) - { - variable_list.q_point = q; - - dealii::Point> q_point_loc = - variable_list.get_q_point_location(); - - // Calculate the residuals - equationLHS(variable_list, q_point_loc, local_element_volume); - } - - // Integrate the residuals and distribute from local to global - variable_list.integrate_and_distribute_change_in_solution_LHS(dst, - currentFieldIndex); - } -} - -template class MatrixFreePDE<2, 1>; -template class MatrixFreePDE<3, 1>; - -template class MatrixFreePDE<2, 2>; -template class MatrixFreePDE<3, 2>; - -template class MatrixFreePDE<3, 3>; -template class MatrixFreePDE<2, 3>; - -template class MatrixFreePDE<3, 4>; -template class MatrixFreePDE<2, 4>; - -template class MatrixFreePDE<3, 5>; -template class MatrixFreePDE<2, 5>; - -template class MatrixFreePDE<3, 6>; -template class MatrixFreePDE<2, 6>; \ No newline at end of file diff --git a/src/core/solvers/computeRHS.cc b/src/core/solvers/computeRHS.cc deleted file mode 100644 index e69b071a..00000000 --- a/src/core/solvers/computeRHS.cc +++ /dev/null @@ -1,136 +0,0 @@ -// computeRHS() method for MatrixFreePDE class - -#include -#include - -// update RHS of each field -template -void -MatrixFreePDE::computeExplicitRHS() -{ - // log time - computing_timer.enter_subsection("matrixFreePDE: computeRHS"); - - // call to integrate and assemble while clearing residual vecotrs - matrixFreeObject.cell_loop(&MatrixFreePDE::getExplicitRHS, - this, - residualSet, - solutionSet, - true); - - // end log - computing_timer.leave_subsection("matrixFreePDE: computeRHS"); -} - -template -void -MatrixFreePDE::getExplicitRHS( - const MatrixFree &data, - std::vector *> &dst, - const std::vector *> &src, - const std::pair &cell_range) const -{ - variableContainer> variable_list( - data, - userInputs.varInfoListExplicitRHS); - - // loop over cells - for (unsigned int cell = cell_range.first; cell < cell_range.second; ++cell) - { - // Initialize, read DOFs, and set evaulation flags for each variable - variable_list.reinit_and_eval(src, cell); - - unsigned int num_q_points = variable_list.get_num_q_points(); - - dealii::VectorizedArray local_element_volume = element_volume[cell]; - - // loop over quadrature points - for (unsigned int q = 0; q < num_q_points; ++q) - { - variable_list.q_point = q; - - dealii::Point> q_point_loc = - variable_list.get_q_point_location(); - - // Calculate the residuals - explicitEquationRHS(variable_list, q_point_loc, local_element_volume); - } - - variable_list.integrate_and_distribute(dst); - } -} - -// update RHS of each field -template -void -MatrixFreePDE::computeNonexplicitRHS() -{ - // log time - computing_timer.enter_subsection("matrixFreePDE: computeRHS"); - - // call to integrate and assemble while clearing residual vecotrs - matrixFreeObject.cell_loop(&MatrixFreePDE::getNonexplicitRHS, - this, - residualSet, - solutionSet, - true); - - // end log - computing_timer.leave_subsection("matrixFreePDE: computeRHS"); -} - -template -void -MatrixFreePDE::getNonexplicitRHS( - const MatrixFree &data, - std::vector *> &dst, - const std::vector *> &src, - const std::pair &cell_range) const -{ - variableContainer> variable_list( - data, - userInputs.varInfoListNonexplicitRHS); - - // loop over cells - for (unsigned int cell = cell_range.first; cell < cell_range.second; ++cell) - { - // Initialize, read DOFs, and set evaulation flags for each variable - variable_list.reinit_and_eval(src, cell); - - unsigned int num_q_points = variable_list.get_num_q_points(); - - dealii::VectorizedArray local_element_volume = element_volume[cell]; - - // loop over quadrature points - for (unsigned int q = 0; q < num_q_points; ++q) - { - variable_list.q_point = q; - - dealii::Point> q_point_loc = - variable_list.get_q_point_location(); - - // Calculate the residuals - nonExplicitEquationRHS(variable_list, q_point_loc, local_element_volume); - } - - variable_list.integrate_and_distribute(dst); - } -} - -template class MatrixFreePDE<2, 1>; -template class MatrixFreePDE<3, 1>; - -template class MatrixFreePDE<2, 2>; -template class MatrixFreePDE<3, 2>; - -template class MatrixFreePDE<3, 3>; -template class MatrixFreePDE<2, 3>; - -template class MatrixFreePDE<3, 4>; -template class MatrixFreePDE<2, 4>; - -template class MatrixFreePDE<3, 5>; -template class MatrixFreePDE<2, 5>; - -template class MatrixFreePDE<3, 6>; -template class MatrixFreePDE<2, 6>; \ No newline at end of file diff --git a/src/core/solvers/setNonlinearEqInitialGuess.cc b/src/core/solvers/setNonlinearEqInitialGuess.cc deleted file mode 100644 index 17e52f16..00000000 --- a/src/core/solvers/setNonlinearEqInitialGuess.cc +++ /dev/null @@ -1,220 +0,0 @@ -// setNonlinearEqInitialGuess() method for MatrixFreePDE class -#include -#include - -#include -#include - -// solve each time increment -template -void -MatrixFreePDE::setNonlinearEqInitialGuess() -{ - // log time - computing_timer.enter_subsection("matrixFreePDE: setNonlinearEqInitialGuess"); - Timer time; - time.start(); - char buffer[200]; - - for (const auto &[fieldIndex, variable] : var_attributes) - { - if ((variable.eq_type == TIME_INDEPENDENT) && variable.is_nonlinear && - userInputs.nonlinear_solver_parameters.getLaplaceInitializationFlag(fieldIndex)) - { - currentFieldIndex = fieldIndex; // Used in computeLaplaceLHS() - - computeLaplaceRHS(fieldIndex); - - for (const auto &it : *valuesDirichletSet[fieldIndex]) - { - if (residualSet[fieldIndex]->in_local_range(it.first)) - { - (*residualSet[fieldIndex])(it.first) = 0.0; - } - } - - // solver controls - double tol_value = NAN; - if (userInputs.linear_solver_parameters.getToleranceType(fieldIndex) == - ABSOLUTE_RESIDUAL) - { - tol_value = - userInputs.linear_solver_parameters.getToleranceValue(fieldIndex); - } - else - { - tol_value = - userInputs.linear_solver_parameters.getToleranceValue(fieldIndex) * - residualSet[fieldIndex]->l2_norm(); - } - - SolverControl solver_control( - userInputs.linear_solver_parameters.getMaxIterations(fieldIndex), - tol_value); - - // Currently the only allowed solver is SolverCG, the SolverType input - // variable is a dummy - SolverCG> solver( - solver_control); - - // solve - try - { - if (fields[fieldIndex].type == SCALAR) - { - dU_scalar = 0.0; - solver.solve(*this, - dU_scalar, - *residualSet[fieldIndex], - IdentityMatrix(solutionSet[fieldIndex]->size())); - } - else - { - dU_vector = 0.0; - solver.solve(*this, - dU_vector, - *residualSet[fieldIndex], - IdentityMatrix(solutionSet[fieldIndex]->size())); - } - } - catch (...) - { - pcout << "\nWarning: implicit solver did not converge as per set " - "tolerances. consider increasing maxSolverIterations or " - "decreasing solverTolerance.\n"; - } - - if (fields[fieldIndex].type == SCALAR) - { - *solutionSet[fieldIndex] += dU_scalar; - } - else - { - *solutionSet[fieldIndex] += dU_vector; - } - - if (currentIncrement % userInputs.skip_print_steps == 0) - { - double dU_norm = NAN; - if (fields[fieldIndex].type == SCALAR) - { - dU_norm = dU_scalar.l2_norm(); - } - else - { - dU_norm = dU_vector.l2_norm(); - } - snprintf(buffer, - sizeof(buffer), - "field '%2s' [laplace solve for initial guess]: initial " - "residual:%12.6e, current residual:%12.6e, nsteps:%u, " - "tolerance criterion:%12.6e, solution: %12.6e, dU: %12.6e\n", - fields[fieldIndex].name.c_str(), - residualSet[fieldIndex]->l2_norm(), - solver_control.last_value(), - solver_control.last_step(), - solver_control.tolerance(), - solutionSet[fieldIndex]->l2_norm(), - dU_norm); - pcout << buffer << "\n"; - } - } - } - - if (currentIncrement % userInputs.skip_print_steps == 0) - { - pcout << "wall time: " << time.wall_time() << "s\n"; - } - // log time - computing_timer.leave_subsection("matrixFreePDE: setNonlinearEqInitialGuess"); -} - -template -void -MatrixFreePDE::computeLaplaceRHS(unsigned int fieldIndex) -{ - // log time - computing_timer.enter_subsection("matrixFreePDE: computeLaplaceRHS"); - - // call to integrate and assemble while clearing residual vecotrs - matrixFreeObject.cell_loop(&MatrixFreePDE::getLaplaceRHS, - this, - *residualSet[fieldIndex], - *solutionSet[fieldIndex], - true); - - // end log - computing_timer.leave_subsection("matrixFreePDE: computeLaplaceRHS"); -} - -template -void -MatrixFreePDE::getLaplaceRHS( - const MatrixFree &data, - dealii::LinearAlgebra::distributed::Vector &dst, - const dealii::LinearAlgebra::distributed::Vector &src, - const std::pair &cell_range) const -{ - FEEvaluation mat(data); - - dealii::EvaluationFlags::EvaluationFlags laplace_flags = - dealii::EvaluationFlags::gradients; - // loop over all "cells" - for (unsigned int cell = cell_range.first; cell < cell_range.second; ++cell) - { - mat.reinit(cell); - mat.read_dof_values(src); - mat.evaluate(laplace_flags); - for (unsigned int q = 0; q < mat.n_q_points; ++q) - { - mat.submit_gradient(mat.get_gradient(q), q); - } - mat.integrate(laplace_flags); - mat.distribute_local_to_global(dst); - } -} - -template -void -MatrixFreePDE::getLaplaceLHS( - const MatrixFree &data, - dealii::LinearAlgebra::distributed::Vector &dst, - const dealii::LinearAlgebra::distributed::Vector &src, - const std::pair &cell_range) const -{ - FEEvaluation mat(data); - - dealii::EvaluationFlags::EvaluationFlags laplace_flags = - dealii::EvaluationFlags::gradients; - // loop over all "cells" - for (unsigned int cell = cell_range.first; cell < cell_range.second; ++cell) - { - mat.reinit(cell); - mat.read_dof_values(src); - mat.evaluate(laplace_flags); - for (unsigned int q = 0; q < mat.n_q_points; ++q) - { - mat.submit_gradient(-mat.get_gradient(q), q); - } - mat.integrate(laplace_flags); - mat.distribute_local_to_global(dst); - } -} - -template class MatrixFreePDE<2, 1>; -template class MatrixFreePDE<3, 1>; - -template class MatrixFreePDE<2, 2>; -template class MatrixFreePDE<3, 2>; - -template class MatrixFreePDE<3, 3>; -template class MatrixFreePDE<2, 3>; - -template class MatrixFreePDE<3, 4>; -template class MatrixFreePDE<2, 4>; - -template class MatrixFreePDE<3, 5>; -template class MatrixFreePDE<2, 5>; - -template class MatrixFreePDE<3, 6>; -template class MatrixFreePDE<2, 6>; \ No newline at end of file diff --git a/src/core/solvers/solve.cc b/src/core/solvers/solve.cc index 48a083e3..e1748d22 100644 --- a/src/core/solvers/solve.cc +++ b/src/core/solvers/solve.cc @@ -1,5 +1,5 @@ -// solve() method for MatrixFreePDE class - +#include +#include #include // solve BVP @@ -7,9 +7,7 @@ template void MatrixFreePDE::solve() { - // log time - computing_timer.enter_subsection("matrixFreePDE: solve"); - pcout << "\nsolving...\n\n"; + conditionalOStreams::pout_base << "\nsolving...\n\n"; // time dependent BVP if (isTimeDependentBVP) @@ -24,22 +22,22 @@ MatrixFreePDE::solve() // For any nonlinear equation, set the initial guess as the solution to // Laplace's equations - generatingInitialGuess = true; - setNonlinearEqInitialGuess(); - generatingInitialGuess = false; + // generatingInitialGuess = true; + // setNonlinearEqInitialGuess(); + // generatingInitialGuess = false; // Do an initial solve to set the elliptic fields - solveIncrement(true); + // solveIncrement(true); // output initial conditions for time dependent BVP if (userInputs.outputTimeStepList[currentOutput] == currentIncrement) { - for (unsigned int fieldIndex = 0; fieldIndex < fields.size(); fieldIndex++) - { - constraintsDirichletSet[fieldIndex]->distribute(*solutionSet[fieldIndex]); - constraintsOtherSet[fieldIndex]->distribute(*solutionSet[fieldIndex]); - solutionSet[fieldIndex]->update_ghost_values(); - } + // for (unsigned int fieldIndex = 0; fieldIndex < fields.size(); fieldIndex++) + // { + // constraintsDirichletSet[fieldIndex]->distribute(*solutionSet[fieldIndex]); + // constraintsOtherSet[fieldIndex]->distribute(*solutionSet[fieldIndex]); + // solutionSet[fieldIndex]->update_ghost_values(); + // } outputResults(); currentOutput++; } @@ -67,9 +65,10 @@ MatrixFreePDE::solve() } // time stepping - pcout << "\nTime stepping parameters: timeStep: " << userInputs.dtValue - << " timeFinal: " << userInputs.finalTime - << " timeIncrements: " << userInputs.totalIncrements << "\n"; + conditionalOStreams::pout_base + << "\nTime stepping parameters: timeStep: " << userInputs.dtValue + << " timeFinal: " << userInputs.finalTime + << " timeIncrements: " << userInputs.totalIncrements << "\n"; // This is the main time-stepping loop for (; currentIncrement <= userInputs.totalIncrements; ++currentIncrement) @@ -78,20 +77,16 @@ MatrixFreePDE::solve() currentTime += userInputs.dtValue; if (currentIncrement % userInputs.skip_print_steps == 0) { - pcout << "\ntime increment:" << currentIncrement - << " time: " << currentTime << "\n"; + conditionalOStreams::pout_base << "\ntime increment:" << currentIncrement + << " time: " << currentTime << "\n"; } // check and perform adaptive mesh refinement if (userInputs.h_adaptivity == true && currentIncrement % userInputs.skip_remeshing_steps == 0) { - computing_timer.enter_subsection("matrixFreePDE: AMR"); - AMR.do_adaptive_refinement(currentIncrement); reinit(); - - computing_timer.leave_subsection("matrixFreePDE: AMR"); } // Update the list of nuclei (if relevant) @@ -106,24 +101,20 @@ MatrixFreePDE::solve() } // solve time increment - solveIncrement(false); + // solveIncrement(false); // Output results to file (on the proper increments) if (userInputs.outputTimeStepList[currentOutput] == currentIncrement) { - for (unsigned int fieldIndex = 0; fieldIndex < fields.size(); fieldIndex++) - { - constraintsDirichletSet[fieldIndex]->distribute( - *solutionSet[fieldIndex]); - constraintsOtherSet[fieldIndex]->distribute(*solutionSet[fieldIndex]); - solutionSet[fieldIndex]->update_ghost_values(); - } + // for (unsigned int fieldIndex = 0; fieldIndex < fields.size(); + // fieldIndex++) + // { + // constraintsDirichletSet[fieldIndex]->distribute( + // *solutionSet[fieldIndex]); + // constraintsOtherSet[fieldIndex]->distribute(*solutionSet[fieldIndex]); + // solutionSet[fieldIndex]->update_ghost_values(); + // } outputResults(); - if (userInputs.print_timing_with_output && - currentIncrement < userInputs.totalIncrements) - { - computing_timer.print_summary(); - } currentOutput++; } @@ -140,33 +131,14 @@ MatrixFreePDE::solve() // time independent BVP else { - generatingInitialGuess = false; + // generatingInitialGuess = false; // solve - solveIncrement(false); + // solveIncrement(false); // output results to file outputResults(); } - - // log time - computing_timer.leave_subsection("matrixFreePDE: solve"); } -template class MatrixFreePDE<2, 1>; -template class MatrixFreePDE<3, 1>; - -template class MatrixFreePDE<2, 2>; -template class MatrixFreePDE<3, 2>; - -template class MatrixFreePDE<3, 3>; -template class MatrixFreePDE<2, 3>; - -template class MatrixFreePDE<3, 4>; -template class MatrixFreePDE<2, 4>; - -template class MatrixFreePDE<3, 5>; -template class MatrixFreePDE<2, 5>; - -template class MatrixFreePDE<3, 6>; -template class MatrixFreePDE<2, 6>; \ No newline at end of file +INSTANTIATE_BI_TEMPLATE(MatrixFreePDE) \ No newline at end of file diff --git a/src/core/solvers/solveIncrement.cc b/src/core/solvers/solveIncrement.cc index 58d91399..6f314593 100644 --- a/src/core/solvers/solveIncrement.cc +++ b/src/core/solvers/solveIncrement.cc @@ -1,570 +1,531 @@ #include +#include + +#include "core/typeEnums.h" #include +#include #include #include +#include // solve each time increment template void -MatrixFreePDE::solveIncrement(bool skip_time_dependent) +MatrixFreePDE::solveIncrement(uint current_increment, + bool skip_time_dependent) { - // log time - computing_timer.enter_subsection("matrixFreePDE: solveIncrements"); - Timer time; - char buffer[200]; + char buffer[200]; - // Get the RHS of the explicit equations - if (hasExplicitEquation && !skip_time_dependent) + // Solve and update auxiliary fields that are not nonlinear (sequentially) + for (const auto &[var_index, variable] : var_attributes) { - computeExplicitRHS(); + if (variable.eq_type == AUXILIARY && !variable.is_nonlinear) + { + auxiliaryOnce(var_index); + } } - // solve for each field - for (unsigned int fieldIndex = 0; fieldIndex < fields.size(); fieldIndex++) + // Solve and update explicit fields (concurrently) + // if (hasExplicitEquation && !skip_time_dependent) + // computeExplicitRHS(); // HERE + for (const auto &[var_index, variable] : var_attributes) { - currentFieldIndex = fieldIndex; // Used in computeLHS() - - // Parabolic (first order derivatives in time) fields - if (fields[fieldIndex].pdetype == EXPLICIT_TIME_DEPENDENT && !skip_time_dependent) + if (variable.eq_type == EXPLICIT_TIME_DEPENDENT) { - updateExplicitSolution(fieldIndex); - + updateExplicitSolution(var_index); // Apply Boundary conditions - applyBCs(fieldIndex); - - // Print update to screen and confirm that solution isn't nan - if (currentIncrement % userInputs.skip_print_steps == 0) - { - double solution_L2_norm = solutionSet[fieldIndex]->l2_norm(); - - snprintf(buffer, - sizeof(buffer), - "field '%2s' [explicit solve]: current solution: " - "%12.6e, current residual:%12.6e\n", - fields[fieldIndex].name.c_str(), - solution_L2_norm, - residualSet[fieldIndex]->l2_norm()); - pcout << buffer; - - if (!numbers::is_finite(solution_L2_norm)) - { - snprintf(buffer, - sizeof(buffer), - "ERROR: field '%s' solution is NAN. exiting.\n\n", - fields[fieldIndex].name.c_str()); - pcout << buffer; - exit(-1); - } - } + applyBCs(var_index); + print_explicit_update(var_index, current_increment); } } - // Now, update the non-explicit variables - // For the time being, this is just the elliptic equations, but implicit - // parabolic and auxilary equations should also be here - if (hasNonExplicitEquation) + // Solve and update linear-equation fields (sequentially) + for (const auto &[var_index, variable] : var_attributes) { - bool nonlinear_iteration_converged = false; - unsigned int nonlinear_iteration_index = 0; - - while (!nonlinear_iteration_converged) + if (variable.eq_type != EXPLICIT_TIME_DEPENDENT && !variable.is_nonlinear) { - nonlinear_iteration_converged = true; + SolverControl solver_control; + // compute_nonexplicit_RHS(var_index); // HERE + linearSolveOnce(var_index, solver_control); // lhs + updateLinearSolution(var_index, solver_control); + print_linear_update(var_index, solver_control); - // Update residualSet for the non-explicitly updated variables - computeNonexplicitRHS(); + // Apply Boundary conditions + applyBCs(var_index); + // print_update(var_index, current_increment); + } + } - for (const auto &[fieldIndex, variable] : var_attributes) + // Solve and update nonlinear-equation fields (concurrently) + bool nonlinear_converged = false; + uint nonlinear_iteration_index = 0; + uint max_nonlinear_iterations = + userInputs.nonlinear_solver_parameters.getMaxIterations(); + SolverControl solver_control; + while (!nonlinear_converged && nonlinear_iteration_index < max_nonlinear_iterations) + { + nonlinear_converged = true; + for (const auto &[var_index, variable] : var_attributes) + { + if (!variable.is_nonlinear) { - currentFieldIndex = fieldIndex; // Used in computeLHS() - - if ((fields[fieldIndex].pdetype == IMPLICIT_TIME_DEPENDENT && - !skip_time_dependent) || - fields[fieldIndex].pdetype == TIME_INDEPENDENT) - { - if (currentIncrement % userInputs.skip_print_steps == 0 && - variable.is_nonlinear) - { - snprintf(buffer, - sizeof(buffer), - "field '%2s' [nonlinear solve]: current " - "solution: %12.6e, current residual:%12.6e\n", - fields[fieldIndex].name.c_str(), - solutionSet[fieldIndex]->l2_norm(), - residualSet[fieldIndex]->l2_norm()); - pcout << buffer; - } - - nonlinear_iteration_converged = - updateImplicitSolution(fieldIndex, nonlinear_iteration_index); - - // Apply Boundary conditions - applyBCs(fieldIndex); - } - else if (fields[fieldIndex].pdetype == AUXILIARY) + continue; + } + if (variable.eq_type == TIME_INDEPENDENT) + { + print_nonlinear_update(var_index, current_increment); + if (!nonlinearIncrement(var_index, solver_control)) { - if (variable.is_nonlinear || nonlinear_iteration_index == 0) - { - // If the equation for this field is nonlinear, save the old - // solution - if (variable.is_nonlinear) - { - if (fields[fieldIndex].type == SCALAR) - { - dU_scalar = *solutionSet[fieldIndex]; - } - else - { - dU_vector = *solutionSet[fieldIndex]; - } - } - - updateExplicitSolution(fieldIndex); - - // Apply Boundary conditions - applyBCs(fieldIndex); - - // Print update to screen - if (currentIncrement % userInputs.skip_print_steps == 0) - { - snprintf(buffer, - sizeof(buffer), - "field '%2s' [auxiliary solve]: current solution: " - "%12.6e, current residual:%12.6e\n", - fields[fieldIndex].name.c_str(), - solutionSet[fieldIndex]->l2_norm(), - residualSet[fieldIndex]->l2_norm()); - pcout << buffer; - } - - // Check to see if this individual variable has converged - if (variable.is_nonlinear) - { - if (userInputs.nonlinear_solver_parameters.getToleranceType( - fieldIndex) == ABSOLUTE_SOLUTION_CHANGE) - { - double diff = NAN; - - if (fields[fieldIndex].type == SCALAR) - { - dU_scalar -= *solutionSet[fieldIndex]; - diff = dU_scalar.l2_norm(); - } - else - { - dU_vector -= *solutionSet[fieldIndex]; - diff = dU_vector.l2_norm(); - } - if (currentIncrement % userInputs.skip_print_steps == 0) - { - snprintf(buffer, - sizeof(buffer), - " field '%2s' [nonlinear solve] current " - "increment: %u, nonlinear " - "iteration: " - "%u, dU: %12.6e\n", - fields[fieldIndex].name.c_str(), - currentIncrement, - nonlinear_iteration_index, - diff); - pcout << buffer; - } - - if (diff > userInputs.nonlinear_solver_parameters - .getToleranceValue(fieldIndex) && - nonlinear_iteration_index < - userInputs.nonlinear_solver_parameters - .getMaxIterations()) - { - nonlinear_iteration_converged = false; - } - } - else - { - AssertThrow( - false, - FeatureNotImplemented( - "Nonlinear solver tolerances besides ABSOLUTE_CHANGE")); - } - } - } + nonlinear_converged = false; } - - // check if solution is nan - if (!numbers::is_finite(solutionSet[fieldIndex]->l2_norm())) + // Apply Boundary conditions (could be placed in nonlinearIncrement) + applyBCs(var_index); + } + else if (variable.eq_type == AUXILIARY) + { + if (!auxiliaryIncrement(var_index, solver_control)) { - snprintf(buffer, - sizeof(buffer), - "ERROR: field '%s' solution is NAN. exiting.\n\n", - fields[fieldIndex].name.c_str()); - pcout << buffer; - exit(-1); + nonlinear_converged = false; } + // Apply Boundary conditions (could be placed in auxiliaryIncrement) + applyBCs(var_index); + } + // check if solution is nan + if (!numbers::is_finite(solutionSet[var_index]->l2_norm())) + { + // snprintf(buffer, + // sizeof(buffer), + // "ERROR: field '%s' solution is NAN. exiting.\n\n", + // fields[var_index].name.c_str()); + conditionalOStreams::pout_base << buffer; + exit(-1); } - - nonlinear_iteration_index++; } + nonlinear_iteration_index++; } - - if (currentIncrement % userInputs.skip_print_steps == 0) - { - pcout << "wall time: " << time.wall_time() << "s\n"; - } - // log time - computing_timer.leave_subsection("matrixFreePDE: solveIncrements"); -} - -// Application of boundary conditions -template -void -MatrixFreePDE::applyBCs(unsigned int fieldIndex) -{ - // Add Neumann BCs - if (fields[fieldIndex].hasNeumannBCs) - { - // Currently commented out because it isn't working yet - // applyNeumannBCs(); - } - - // Set the Dirichelet values (hanging node constraints don't need to be distributed - // every time step, only at output) - if (fields[fieldIndex].hasDirichletBCs) + if (!nonlinear_converged) { - // Apply non-uniform Dirlichlet_BCs to the current field - if (fields[fieldIndex].hasnonuniformDirichletBCs) - { - DoFHandler *dof_handler = nullptr; - dof_handler = dofHandlersSet_nonconst.at(currentFieldIndex); - IndexSet *locally_relevant_dofs = nullptr; - locally_relevant_dofs = locally_relevant_dofsSet_nonconst.at(currentFieldIndex); - locally_relevant_dofs->clear(); - DoFTools::extract_locally_relevant_dofs(*dof_handler, *locally_relevant_dofs); - AffineConstraints *constraintsDirichlet = nullptr; - constraintsDirichlet = constraintsDirichletSet_nonconst.at(currentFieldIndex); - constraintsDirichlet->clear(); - constraintsDirichlet->reinit(*locally_relevant_dofs); - applyDirichletBCs(); - constraintsDirichlet->close(); - } - // Distribute for Uniform or Non-Uniform Dirichlet BCs - constraintsDirichletSet[fieldIndex]->distribute(*solutionSet[fieldIndex]); + conditionalOStreams::pout_base << "\nWarning: nonlinear solver did not converge as " + "per set tolerances. consider increasing the " + "maximum number of iterations or decreasing the " + "solver tolerance.\n"; } - solutionSet[fieldIndex]->update_ghost_values(); } -// Explicit time step for matrixfree solve template void -MatrixFreePDE::updateExplicitSolution(unsigned int fieldIndex) +MatrixFreePDE::linearSolveOnce(unsigned int var_index, + SolverControl &solver_control) { - // Explicit-time step each DOF - // Takes advantage of knowledge that the length of solutionSet and residualSet - // is an integer multiple of the length of invM for vector variables - if (fields[fieldIndex].type == SCALAR) - { - unsigned int invM_size = invMscalar.locally_owned_size(); - for (unsigned int dof = 0; dof < solutionSet[fieldIndex]->locally_owned_size(); - ++dof) - { - solutionSet[fieldIndex]->local_element(dof) = - invMscalar.local_element(dof % invM_size) * - residualSet[fieldIndex]->local_element(dof); - } - } - else if (fields[fieldIndex].type == VECTOR) - { - unsigned int invM_size = invMvector.locally_owned_size(); - for (unsigned int dof = 0; dof < solutionSet[fieldIndex]->locally_owned_size(); - ++dof) - { - solutionSet[fieldIndex]->local_element(dof) = - invMvector.local_element(dof % invM_size) * - residualSet[fieldIndex]->local_element(dof); - } - } -} - -template -bool -MatrixFreePDE::updateImplicitSolution(unsigned int fieldIndex, - unsigned int nonlinear_iteration_index) -{ - char buffer[200]; - - // Assume convergence criterion is met, unless otherwise proven later on. - bool nonlinear_iteration_converged = true; - // Apply Dirichlet BC's. This clears the residual where we want to apply Dirichlet BCs, // otherwise the solver sees a positive residual - constraintsDirichletSet[fieldIndex]->set_zero(*residualSet[fieldIndex]); + constraintsDirichletSet[var_index]->set_zero(*residualSet[var_index]); // Grab solver controls double tol_value = NAN; - if (userInputs.linear_solver_parameters.getToleranceType(fieldIndex) == + if (userInputs.linear_solver_parameters.getToleranceType(var_index) == ABSOLUTE_RESIDUAL) { - tol_value = userInputs.linear_solver_parameters.getToleranceValue(fieldIndex); + tol_value = userInputs.linear_solver_parameters.getToleranceValue(var_index); } else { - tol_value = userInputs.linear_solver_parameters.getToleranceValue(fieldIndex) * - residualSet[fieldIndex]->l2_norm(); + tol_value = userInputs.linear_solver_parameters.getToleranceValue(var_index) * + residualSet[var_index]->l2_norm(); } - - SolverControl solver_control(userInputs.linear_solver_parameters.getMaxIterations( - fieldIndex), - tol_value); - - // Currently the only allowed solver is SolverCG, the - // SolverType input variable is a dummy + solver_control = + SolverControl(userInputs.linear_solver_parameters.getMaxIterations(var_index), + tol_value); SolverCG> solver(solver_control); - // Solve + // Solve the linear system try { - if (fields[fieldIndex].type == SCALAR) - { - dU_scalar = 0.0; - solver.solve(*this, - dU_scalar, - *residualSet[fieldIndex], - IdentityMatrix(solutionSet[fieldIndex]->size())); - } - else - { - dU_vector = 0.0; - solver.solve(*this, - dU_vector, - *residualSet[fieldIndex], - IdentityMatrix(solutionSet[fieldIndex]->size())); - } + // if (fields[var_index].type == SCALAR) + // { + // dU_scalar = 0.0; + // solver.solve(*this, + // dU_scalar, + // *residualSet[var_index], + // IdentityMatrix(solutionSet[var_index]->size())); + // } + // else + // { + // dU_vector = 0.0; + // solver.solve(*this, + // dU_vector, + // *residualSet[var_index], + // IdentityMatrix(solutionSet[var_index]->size())); + // } } catch (...) { - pcout << "\nWarning: linear solver did not converge as " - "per set tolerances. consider increasing the " - "maximum number of iterations or decreasing the " - "solver tolerance.\n"; + conditionalOStreams::pout_base << "\nWarning: linear solver did not converge as " + "per set tolerances. consider increasing the " + "maximum number of iterations or decreasing the " + "solver tolerance.\n"; } +} - if (var_attributes.at(fieldIndex).is_nonlinear) - { - // Now that we have the calculated change in the solution, - // we need to select a damping coefficient - double damping_coefficient = NAN; - - if (userInputs.nonlinear_solver_parameters.getBacktrackDampingFlag(fieldIndex)) - { - dealii::LinearAlgebra::distributed::Vector solutionSet_old = - *solutionSet[fieldIndex]; - double residual_old = residualSet[fieldIndex]->l2_norm(); - - damping_coefficient = 1.0; - bool damping_coefficient_found = false; - while (!damping_coefficient_found) - { - if (fields[fieldIndex].type == SCALAR) - { - solutionSet[fieldIndex]->sadd(1.0, damping_coefficient, dU_scalar); - } - else - { - solutionSet[fieldIndex]->sadd(1.0, damping_coefficient, dU_vector); - } +template +void +MatrixFreePDE::auxiliaryOnce(unsigned int var_index, uint current_increment) +{ + // computeNonexplicitRHS(var_index); // HERE + updateExplicitSolution(var_index); - computeNonexplicitRHS(); + // Apply Boundary conditions + applyBCs(var_index); - for (const auto &it : *valuesDirichletSet[fieldIndex]) - { - if (residualSet[fieldIndex]->in_local_range(it.first)) - { - (*residualSet[fieldIndex])(it.first) = 0.0; - } - } + // Print update to screen + if (current_increment % userInputs.skip_print_steps == 0) + { + // char buffer[200]; + // snprintf(buffer, + // sizeof(buffer), + // "field '%2s' [auxiliary solve]: current solution: " + // "%12.6e, current residual:%12.6e\n", + // fields[var_index].name.c_str(), + // solutionSet[var_index]->l2_norm(), + // residualSet[var_index]->l2_norm()); + // conditionalOStreams::pout_base << buffer; + } +} - double residual_new = residualSet[fieldIndex]->l2_norm(); +// Nonlinear increment for time-independent fields +template +bool +MatrixFreePDE::nonlinearIncrement(unsigned int var_index, + SolverControl &solver_control, + uint current_increment) +{ + linearSolveOnce(var_index, solver_control); // lhs - if (currentIncrement % userInputs.skip_print_steps == 0) - { - pcout << " Old residual: " << residual_old - << " Damping Coeff: " << damping_coefficient - << " New Residual: " << residual_new << "\n"; - } + uint nonlinear_iteration_index = 0; + double damping_coefficient = NAN; + if (userInputs.nonlinear_solver_parameters.getBacktrackDampingFlag(var_index)) + { + dealii::LinearAlgebra::distributed::Vector solutionSet_old = + *solutionSet[var_index]; + double residual_old = residualSet[var_index]->l2_norm(); - // An improved approach would use the - // Armijo–Goldstein condition to ensure a - // sufficent decrease in the residual. This way is - // just scales the residual. - if ((residual_new < - (residual_old * userInputs.nonlinear_solver_parameters - .getBacktrackResidualDecreaseCoeff(fieldIndex))) || - damping_coefficient < 1.0e-4) - { - damping_coefficient_found = true; - } - else + damping_coefficient = 1.0; + bool damping_coefficient_found = false; + while (!damping_coefficient_found) + { + // if (fields[var_index].type == SCALAR) + // { + // solutionSet[var_index]->sadd(1.0, damping_coefficient, dU_scalar); + // } + // else + // { + // solutionSet[var_index]->sadd(1.0, damping_coefficient, dU_vector); + // } + + // computeNonexplicitRHS(); // HERE + + for (const auto &it : *valuesDirichletSet[var_index]) + { + if (residualSet[var_index]->in_local_range(it.first)) { - damping_coefficient *= - userInputs.nonlinear_solver_parameters.getBacktrackStepModifier( - fieldIndex); - *solutionSet[fieldIndex] = solutionSet_old; + (*residualSet[var_index])(it.first) = 0.0; } } - } - else - { - damping_coefficient = - userInputs.nonlinear_solver_parameters.getDefaultDampingCoefficient( - fieldIndex); - if (fields[fieldIndex].type == SCALAR) - { - solutionSet[fieldIndex]->sadd(1.0, damping_coefficient, dU_scalar); - } - else - { - solutionSet[fieldIndex]->sadd(1.0, damping_coefficient, dU_vector); - } - } + double residual_new = residualSet[var_index]->l2_norm(); - if (currentIncrement % userInputs.skip_print_steps == 0) - { - double dU_norm = NAN; - if (fields[fieldIndex].type == SCALAR) - { - dU_norm = dU_scalar.l2_norm(); - } - else + if (current_increment % userInputs.skip_print_steps == 0) { - dU_norm = dU_vector.l2_norm(); + conditionalOStreams::pout_base << "\tOld residual: " << residual_old + << " Damping Coeff: " << damping_coefficient + << " New Residual: " << residual_new << "\n"; } - snprintf(buffer, - sizeof(buffer), - "field '%2s' [linear solve]: initial " - "residual:%12.6e, current residual:%12.6e, " - "nsteps:%u, tolerance criterion:%12.6e, " - "solution: %12.6e, dU: %12.6e\n", - fields[fieldIndex].name.c_str(), - residualSet[fieldIndex]->l2_norm(), - solver_control.last_value(), - solver_control.last_step(), - solver_control.tolerance(), - solutionSet[fieldIndex]->l2_norm(), - dU_norm); - pcout << buffer; - } - // Check to see if this individual variable has converged - if (userInputs.nonlinear_solver_parameters.getToleranceType(fieldIndex) == - ABSOLUTE_SOLUTION_CHANGE) - { - double diff = NAN; - - if (fields[fieldIndex].type == SCALAR) + // An improved approach would use the + // Armijo–Goldstein condition to ensure a + // sufficent decrease in the residual. This way is + // just scales the residual. + if ((residual_new < + (residual_old * + userInputs.nonlinear_solver_parameters.getBacktrackResidualDecreaseCoeff( + var_index))) || + damping_coefficient < 1.0e-4) { - diff = dU_scalar.l2_norm(); + damping_coefficient_found = true; } else { - diff = dU_vector.l2_norm(); - } - if (currentIncrement % userInputs.skip_print_steps == 0) - { - snprintf(buffer, - sizeof(buffer), - " field '%2s' [nonlinear solve] current increment: %u, nonlinear " - "iteration: " - "%u, dU: %12.6e\n", - fields[fieldIndex].name.c_str(), - currentIncrement, - nonlinear_iteration_index, - diff); - pcout << buffer; - } - - if (diff > - userInputs.nonlinear_solver_parameters.getToleranceValue(fieldIndex) && - nonlinear_iteration_index < - userInputs.nonlinear_solver_parameters.getMaxIterations()) - { - nonlinear_iteration_converged = false; - } - else if (diff > - userInputs.nonlinear_solver_parameters.getToleranceValue(fieldIndex)) - { - pcout << "\nWarning: nonlinear solver did not converge as " - "per set tolerances. consider increasing the " - "maximum number of iterations or decreasing the " - "solver tolerance.\n"; + damping_coefficient *= + userInputs.nonlinear_solver_parameters.getBacktrackStepModifier( + var_index); + *solutionSet[var_index] = solutionSet_old; } } - else - { - AssertThrow(false, - FeatureNotImplemented( - "Nonlinear solver tolerances besides ABSOLUTE_CHANGE")); - } } else { - if (nonlinear_iteration_index == 0) - { - if (fields[fieldIndex].type == SCALAR) - { - *solutionSet[fieldIndex] += dU_scalar; - } - else - { - *solutionSet[fieldIndex] += dU_vector; - } - - if (currentIncrement % userInputs.skip_print_steps == 0) - { - double dU_norm = NAN; - if (fields[fieldIndex].type == SCALAR) - { - dU_norm = dU_scalar.l2_norm(); - } - else - { - dU_norm = dU_vector.l2_norm(); - } - snprintf(buffer, - sizeof(buffer), - "field '%2s' [linear solve]: initial " - "residual:%12.6e, current residual:%12.6e, " - "nsteps:%u, tolerance criterion:%12.6e, " - "solution: %12.6e, dU: %12.6e\n", - fields[fieldIndex].name.c_str(), - residualSet[fieldIndex]->l2_norm(), - solver_control.last_value(), - solver_control.last_step(), - solver_control.tolerance(), - solutionSet[fieldIndex]->l2_norm(), - dU_norm); - pcout << buffer; - } - } + damping_coefficient = + userInputs.nonlinear_solver_parameters.getDefaultDampingCoefficient(var_index); + + // if (fields[var_index].type == SCALAR) + // { + // solutionSet[var_index]->sadd(1.0, damping_coefficient, dU_scalar); + // } + // else + // { + // solutionSet[var_index]->sadd(1.0, damping_coefficient, dU_vector); + // } } - return nonlinear_iteration_converged; + // Print linear???? + print_linear_update(var_index, solver_control); + + // Check to see if this individual variable has converged + bool nonlinear_converged = false; + AssertThrow(userInputs.nonlinear_solver_parameters.getToleranceType(var_index) == + ABSOLUTE_SOLUTION_CHANGE, + FeatureNotImplemented( + "Nonlinear solver tolerances besides ABSOLUTE_CHANGE")); + double diff = NAN; + // if (fields[var_index].type == SCALAR) + // { + // diff = dU_scalar.l2_norm(); + // } + // else + // { + // diff = dU_vector.l2_norm(); + // } + print_nonlinear_status(var_index, current_increment, nonlinear_iteration_index, diff); + return !(diff > userInputs.nonlinear_solver_parameters.getToleranceValue(var_index)); +} + +// Nonlinear increment for auxiliary fields +template +bool +MatrixFreePDE::auxiliaryIncrement(unsigned int var_index, + uint current_increment) +{ + // Save the old solution + // if (fields[var_index].type == SCALAR) + // { + // dU_scalar = *solutionSet[var_index]; + // } + // else + // { + // dU_vector = *solutionSet[var_index]; + // } + + auxiliaryOnce(var_index); + + // Check to see if this individual variable has converged + AssertThrow(userInputs.nonlinear_solver_parameters.getToleranceType(var_index) == + ABSOLUTE_SOLUTION_CHANGE, + FeatureNotImplemented( + "Nonlinear solver tolerances besides ABSOLUTE_CHANGE")); + double diff = NAN; + // if (fields[var_index].type == SCALAR) + // { + // dU_scalar -= *solutionSet[var_index]; + // diff = dU_scalar.l2_norm(); + // } + // else + // { + // dU_vector -= *solutionSet[var_index]; + // diff = dU_vector.l2_norm(); + // } + print_nonlinear_update(var_index, current_increment); + return !(diff > userInputs.nonlinear_solver_parameters.getToleranceValue(var_index)); +} + +// Application of boundary conditions +template +void +MatrixFreePDE::applyBCs(unsigned int var_index) +{ + // Add Neumann BCs + // if (fields[var_index].hasNeumannBCs) + // { + // // Currently commented out because it isn't working yet + // // applyNeumannBCs(); + // } + + // Set the Dirichelet values (hanging node constraints don't need to be distributed + // every time step, only at output) + // if (fields[var_index].hasDirichletBCs) + // { + // // Apply non-uniform Dirlichlet_BCs to the current field + // if (fields[var_index].hasnonuniformDirichletBCs) + // { + // // DoFHandler *dof_handler = + // dofHandlersSet_nonconst.at(currentvar_index); + // // IndexSet *locally_relevant_dofs = + // // locally_relevant_dofsSet_nonconst.at(currentvar_index); + // // locally_relevant_dofs->clear(); + // // DoFTools::extract_locally_relevant_dofs(*dof_handler, + // // *locally_relevant_dofs); AffineConstraints *constraintsDirichlet = + // // constraintsDirichletSet_nonconst.at(currentvar_index); + // // constraintsDirichlet->clear(); + // // constraintsDirichlet->reinit(*locally_relevant_dofs); + // // applyDirichletBCs(); + // // constraintsDirichlet->close(); + // } + // // Distribute for Uniform or Non-Uniform Dirichlet BCs + // constraintsDirichletSet[var_index]->distribute(*solutionSet[var_index]); + // } + solutionSet[var_index]->update_ghost_values(); +} + +// Explicit time step for matrixfree solve +template +void +MatrixFreePDE::updateExplicitSolution(unsigned int var_index) +{ + // Explicit-time step each DOF + // Takes advantage of knowledge that the length of solutionSet and residualSet + // is an integer multiple of the length of invM for vector variables + // if (fields[var_index].type == SCALAR) + // { + // unsigned int invM_size = invMscalar.locally_owned_size(); + // for (unsigned int dof = 0; dof < solutionSet[var_index]->locally_owned_size(); + // ++dof) + // { + // solutionSet[var_index]->local_element(dof) = + // invMscalar.local_element(dof % invM_size) * + // residualSet[var_index]->local_element(dof); + // } + // } + // else if (fields[var_index].type == VECTOR) + // { + // unsigned int invM_size = invMvector.locally_owned_size(); + // for (unsigned int dof = 0; dof < solutionSet[var_index]->locally_owned_size(); + // ++dof) + // { + // solutionSet[var_index]->local_element(dof) = + // invMvector.local_element(dof % invM_size) * + // residualSet[var_index]->local_element(dof); + // } + // } +} + +template +void +MatrixFreePDE::updateLinearSolution(unsigned int var_index) +{ + // if (fields[var_index].type == SCALAR) + // { + // *solutionSet[var_index] += dU_scalar; + // } + // else + // { + // *solutionSet[var_index] += dU_vector; + // } } -template class MatrixFreePDE<2, 1>; -template class MatrixFreePDE<3, 1>; +template +void +MatrixFreePDE::print_explicit_update(uint var_index, uint current_increment) +{ + char buffer[200]; -template class MatrixFreePDE<2, 2>; -template class MatrixFreePDE<3, 2>; + if (current_increment % userInputs.skip_print_steps == 0) + { + double solution_L2_norm = solutionSet[var_index]->l2_norm(); + + // snprintf(buffer, + // sizeof(buffer), + // "field '%2s' [explicit solve]: current solution: " + // "%12.6e, current residual:%12.6e\n", + // fields[var_index].name.c_str(), + // solution_L2_norm, + // residualSet[var_index]->l2_norm()); + // conditionalOStreams::pout_base << buffer; + + if (!numbers::is_finite(solution_L2_norm)) + { + // snprintf(buffer, + // sizeof(buffer), + // "ERROR: field '%s' solution is NAN. exiting.\n\n", + // fields[var_index].name.c_str()); + // conditionalOStreams::pout_base << buffer; + exit(-1); + } + } +} -template class MatrixFreePDE<3, 3>; -template class MatrixFreePDE<2, 3>; +template +void +MatrixFreePDE::print_linear_update(uint var_index, + uint current_increment, + SolverControl &solver_control) +{ + if (current_increment % userInputs.skip_print_steps == 0) + { + double dU_norm = NAN; + // if (fields[var_index].type == SCALAR) + // { + // dU_norm = dU_scalar.l2_norm(); + // } + // else + // { + // dU_norm = dU_vector.l2_norm(); + // } + char buffer[200]; + // snprintf(buffer, + // sizeof(buffer), + // "field '%2s' [linear solve]: initial " + // "residual:%12.6e, current residual:%12.6e, " + // "nsteps:%u, tolerance criterion:%12.6e, " + // "solution: %12.6e, dU: %12.6e\n", + // fields[var_index].name.c_str(), + // residualSet[var_index]->l2_norm(), + // solver_control.last_value(), + // solver_control.last_step(), + // solver_control.tolerance(), + // solutionSet[var_index]->l2_norm(), + // dU_norm); + conditionalOStreams::pout_base << buffer; + } +} -template class MatrixFreePDE<3, 4>; -template class MatrixFreePDE<2, 4>; +template +void +MatrixFreePDE::print_nonlinear_update(uint var_index, uint current_increment) +{ + char buffer[200]; -template class MatrixFreePDE<3, 5>; -template class MatrixFreePDE<2, 5>; + if (!(current_increment % userInputs.skip_print_steps)) + { + // snprintf(buffer, + // sizeof(buffer), + // "field '%2s' [nonlinear solve]: current " + // "solution: %12.6e, current residual:%12.6e\n", + // fields[var_index].name.c_str(), + // solutionSet[var_index]->l2_norm(), + // residualSet[var_index]->l2_norm()); + conditionalOStreams::pout_base << buffer; + } +} -template class MatrixFreePDE<3, 6>; -template class MatrixFreePDE<2, 6>; \ No newline at end of file +template +void +MatrixFreePDE::print_nonlinear_status(uint var_index, + uint current_increment, + uint nonlinear_iteration_index, + double diff) +{ + if (!(current_increment % userInputs.skip_print_steps)) + { + char buffer[200]; + // snprintf(buffer, + // sizeof(buffer), + // " field '%2s' [nonlinear solve] current increment: %u, nonlinear " + // "iteration: " + // "%u, dU: %12.6e\n", + // fields[var_index].name.c_str(), + // current_increment, + // nonlinear_iteration_index, + // diff); + conditionalOStreams::pout_base << buffer; + } +} \ No newline at end of file diff --git a/src/core/userInputParameters.cc b/src/core/userInputParameters.cc index 10e50357..78799bc9 100644 --- a/src/core/userInputParameters.cc +++ b/src/core/userInputParameters.cc @@ -1,6 +1,7 @@ #include #include +#include #include #include @@ -75,7 +76,7 @@ userInputParameters::loadVariableAttributes() varInfoListExplicitRHS.reserve(num_var_explicit_RHS); for (const auto &[index, variable] : var_attributes) { - variable_info varInfo {}; + variableInfo varInfo {}; varInfo.evaluation_flags = variable.eval_flags_explicit_RHS; @@ -83,10 +84,7 @@ userInputParameters::loadVariableAttributes() varInfo.global_var_index = index; - varInfo.var_needed = - !static_cast(varInfo.evaluation_flags & dealii::EvaluationFlags::nothing); - - varInfo.is_scalar = variable.var_type == SCALAR; + varInfo.var_type = variable.var_type; varInfoListExplicitRHS.push_back(varInfo); } @@ -104,7 +102,7 @@ userInputParameters::loadVariableAttributes() varInfoListNonexplicitRHS.reserve(num_var_nonexplicit_RHS); for (const auto &[index, variable] : var_attributes) { - variable_info varInfo {}; + variableInfo varInfo {}; varInfo.evaluation_flags = variable.eval_flags_nonexplicit_RHS; @@ -112,10 +110,7 @@ userInputParameters::loadVariableAttributes() varInfo.global_var_index = index; - varInfo.var_needed = - !static_cast(varInfo.evaluation_flags & dealii::EvaluationFlags::nothing); - - varInfo.is_scalar = variable.var_type == SCALAR; + varInfo.var_type = variable.var_type; varInfoListNonexplicitRHS.push_back(varInfo); } @@ -134,7 +129,7 @@ userInputParameters::loadVariableAttributes() varInfoListLHS.reserve(num_var_LHS); for (const auto &[index, variable] : var_attributes) { - variable_info varInfo {}; + variableInfo varInfo {}; varInfo.evaluation_flags = variable.eval_flags_nonexplicit_LHS; @@ -142,10 +137,7 @@ userInputParameters::loadVariableAttributes() varInfo.global_var_index = index; - varInfo.var_needed = - !static_cast(varInfo.evaluation_flags & dealii::EvaluationFlags::nothing); - - varInfo.is_scalar = variable.var_type == SCALAR; + varInfo.var_type = variable.var_type; varInfoListLHS.push_back(varInfo); } @@ -153,7 +145,7 @@ userInputParameters::loadVariableAttributes() varChangeInfoListLHS.reserve(num_var_LHS); for (const auto &[index, variable] : var_attributes) { - variable_info varInfo {}; + variableInfo varInfo {}; varInfo.evaluation_flags = variable.eval_flags_change_nonexplicit_LHS; @@ -162,10 +154,7 @@ userInputParameters::loadVariableAttributes() varInfo.global_var_index = index; - varInfo.var_needed = - !static_cast(varInfo.evaluation_flags & dealii::EvaluationFlags::nothing); - - varInfo.is_scalar = variable.var_type == SCALAR; + varInfo.var_type = variable.var_type; varChangeInfoListLHS.push_back(varInfo); } @@ -175,16 +164,13 @@ userInputParameters::loadVariableAttributes() pp_baseVarInfoList.reserve(var_attributes.size()); for (const auto &[index, variable] : var_attributes) { - variable_info varInfo {}; + variableInfo varInfo {}; varInfo.evaluation_flags = variable.eval_flags_postprocess; varInfo.global_var_index = index; - varInfo.var_needed = - !static_cast(varInfo.evaluation_flags & dealii::EvaluationFlags::nothing); - - varInfo.is_scalar = variable.var_type == SCALAR; + varInfo.var_type = variable.var_type; pp_baseVarInfoList.push_back(varInfo); } @@ -208,15 +194,13 @@ userInputParameters::loadVariableAttributes() pp_varInfoList.reserve(pp_attributes.size()); for (const auto &[pp_index, pp_variable] : pp_attributes) { - variable_info varInfo {}; - - varInfo.var_needed = true; + variableInfo varInfo {}; varInfo.residual_flags = pp_variable.eval_flags_residual_postprocess; varInfo.global_var_index = pp_index; - varInfo.is_scalar = pp_variable.var_type == SCALAR; + varInfo.var_type = pp_variable.var_type; pp_varInfoList.push_back(varInfo); } @@ -541,15 +525,13 @@ userInputParameters::assign_nonlinear_solve_parameters( } else { - if (dealii::Utilities::MPI::this_mpi_process(MPI_COMM_WORLD) == 0) - { - std::cout << "PRISMS-PF Warning: Laplace's equation is only used " - "to generate the initial guess for time independent " - "equations. The equation for variable " - << variable.name - << " is not a time independent equation. No initial " - "guess is needed for this equation.\n"; - } + conditionalOStreams::pout_base + << "PRISMS-PF Warning: Laplace's equation is only used " + "to generate the initial guess for time independent " + "equations. The equation for variable " + << variable.name + << " is not a time independent equation. No initial " + "guess is needed for this equation.\n"; } nonlinear_solver_parameters.loadParameters(index, @@ -603,13 +585,12 @@ userInputParameters::assign_output_parameters( if ((output_file_type == "vtk") && (!output_vtu_per_process)) { output_vtu_per_process = true; - if (dealii::Utilities::MPI::this_mpi_process(MPI_COMM_WORLD) == 0) - { - std::cout << "PRISMS-PF Warning: 'Output file type' given as 'vtk' and " - "'Output separate files per process' given as 'false'. Shared " - "output files are not supported for the vtk output format. " - "Separate files per process will be created.\n"; - } + + conditionalOStreams::pout_base + << "PRISMS-PF Warning: 'Output file type' given as 'vtk' and " + "'Output separate files per process' given as 'false'. Shared " + "output files are not supported for the vtk output format. " + "Separate files per process will be created.\n"; } print_timing_with_output = @@ -1319,7 +1300,7 @@ userInputParameters::load_model_constants( } template -dealii::Tensor<2, 2 * dim - 1 + dim / 3> +dealii::Tensor<2, (2 * dim) - 1 + (dim / 3)> userInputParameters::get_Cij_tensor(std::vector elastic_constants, const std::string &elastic_const_symmetry) const { @@ -1362,34 +1343,29 @@ userInputParameters::get_Cij_tensor(std::vector elastic_constants, } } - dealii::ConditionalOStream pcout(std::cout, - dealii::Utilities::MPI::this_mpi_process( - MPI_COMM_WORLD) == 0); - - return getCIJMatrix(mat_model, elastic_constants, pcout); + return getCIJMatrix(mat_model, elastic_constants); } template -dealii::Tensor<2, 2 * dim - 1 + dim / 3> -userInputParameters::getCIJMatrix(const elasticityModel model, - const std::vector &constants, - dealii::ConditionalOStream &pcout) const +dealii::Tensor<2, (2 * dim) - 1 + (dim / 3)> +userInputParameters::getCIJMatrix(const elasticityModel model, + const std::vector &constants) const { // CIJ.fill(0.0); dealii::Tensor<2, 2 * dim - 1 + dim / 3> CIJ; - pcout << "Reading material model:"; + conditionalOStreams::pout_base << "Reading material model:"; switch (dim) { case 1: { - pcout << " 1D "; + conditionalOStreams::pout_base << " 1D "; // 1D models switch (model) { case ISOTROPIC: { - pcout << " ISOTROPIC \n"; + conditionalOStreams::pout_base << " ISOTROPIC \n"; CIJ[0][0] = constants[0]; break; } @@ -1405,13 +1381,13 @@ userInputParameters::getCIJMatrix(const elasticityModel model, } case 2: { - pcout << " 2D "; + conditionalOStreams::pout_base << " 2D "; // 2D models switch (model) { case ISOTROPIC: { - pcout << " ISOTROPIC \n"; + conditionalOStreams::pout_base << " ISOTROPIC \n"; const double E = constants[0]; const double nu = constants[1]; const double mu = E / (2 * (1 + nu)); @@ -1424,7 +1400,7 @@ userInputParameters::getCIJMatrix(const elasticityModel model, } case ANISOTROPIC: { - pcout << " ANISOTROPIC \n"; + conditionalOStreams::pout_base << " ANISOTROPIC \n"; CIJ[0][0] = constants[0]; // C11 CIJ[1][1] = constants[1]; // C22 CIJ[2][2] = constants[2]; // C33 @@ -1445,13 +1421,13 @@ userInputParameters::getCIJMatrix(const elasticityModel model, } case 3: { - pcout << " 3D "; + conditionalOStreams::pout_base << " 3D "; // 3D models switch (model) { case ISOTROPIC: { - pcout << " ISOTROPIC \n"; + conditionalOStreams::pout_base << " ISOTROPIC \n"; const double E = constants[0]; const double nu = constants[1]; const double mu = E / (2 * (1 + nu)); @@ -1469,7 +1445,7 @@ userInputParameters::getCIJMatrix(const elasticityModel model, } case TRANSVERSE: { - pcout << " TRANSVERSE \n"; + conditionalOStreams::pout_base << " TRANSVERSE \n"; CIJ[0][0] = constants[0]; // C11 CIJ[1][1] = constants[0]; // C11 CIJ[2][2] = constants[1]; // C33 @@ -1483,7 +1459,7 @@ userInputParameters::getCIJMatrix(const elasticityModel model, } case ORTHOTROPIC: { - pcout << " ORTHOTROPIC \n"; + conditionalOStreams::pout_base << " ORTHOTROPIC \n"; CIJ[0][0] = constants[0]; // C11 CIJ[1][1] = constants[1]; // C22 CIJ[2][2] = constants[2]; // C33 @@ -1497,7 +1473,7 @@ userInputParameters::getCIJMatrix(const elasticityModel model, } case ANISOTROPIC: { - pcout << " ANISOTROPIC \n"; + conditionalOStreams::pout_base << " ANISOTROPIC \n"; CIJ[0][0] = constants[0]; // C11 CIJ[1][1] = constants[1]; // C22 CIJ[2][2] = constants[2]; // C33 @@ -1538,18 +1514,18 @@ userInputParameters::getCIJMatrix(const elasticityModel model, } } // print CIJ to terminal - pcout << "Elasticity matrix (Voigt notation):\n"; - constexpr unsigned int voight_matrix_size = 2 * dim - 1 + dim / 3; + conditionalOStreams::pout_base << "Elasticity matrix (Voigt notation):\n"; + constexpr unsigned int voight_matrix_size = (2 * dim) - 1 + (dim / 3); for (unsigned int i = 0; i < voight_matrix_size; i++) { for (unsigned int j = 0; j < voight_matrix_size; j++) { - pcout << std::setw(8) << std::setprecision(3) << std::scientific << CIJ[i][j] - << " "; + conditionalOStreams::pout_base << std::setw(8) << std::setprecision(3) + << std::scientific << CIJ[i][j] << " "; } - pcout << "\n"; + conditionalOStreams::pout_base << "\n"; } - pcout << "\n"; + conditionalOStreams::pout_base << "\n"; return CIJ; } diff --git a/src/core/variableContainer.cc b/src/core/variableContainer.cc index 7d5bcf45..bacf3527 100644 --- a/src/core/variableContainer.cc +++ b/src/core/variableContainer.cc @@ -1,120 +1,60 @@ +#include +#include + +#include #include +#include #include +#include +#include -template -variableContainer::variableContainer( - const dealii::MatrixFree &data, - const std::vector &_varInfoList, - const std::vector &_varChangeInfoList) - : varInfoList(_varInfoList) - , varChangeInfoList(_varChangeInfoList) - , num_var(varInfoList.size()) -{ - for (unsigned int i = 0; i < num_var; i++) - { - const auto &var_info = varInfoList[i]; - const auto &var_change_info = varChangeInfoList[i]; +template +variableContainer::variableContainer( + const dealii::MatrixFree &data) - if (var_info.var_needed) - { - const unsigned int var_index = var_info.global_var_index; - - if (var_info.is_scalar) - { - scalar_vars_map.emplace(var_index, - std::make_unique(data, i)); - } - else - { - vector_vars_map.emplace(var_index, - std::make_unique(data, i)); - } - } +{ + // Set all flags here for testing purposes. Also only one scalar field + variableInfo var_info {0, + fieldType::SCALAR, + dealii::EvaluationFlags::gradients, + dealii::EvaluationFlags::gradients}; - if (var_change_info.var_needed) - { - const unsigned int var_index = var_change_info.global_var_index; - - if (var_change_info.is_scalar) - { - scalar_change_in_vars_map.emplace(var_index, - std::make_unique(data, i)); - } - else - { - vector_change_in_vars_map.emplace(var_index, - std::make_unique(data, i)); - } - } - } -} + varInfoList.emplace_back(var_info); + num_var = varInfoList.size(); -template -variableContainer::variableContainer( - const dealii::MatrixFree &data, - const std::vector &_varInfoList) - : varInfoList(_varInfoList) - , num_var(varInfoList.size()) -{ for (unsigned int i = 0; i < num_var; i++) { const auto &var_info = varInfoList[i]; - if (!var_info.var_needed) - { - continue; - } - const unsigned int var_index = var_info.global_var_index; - if (var_info.is_scalar) + if (var_info.var_type == fieldType::SCALAR) { scalar_vars_map.emplace(var_index, std::make_unique(data, i)); } - else + else if (var_info.var_type == fieldType::VECTOR) { vector_vars_map.emplace(var_index, std::make_unique(data, i)); } } } -// Variant of the constructor where it reads from a fixed index of "data", used -// for post-processing -template -variableContainer::variableContainer( - const dealii::MatrixFree &data, - const std::vector &_varInfoList, - const unsigned int &fixed_index) - : varInfoList(_varInfoList) - , num_var(varInfoList.size()) +template +void +variableContainer::access_valid( + [[maybe_unused]] const unsigned int &global_variable_index, + [[maybe_unused]] const EvalFlags &flag) const { - for (unsigned int i = 0; i < num_var; i++) - { - const auto &var_info = varInfoList[i]; - - if (!var_info.var_needed) - { - continue; - } - - const unsigned int var_index = var_info.global_var_index; - - if (var_info.is_scalar) - { - scalar_vars_map.emplace(var_index, - std::make_unique(data, fixed_index)); - } - else - { - vector_vars_map.emplace(var_index, - std::make_unique(data, fixed_index)); - } - } + Assert(varInfoList.at(global_variable_index).evaluation_flags & flag, + dealii::ExcMessage("PRISMS-PF Error: Attempted access of a variable that was " + "not marked as needed in 'equations.cc'. The attempted " + "access was for variable with index " + + std::to_string(global_variable_index) + ".")); } -template +template unsigned int -variableContainer::get_num_q_points() const +variableContainer::get_num_q_points() const { if (!scalar_vars_map.empty()) { @@ -124,26 +64,17 @@ variableContainer::get_num_q_points() const { return vector_vars_map.begin()->second->n_q_points; } - else if (!scalar_change_in_vars_map.empty()) - { - return scalar_change_in_vars_map.begin()->second->n_q_points; - } - else if (!vector_change_in_vars_map.empty()) - { - return vector_change_in_vars_map.begin()->second->n_q_points; - } - else - { - AssertThrow(false, - dealii::ExcMessage( - "PRISMS-PF Error: When trying to access the number of quadrature " - "points, all FEEvaluation object containers were empty.")); - } + + AssertThrow(false, + dealii::ExcMessage( + "PRISMS-PF Error: When trying to access the number of quadrature " + "points, all FEEvaluation object containers were empty.")); + return 0; } -template -dealii::Point -variableContainer::get_q_point_location() const +template +dealii::Point> +variableContainer::get_q_point_location() const { if (!scalar_vars_map.empty()) { @@ -153,493 +84,160 @@ variableContainer::get_q_point_location() const { return vector_vars_map.begin()->second->quadrature_point(q_point); } - else if (!scalar_change_in_vars_map.empty()) - { - return scalar_change_in_vars_map.begin()->second->quadrature_point(q_point); - } - else if (!vector_change_in_vars_map.empty()) - { - return vector_change_in_vars_map.begin()->second->quadrature_point(q_point); - } - else - { - AssertThrow(false, - dealii::ExcMessage( - "PRISMS-PF Error: When trying to access the quadrature point " - "location, all FEEvaluation object containers were empty.")); - } + + AssertThrow(false, + dealii::ExcMessage( + "PRISMS-PF Error: When trying to access the quadrature point " + "location, all FEEvaluation object containers were empty.")); + return dealii::Point>(); } -template +template void -variableContainer::reinit_and_eval( - const std::vector *> &src, - unsigned int cell) +variableContainer::reinit_and_eval( + const dealii::LinearAlgebra::distributed::Vector &src, + unsigned int cell) { for (unsigned int i = 0; i < num_var; i++) { - const auto &var_info = varInfoList[i]; - - if (!var_info.var_needed) - { - continue; - } - + const auto &var_info = varInfoList[i]; const unsigned int var_index = var_info.global_var_index; - if (var_info.is_scalar) + if (var_info.var_type == fieldType::SCALAR) { - auto *scalar_FEEval_ptr = scalar_vars_map[var_index].get(); + auto *scalar_FEEval_ptr = scalar_vars_map.at(var_index).get(); scalar_FEEval_ptr->reinit(cell); - scalar_FEEval_ptr->read_dof_values(*src[i]); + scalar_FEEval_ptr->read_dof_values(src); scalar_FEEval_ptr->evaluate(var_info.evaluation_flags); } - else + else if (var_info.var_type == fieldType::VECTOR) { - auto *vector_FEEval_ptr = vector_vars_map[var_index].get(); + auto *vector_FEEval_ptr = vector_vars_map.at(var_index).get(); vector_FEEval_ptr->reinit(cell); - vector_FEEval_ptr->read_dof_values(*src[i]); + vector_FEEval_ptr->read_dof_values(src); vector_FEEval_ptr->evaluate(var_info.evaluation_flags); } } } -/** - * This is specialized for the LHS where a change in solution is needed. The RHS method - * takes the src as a vector of dealii::LinearAlgebra::distributed::Vectors. - */ -template -void -variableContainer::reinit_and_eval_change_in_solution( - const dealii::LinearAlgebra::distributed::Vector &src, - unsigned int cell, - unsigned int var_being_solved) -{ - if (varChangeInfoList[var_being_solved].is_scalar) - { - auto *scalar_FEEval_ptr = scalar_change_in_vars_map[var_being_solved].get(); - scalar_FEEval_ptr->reinit(cell); - scalar_FEEval_ptr->read_dof_values(src); - scalar_FEEval_ptr->evaluate(varChangeInfoList[var_being_solved].evaluation_flags); - } - else - { - auto *vector_FEEval_ptr = vector_change_in_vars_map[var_being_solved].get(); - vector_FEEval_ptr->reinit(cell); - vector_FEEval_ptr->read_dof_values(src); - vector_FEEval_ptr->evaluate(varChangeInfoList[var_being_solved].evaluation_flags); - } -} - -template +template void -variableContainer::reinit(unsigned int cell) +variableContainer::integrate_and_distribute( + dealii::LinearAlgebra::distributed::Vector &dst) { for (unsigned int i = 0; i < num_var; i++) { const auto &var_info = varInfoList[i]; - if (!var_info.var_needed) - { - continue; - } - const unsigned int var_index = var_info.global_var_index; - if (var_info.is_scalar) + if (var_info.var_type == fieldType::SCALAR) { - scalar_vars_map[var_index]->reinit(cell); + auto *scalar_FEEval_ptr = scalar_vars_map.at(var_index).get(); + scalar_FEEval_ptr->integrate_scatter(var_info.residual_flags, dst); } - else + else if (var_info.var_type == fieldType::VECTOR) { - vector_vars_map[var_index]->reinit(cell); + auto *vector_FEEval_ptr = vector_vars_map.at(var_index).get(); + vector_FEEval_ptr->integrate_scatter(var_info.residual_flags, dst); } } } -template -void -variableContainer::integrate_and_distribute( - std::vector *> &dst) -{ - for (unsigned int i = 0; i < num_var; i++) - { - const auto &var_info = varInfoList[i]; - - if (var_info.residual_flags & dealii::EvaluationFlags::nothing) - { - continue; - } - - const unsigned int var_index = var_info.global_var_index; - - if (var_info.is_scalar) - { - auto *scalar_FEEval_ptr = scalar_vars_map[var_index].get(); - scalar_FEEval_ptr->integrate(var_info.residual_flags); - scalar_FEEval_ptr->distribute_local_to_global(*dst[i]); - } - else - { - auto *vector_FEEval_ptr = vector_vars_map[var_index].get(); - vector_FEEval_ptr->integrate(var_info.residual_flags); - vector_FEEval_ptr->distribute_local_to_global(*dst[i]); - } - } -} - -template -void -variableContainer::integrate_and_distribute_change_in_solution_LHS( - dealii::LinearAlgebra::distributed::Vector &dst, - const unsigned int var_being_solved) -{ - // integrate - if (varChangeInfoList[var_being_solved].is_scalar) - { - auto *scalar_FEEval_ptr = scalar_change_in_vars_map[var_being_solved].get(); - scalar_FEEval_ptr->integrate(varChangeInfoList[var_being_solved].residual_flags); - scalar_FEEval_ptr->distribute_local_to_global(dst); - } - else - { - auto *vector_FEEval_ptr = vector_change_in_vars_map[var_being_solved].get(); - vector_FEEval_ptr->integrate(varChangeInfoList[var_being_solved].residual_flags); - vector_FEEval_ptr->distribute_local_to_global(dst); - } -} - -// Need to add index checking to these functions so that an error is thrown if -// the index wasn't set -template -T -variableContainer::get_scalar_value( +template +dealii::VectorizedArray +variableContainer::get_scalar_value( unsigned int global_variable_index) const { - if (varInfoList[global_variable_index].evaluation_flags & - dealii::EvaluationFlags::values) - { - return scalar_vars_map.at(global_variable_index)->get_value(q_point); - } - else - { - std::cerr << "PRISMS-PF Error: Attempted access of a variable value that " - "was not marked as needed in 'equations.cc'. The attempted " - "access was for variable with index " - << global_variable_index << " .\n"; - abort(); - } -} + access_valid(global_variable_index, EvalFlags::values); -template -dealii::Tensor<1, dim, T> -variableContainer::get_scalar_gradient( - unsigned int global_variable_index) const -{ - if (varInfoList[global_variable_index].evaluation_flags & - dealii::EvaluationFlags::gradients) - { - return scalar_vars_map.at(global_variable_index)->get_gradient(q_point); - } - else - { - std::cerr << "PRISMS-PF Error: Attempted access of a variable gradient " - "that was not marked as needed in 'equations.cc'. The " - "attempted access was for variable with index " - << global_variable_index << " .\n"; - abort(); - } + return scalar_vars_map.at(global_variable_index)->get_value(q_point); } -template -dealii::Tensor<2, dim, T> -variableContainer::get_scalar_hessian( +template +dealii::Tensor<1, dim, dealii::VectorizedArray> +variableContainer::get_scalar_gradient( unsigned int global_variable_index) const { - if (varInfoList[global_variable_index].evaluation_flags & - dealii::EvaluationFlags::hessians) - { - return scalar_vars_map.at(global_variable_index)->get_hessian(q_point); - } - else - { - std::cerr << "PRISMS-PF Error: Attempted access of a variable hessian " - "that was not marked as needed in 'equations.cc'. The " - "attempted access was for variable with index " - << global_variable_index << " .\n"; - abort(); - } -} + access_valid(global_variable_index, EvalFlags::gradients); -template -dealii::Tensor<1, dim, T> -variableContainer::get_vector_value( - unsigned int global_variable_index) const -{ - if (varInfoList[global_variable_index].evaluation_flags & - dealii::EvaluationFlags::values) - { - return vector_vars_map.at(global_variable_index)->get_value(q_point); - } - else - { - std::cerr << "PRISMS-PF Error: Attempted access of a variable value that " - "was not marked as needed in 'equations.cc'. The attempted " - "access was for variable with index " - << global_variable_index << " .\n"; - abort(); - } + return scalar_vars_map.at(global_variable_index)->get_gradient(q_point); } -template -dealii::Tensor<2, dim, T> -variableContainer::get_vector_gradient( +template +dealii::Tensor<2, dim, dealii::VectorizedArray> +variableContainer::get_scalar_hessian( unsigned int global_variable_index) const { - if (varInfoList[global_variable_index].evaluation_flags & - dealii::EvaluationFlags::gradients) - { - return vector_vars_map.at(global_variable_index)->get_gradient(q_point); - } - else - { - std::cerr << "PRISMS-PF Error: Attempted access of a variable gradient " - "that was not marked as needed in 'equations.cc'. The " - "attempted access was for variable with index " - << global_variable_index << " .\n"; - abort(); - } -} + access_valid(global_variable_index, EvalFlags::hessians); -template -dealii::Tensor<3, dim, T> -variableContainer::get_vector_hessian( - unsigned int global_variable_index) const -{ - if (varInfoList[global_variable_index].evaluation_flags & - dealii::EvaluationFlags::hessians) - { - return vector_vars_map.at(global_variable_index)->get_hessian(q_point); - } - else - { - std::cerr << "PRISMS-PF Error: Attempted access of a variable hessian " - "that was not marked as needed in 'equations.cc'. The " - "attempted access was for variable with index " - << global_variable_index << " .\n"; - abort(); - } + return scalar_vars_map.at(global_variable_index)->get_hessian(q_point); } -// Need to add index checking to these functions so that an error is thrown if -// the index wasn't set -template -T -variableContainer::get_change_in_scalar_value( +template +dealii::Tensor<1, dim, dealii::VectorizedArray> +variableContainer::get_vector_value( unsigned int global_variable_index) const { - if (varChangeInfoList[global_variable_index].evaluation_flags & - dealii::EvaluationFlags::values) - { - return scalar_change_in_vars_map.at(global_variable_index)->get_value(q_point); - } - else - { - std::cerr << "PRISMS-PF Error: Attempted access of a change in variable " - "value that was not marked as needed in 'equations.cc'. The " - "attempted access was for variable with index " - << global_variable_index << " .\n"; - abort(); - } -} - -template -dealii::Tensor<1, dim, T> -variableContainer::get_change_in_scalar_gradient( - unsigned int global_variable_index) const -{ - if (varChangeInfoList[global_variable_index].evaluation_flags & - dealii::EvaluationFlags::gradients) - { - return scalar_change_in_vars_map.at(global_variable_index)->get_gradient(q_point); - } - else - { - std::cerr << "PRISMS-PF Error: Attempted access of a change in variable " - "gradient that was not marked as needed in 'equations.cc'. " - "The attempted access was for variable with index " - << global_variable_index << " .\n"; - abort(); - } -} + access_valid(global_variable_index, EvalFlags::values); -template -dealii::Tensor<2, dim, T> -variableContainer::get_change_in_scalar_hessian( - unsigned int global_variable_index) const -{ - if (varChangeInfoList[global_variable_index].evaluation_flags & - dealii::EvaluationFlags::hessians) - { - return scalar_change_in_vars_map.at(global_variable_index)->get_hessian(q_point); - } - else - { - std::cerr << "PRISMS-PF Error: Attempted access of a change in variable " - "hessian that was not marked as needed in 'equations.cc'. " - "The attempted access was for variable with index " - << global_variable_index << " .\n"; - abort(); - } + return vector_vars_map.at(global_variable_index)->get_value(q_point); } -template -dealii::Tensor<1, dim, T> -variableContainer::get_change_in_vector_value( +template +dealii::Tensor<2, dim, dealii::VectorizedArray> +variableContainer::get_vector_gradient( unsigned int global_variable_index) const { - if (varChangeInfoList[global_variable_index].evaluation_flags & - dealii::EvaluationFlags::values) - { - return vector_change_in_vars_map.at(global_variable_index)->get_value(q_point); - } - else - { - std::cerr << "PRISMS-PF Error: Attempted access of a change in variable " - "value that was not marked as needed in 'equations.cc'. The " - "attempted access was for variable with index " - << global_variable_index << " .\n"; - abort(); - } -} + access_valid(global_variable_index, EvalFlags::gradients); -template -dealii::Tensor<2, dim, T> -variableContainer::get_change_in_vector_gradient( - unsigned int global_variable_index) const -{ - if (varChangeInfoList[global_variable_index].evaluation_flags & - dealii::EvaluationFlags::gradients) - { - return vector_change_in_vars_map.at(global_variable_index)->get_gradient(q_point); - } - else - { - std::cerr << "PRISMS-PF Error: Attempted access of a change in variable " - "gradient that was not marked as needed in 'equations.cc'. " - "The attempted access was for variable with index " - << global_variable_index << " .\n"; - abort(); - } + return vector_vars_map.at(global_variable_index)->get_gradient(q_point); } -template -dealii::Tensor<3, dim, T> -variableContainer::get_change_in_vector_hessian( +template +dealii::Tensor<3, dim, dealii::VectorizedArray> +variableContainer::get_vector_hessian( unsigned int global_variable_index) const { - if (varChangeInfoList[global_variable_index].evaluation_flags & - dealii::EvaluationFlags::hessians) - { - return vector_change_in_vars_map.at(global_variable_index)->get_hessian(q_point); - } - else - { - std::cerr << "PRISMS-PF Error: Attempted access of a change in variable " - "hessian that was not marked as needed in 'equations.cc'. " - "The attempted access was for variable with index " - << global_variable_index << " .\n"; - abort(); - } -} + access_valid(global_variable_index, EvalFlags::hessians); -// The methods to set the residual terms -template -void -variableContainer::set_scalar_value_term_RHS( - unsigned int global_variable_index, - T val) -{ - scalar_vars_map[global_variable_index]->submit_value(val, q_point); + return vector_vars_map.at(global_variable_index)->get_hessian(q_point); } -template +template void -variableContainer::set_scalar_gradient_term_RHS( - unsigned int global_variable_index, - dealii::Tensor<1, dim, T> grad) +variableContainer::set_scalar_value_term( + unsigned int global_variable_index, + dealii::VectorizedArray val) { - scalar_vars_map[global_variable_index]->submit_gradient(grad, q_point); + scalar_vars_map.at(global_variable_index)->submit_value(val, q_point); } -template +template void -variableContainer::set_vector_value_term_RHS( - unsigned int global_variable_index, - dealii::Tensor<1, dim, T> val) +variableContainer::set_scalar_gradient_term( + unsigned int global_variable_index, + dealii::Tensor<1, dim, dealii::VectorizedArray> grad) { - vector_vars_map[global_variable_index]->submit_value(val, q_point); + scalar_vars_map.at(global_variable_index)->submit_gradient(grad, q_point); } -template +template void -variableContainer::set_vector_gradient_term_RHS( - unsigned int global_variable_index, - dealii::Tensor<2, dim, T> grad) +variableContainer::set_vector_value_term( + unsigned int global_variable_index, + dealii::Tensor<1, dim, dealii::VectorizedArray> val) { - vector_vars_map[global_variable_index]->submit_gradient(grad, q_point); + vector_vars_map.at(global_variable_index)->submit_value(val, q_point); } -template +template void -variableContainer::set_scalar_value_term_LHS( - unsigned int global_variable_index, - T val) +variableContainer::set_vector_gradient_term( + unsigned int global_variable_index, + dealii::Tensor<2, dim, dealii::VectorizedArray> grad) { - scalar_change_in_vars_map[global_variable_index]->submit_value(val, q_point); + vector_vars_map.at(global_variable_index)->submit_gradient(grad, q_point); } -template -void -variableContainer::set_scalar_gradient_term_LHS( - unsigned int global_variable_index, - dealii::Tensor<1, dim, T> grad) -{ - scalar_change_in_vars_map[global_variable_index]->submit_gradient(grad, q_point); -} - -template -void -variableContainer::set_vector_value_term_LHS( - unsigned int global_variable_index, - dealii::Tensor<1, dim, T> val) -{ - vector_change_in_vars_map[global_variable_index]->submit_value(val, q_point); -} - -template -void -variableContainer::set_vector_gradient_term_LHS( - unsigned int global_variable_index, - dealii::Tensor<2, dim, T> grad) -{ - vector_change_in_vars_map[global_variable_index]->submit_gradient(grad, q_point); -} - -template class variableContainer<2, 1, dealii::VectorizedArray>; -template class variableContainer<3, 1, dealii::VectorizedArray>; - -template class variableContainer<2, 2, dealii::VectorizedArray>; -template class variableContainer<3, 2, dealii::VectorizedArray>; - -template class variableContainer<2, 3, dealii::VectorizedArray>; -template class variableContainer<3, 3, dealii::VectorizedArray>; - -template class variableContainer<2, 4, dealii::VectorizedArray>; -template class variableContainer<3, 4, dealii::VectorizedArray>; - -template class variableContainer<2, 5, dealii::VectorizedArray>; -template class variableContainer<3, 5, dealii::VectorizedArray>; - -template class variableContainer<2, 6, dealii::VectorizedArray>; -template class variableContainer<3, 6, dealii::VectorizedArray>; \ No newline at end of file +INSTANTIATE_TRI_TEMPLATE(variableContainer) \ No newline at end of file diff --git a/src/grains/FloodFiller.cc b/src/grains/FloodFiller.cc index 32e79ac9..b9deb929 100644 --- a/src/grains/FloodFiller.cc +++ b/src/grains/FloodFiller.cc @@ -1,3 +1,4 @@ +#include #include template @@ -410,21 +411,4 @@ FloodFiller::mergeSplitGrains(std::vector> &grain_set } } -// Template instantiations -template class FloodFiller<2, 1>; -template class FloodFiller<3, 1>; - -template class FloodFiller<2, 2>; -template class FloodFiller<3, 2>; - -template class FloodFiller<2, 3>; -template class FloodFiller<3, 3>; - -template class FloodFiller<2, 4>; -template class FloodFiller<3, 4>; - -template class FloodFiller<2, 5>; -template class FloodFiller<3, 5>; - -template class FloodFiller<2, 6>; -template class FloodFiller<3, 6>; +INSTANTIATE_BI_TEMPLATE(FloodFiller) \ No newline at end of file diff --git a/src/grains/SimplifiedGrainRepresentation.cc b/src/grains/SimplifiedGrainRepresentation.cc index b5e31e65..99f31674 100644 --- a/src/grains/SimplifiedGrainRepresentation.cc +++ b/src/grains/SimplifiedGrainRepresentation.cc @@ -1,3 +1,4 @@ +#include #include // ============================================================================ @@ -174,15 +175,11 @@ SimplifiedGrainManipulator::reassignGrains( if ((sum_radii + 2.0 * buffer_distance > center_distance) and (order_parameter_other == order_parameter_base)) { - if (dealii::Utilities::MPI::this_mpi_process(MPI_COMM_WORLD) == 0) - { - std::cout << "Found overlap between grain " - << grain_representations.at(g_base).getGrainId() - << " and grain " - << grain_representations.at(g_other).getGrainId() - << " with order parameter " << order_parameter_base - << std::endl; - } + conditionalOStreams::pout_base + << "Found overlap between grain " + << grain_representations.at(g_base).getGrainId() << " and grain " + << grain_representations.at(g_other).getGrainId() + << " with order parameter " << order_parameter_base << std::endl; grain_representations.at(g_base).setDistanceToNeighbor( center_distance - sum_radii); diff --git a/src/grains/reassignGrains.cc b/src/grains/reassignGrains.cc index 8c09cd2a..3990be62 100644 --- a/src/grains/reassignGrains.cc +++ b/src/grains/reassignGrains.cc @@ -1,3 +1,5 @@ +#include +#include #include #include #include @@ -8,10 +10,7 @@ template void MatrixFreePDE::reassignGrains() { - // log time - computing_timer.enter_subsection("matrixFreePDE: reassignGrains"); - - pcout << "Reassigning grains...\n"; + conditionalOStreams::pout_base << "Reassigning grains...\n"; // Get the index of the first scalar field (used to get the FE object and // DOFHandler) @@ -32,30 +31,30 @@ MatrixFreePDE::reassignGrains() std::vector> grain_sets; unsigned int op_list_index = 0; - for (unsigned int fieldIndex = 0; fieldIndex < fields.size(); fieldIndex++) - { - if (op_list_index < userInputs.variables_for_remapping.size()) - { - if (fieldIndex == userInputs.variables_for_remapping.at(op_list_index)) - { - op_list_index++; - - std::vector> single_OP_grain_sets; - flood_filler.calcGrainSets(*FESet.at(scalar_field_index), - *dofHandlersSet_nonconst.at(scalar_field_index), - solutionSet.at(fieldIndex), - userInputs.order_parameter_threshold, - 1.0 + userInputs.order_parameter_threshold, - 0, - fieldIndex, - single_OP_grain_sets); - - grain_sets.insert(grain_sets.end(), - single_OP_grain_sets.begin(), - single_OP_grain_sets.end()); - } - } - } + // for (unsigned int fieldIndex = 0; fieldIndex < fields.size(); fieldIndex++) + // { + // if (op_list_index < userInputs.variables_for_remapping.size()) + // { + // if (fieldIndex == userInputs.variables_for_remapping.at(op_list_index)) + // { + // op_list_index++; + + // std::vector> single_OP_grain_sets; + // flood_filler.calcGrainSets(*FESet.at(scalar_field_index), + // *dofHandlersSet_nonconst.at(scalar_field_index), + // solutionSet.at(fieldIndex), + // userInputs.order_parameter_threshold, + // 1.0 + userInputs.order_parameter_threshold, + // 0, + // fieldIndex, + // single_OP_grain_sets); + + // grain_sets.insert(grain_sets.end(), + // single_OP_grain_sets.begin(), + // single_OP_grain_sets.end()); + // } + // } + // } // Set the grain indices to unique values for (unsigned int g = 0; g < grain_sets.size(); g++) @@ -71,10 +70,11 @@ MatrixFreePDE::reassignGrains() SimplifiedGrainRepresentation simplified_grain_representation( grain_sets.at(g)); - pcout << "Grain: " << simplified_grain_representation.getGrainId() << " " - << simplified_grain_representation.getOrderParameterId() - << " Center: " << simplified_grain_representation.getCenter()(0) << " " - << simplified_grain_representation.getCenter()(1) << std::endl; + conditionalOStreams::pout_base + << "Grain: " << simplified_grain_representation.getGrainId() << " " + << simplified_grain_representation.getOrderParameterId() + << " Center: " << simplified_grain_representation.getCenter()(0) << " " + << simplified_grain_representation.getCenter()(1) << std::endl; simplified_grain_representations.push_back(simplified_grain_representation); } @@ -93,10 +93,11 @@ MatrixFreePDE::reassignGrains() for (unsigned int g = 0; g < this->simplified_grain_representations.size(); g++) { - pcout << "Grain: " << simplified_grain_representations[g].getGrainId() << " " - << simplified_grain_representations[g].getOrderParameterId() - << " Center: " << simplified_grain_representations[g].getCenter()(0) << " " - << simplified_grain_representations[g].getCenter()(1) << std::endl; + conditionalOStreams::pout_base + << "Grain: " << simplified_grain_representations[g].getGrainId() << " " + << simplified_grain_representations[g].getOrderParameterId() + << " Center: " << simplified_grain_representations[g].getCenter()(0) << " " + << simplified_grain_representations[g].getCenter()(1) << std::endl; } OrderParameterRemapper order_parameter_remapper; @@ -105,26 +106,7 @@ MatrixFreePDE::reassignGrains() *dofHandlersSet_nonconst.at(scalar_field_index), FESet.at(scalar_field_index)->dofs_per_cell); - pcout << "Reassigning grains completed.\n\n"; - - // end log - computing_timer.leave_subsection("matrixFreePDE: reassignGrains"); + conditionalOStreams::pout_base << "Reassigning grains completed.\n\n"; } -template class MatrixFreePDE<2, 1>; -template class MatrixFreePDE<3, 1>; - -template class MatrixFreePDE<2, 2>; -template class MatrixFreePDE<3, 2>; - -template class MatrixFreePDE<3, 3>; -template class MatrixFreePDE<2, 3>; - -template class MatrixFreePDE<3, 4>; -template class MatrixFreePDE<2, 4>; - -template class MatrixFreePDE<3, 5>; -template class MatrixFreePDE<2, 5>; - -template class MatrixFreePDE<3, 6>; -template class MatrixFreePDE<2, 6>; \ No newline at end of file +INSTANTIATE_BI_TEMPLATE(MatrixFreePDE) \ No newline at end of file diff --git a/src/nucleation/nucleation.cc b/src/nucleation/nucleation.cc index 55deab9f..210abe2a 100644 --- a/src/nucleation/nucleation.cc +++ b/src/nucleation/nucleation.cc @@ -1,6 +1,7 @@ -// Methods in MatrixFreePDE to update the list of nuclei #include +#include #include +#include #include #include #include @@ -23,15 +24,15 @@ MatrixFreePDE::updateNucleiList() userInputs.dtValue * (double) currentIncrement <= userInputs.nucleation_end_time) { - computing_timer.enter_subsection("matrixFreePDE: nucleation"); // Apply constraints - for (unsigned int fieldIndex = 0; fieldIndex < fields.size(); fieldIndex++) - { - constraintsDirichletSet[fieldIndex]->distribute( - *solutionSet[fieldIndex]); - constraintsOtherSet[fieldIndex]->distribute(*solutionSet[fieldIndex]); - solutionSet[fieldIndex]->update_ghost_values(); - } + // for (unsigned int fieldIndex = 0; fieldIndex < fields.size(); + // fieldIndex++) + // { + // constraintsDirichletSet[fieldIndex]->distribute( + // *solutionSet[fieldIndex]); + // constraintsOtherSet[fieldIndex]->distribute(*solutionSet[fieldIndex]); + // solutionSet[fieldIndex]->update_ghost_values(); + // } std::vector> new_nuclei; if (currentIncrement == 1 && !userInputs.evolution_before_nucleation) @@ -70,7 +71,6 @@ MatrixFreePDE::updateNucleiList() { refineMeshNearNuclei(new_nuclei); } - computing_timer.leave_subsection("matrixFreePDE: nucleation"); } } } @@ -87,10 +87,12 @@ MatrixFreePDE::getNewNuclei() std::vector> newnuclei; // Get list of prospective new nuclei for the local processor - pcout << "Nucleation attempt for increment " << currentIncrement << "\n"; + conditionalOStreams::pout_base << "Nucleation attempt for increment " + << currentIncrement << "\n"; getLocalNucleiList(newnuclei); - pcout << "nucleation attempt! " << currentTime << " " << currentIncrement << "\n"; + conditionalOStreams::pout_base << "nucleation attempt! " << currentTime << " " + << currentIncrement << "\n"; // Generate global list of new nuclei and resolve conflicts between new nuclei parallelNucleationList new_nuclei_parallel(newnuclei); @@ -587,20 +589,4 @@ MatrixFreePDE::weightedDistanceFromNucleusCenter( return weighted_dist; } -template class MatrixFreePDE<2, 1>; -template class MatrixFreePDE<3, 1>; - -template class MatrixFreePDE<2, 2>; -template class MatrixFreePDE<3, 2>; - -template class MatrixFreePDE<3, 3>; -template class MatrixFreePDE<2, 3>; - -template class MatrixFreePDE<3, 4>; -template class MatrixFreePDE<2, 4>; - -template class MatrixFreePDE<3, 5>; -template class MatrixFreePDE<2, 5>; - -template class MatrixFreePDE<3, 6>; -template class MatrixFreePDE<2, 6>; \ No newline at end of file +INSTANTIATE_BI_TEMPLATE(MatrixFreePDE) \ No newline at end of file diff --git a/src/utilities/CMakeLists.txt b/src/utilities/CMakeLists.txt index be0c80f3..e8104c5d 100644 --- a/src/utilities/CMakeLists.txt +++ b/src/utilities/CMakeLists.txt @@ -1,5 +1,4 @@ # Manually specify files to be included list(APPEND PRISMS_PF_SOURCE_FILES - ${CMAKE_CURRENT_SOURCE_DIR}/utilities.cc ) set(PRISMS_PF_SOURCE_FILES ${PRISMS_PF_SOURCE_FILES} PARENT_SCOPE) \ No newline at end of file diff --git a/src/utilities/utilities.cc b/src/utilities/utilities.cc deleted file mode 100644 index 81ee5234..00000000 --- a/src/utilities/utilities.cc +++ /dev/null @@ -1,37 +0,0 @@ -// utility functions for the MatrixFreePDE class - -#include - -// return index of given field name if exists, else throw error -template -unsigned int -MatrixFreePDE::getFieldIndex(std::string _name) -{ - for (const auto &field : fields) - { - if (field.name.compare(_name) == 0) - { - return field.index; - } - } - pcout << "\nutilities.h: field '" << _name.c_str() << "' not initialized\n"; - exit(-1); -} - -template class MatrixFreePDE<2, 1>; -template class MatrixFreePDE<3, 1>; - -template class MatrixFreePDE<2, 2>; -template class MatrixFreePDE<3, 2>; - -template class MatrixFreePDE<3, 3>; -template class MatrixFreePDE<2, 3>; - -template class MatrixFreePDE<3, 4>; -template class MatrixFreePDE<2, 4>; - -template class MatrixFreePDE<3, 5>; -template class MatrixFreePDE<2, 5>; - -template class MatrixFreePDE<3, 6>; -template class MatrixFreePDE<2, 6>; \ No newline at end of file