From 09016bbed5b1abdc88dd39802a62722f58d893ab Mon Sep 17 00:00:00 2001 From: DavidSCN Date: Sat, 13 Feb 2021 12:13:25 +0100 Subject: [PATCH 01/18] Move adapter files to separate directory --- {adapter => include/adapter}/adapter.h | 3 +-- adapter/time.h => include/adapter/time_handler.h | 6 +++--- 2 files changed, 4 insertions(+), 5 deletions(-) rename {adapter => include/adapter}/adapter.h (99%) rename adapter/time.h => include/adapter/time_handler.h (95%) diff --git a/adapter/adapter.h b/include/adapter/adapter.h similarity index 99% rename from adapter/adapter.h rename to include/adapter/adapter.h index 391d6d9..423ac89 100644 --- a/adapter/adapter.h +++ b/include/adapter/adapter.h @@ -9,10 +9,9 @@ #include #include +#include #include -#include "time.h" - namespace Adapter { using namespace dealii; diff --git a/adapter/time.h b/include/adapter/time_handler.h similarity index 95% rename from adapter/time.h rename to include/adapter/time_handler.h index 80306ba..3e7e486 100644 --- a/adapter/time.h +++ b/include/adapter/time_handler.h @@ -1,5 +1,5 @@ -#ifndef TIME_H -#define TIME_H +#ifndef TIME_HANDLER_H +#define TIME_HANDLER_H #include @@ -81,4 +81,4 @@ namespace Adapter const double delta_t; }; } // namespace Adapter -#endif // TIME_H +#endif // TIME_HANDLER_H From 5e7c9d0db7812edf459f2b34e0d40ffe8edcbb65 Mon Sep 17 00:00:00 2001 From: DavidSCN Date: Sat, 13 Feb 2021 14:08:22 +0100 Subject: [PATCH 02/18] Add a cc file and CMakeLists file --- CMakeLists.txt | 76 ++++++++++++++++++++++++++++++++++++++++ dealii-precice.cc | 88 +++++++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 164 insertions(+) create mode 100644 CMakeLists.txt create mode 100644 dealii-precice.cc diff --git a/CMakeLists.txt b/CMakeLists.txt new file mode 100644 index 0000000..97e278a --- /dev/null +++ b/CMakeLists.txt @@ -0,0 +1,76 @@ +CMAKE_MINIMUM_REQUIRED(VERSION 3.1.0) + +FIND_PACKAGE(deal.II 9.1 QUIET + HINTS ${deal.II_DIR} ${DEAL_II_DIR} ../ ../../ $ENV{DEAL_II_DIR} + ) +IF(NOT ${deal.II_FOUND}) + MESSAGE(FATAL_ERROR "\n" + "*** Could not locate deal.II. ***\n\n" + "You may want to either pass a flag -DDEAL_II_DIR=/path/to/deal.II to cmake\n" + "or set an environment variable \"DEAL_II_DIR\" that contains this path." + ) +ENDIF() + +# Enable a switchable dimension choice +IF (NOT DEFINED DIM) + SET(DIM 2) +ENDIF() +ADD_DEFINITIONS(-DDIM=${DIM}) + +# Set the target and the target source +SET( TARGET "dealii-precice" ) +SET( TARGET_SRC ${TARGET}.cc ) + +DEAL_II_INITIALIZE_CACHED_VARIABLES() + +IF(NOT CMAKE_BUILD_TYPE) + SET(CMAKE_BUILD_TYPE "RELEASE") + MESSAGE(STATUS "No build type specified. Building in ${CMAKE_BUILD_TYPE} mode.") +ENDIF() + +# Set the include directory and the name of the project +PROJECT(${TARGET}) +INCLUDE_DIRECTORIES(${CMAKE_SOURCE_DIR}/include) + +# Set the executable +ADD_EXECUTABLE(${TARGET} ${TARGET_SRC}) + +# Use deal.II macros for setting up the target +DEAL_II_SETUP_TARGET(${TARGET}) +DEAL_II_INITIALIZE_CACHED_VARIABLES() + +# Query the git information and set it in the source +DEAL_II_QUERY_GIT_INFORMATION() +SET_PROPERTY(TARGET ${TARGET} APPEND PROPERTY COMPILE_DEFINITIONS + GIT_BRANCH="${GIT_BRANCH}" + GIT_REVISION="${GIT_REVISION}" + GIT_SHORTREV="${GIT_SHORTREV}") + + +FIND_PACKAGE(precice REQUIRED + HINTS ${PRECICE_DIR} $ENV{precice_DIR}) +TARGET_LINK_LIBRARIES(${TARGET} precice::precice) + +MESSAGE(STATUS "Using the preCICE version found at ${precice_CONFIG}") + +# add the individual solver +ADD_SUBDIRECTORY(source) + +# ...and link them +TARGET_LINK_LIBRARIES(${TARGET} nonlinear_elasticity linear_elasticity) + + +# +# Custom "debug" and "release" make targets: +# +ADD_CUSTOM_TARGET(debug +COMMAND ${CMAKE_COMMAND} -DCMAKE_BUILD_TYPE=Debug ${CMAKE_SOURCE_DIR} +COMMAND ${CMAKE_COMMAND} --build ${CMAKE_BINARY_DIR} --target all +COMMENT "Switch CMAKE_BUILD_TYPE to Debug" +) + +ADD_CUSTOM_TARGET(release +COMMAND ${CMAKE_COMMAND} -DCMAKE_BUILD_TYPE=Release ${CMAKE_SOURCE_DIR} +COMMAND ${CMAKE_COMMAND} --build ${CMAKE_BINARY_DIR} --target all +COMMENT "Switch CMAKE_BUILD_TYPE to Release" +) diff --git a/dealii-precice.cc b/dealii-precice.cc new file mode 100644 index 0000000..15a08d1 --- /dev/null +++ b/dealii-precice.cc @@ -0,0 +1,88 @@ +//#include +#include "source/nonlinear_elasticity/include/nonlinear_elasticity.h" + + +int +main(int argc, char **argv) +{ + using namespace dealii; + +#ifdef DEAL_II_WITH_MPI + Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv); +#endif + + try + { + deallog.depth_console(0); + static const unsigned int n_threads = MultithreadInfo::n_threads(); + + // Query adapter and deal.II info + const std::string adapter_info = + GIT_SHORTREV == std::string("") ? + "unknown" : + (GIT_SHORTREV + std::string(" on branch ") + "GIT_BRANCH"); + const std::string dealii_info = + DEAL_II_GIT_SHORTREV == std::string("") ? + "unknown" : + (DEAL_II_GIT_SHORTREV + std::string(" on branch ") + + DEAL_II_GIT_BRANCH); + + std::cout + << "-----------------------------------------------------------------------------" + << std::endl + << "-- . running with " << n_threads << " thread" + << (n_threads == 1 ? "" : "s") << std::endl; + + std::cout << "-- . adapter revision " << adapter_info << std::endl; + std::cout << "-- . deal.II " << DEAL_II_PACKAGE_VERSION + << " (revision " << dealii_info << ")" << std::endl; + std::cout + << "-----------------------------------------------------------------------------" + << std::endl + << std::endl; + + + std::string parameter_file; + if (argc > 1) + parameter_file = argv[1]; + else + parameter_file = "nonlinear_elasticity.prm"; + + // Extract case path for the output directory + size_t pos = parameter_file.find_last_of("/"); + std::string case_path = + std::string::npos == pos ? "" : parameter_file.substr(0, pos + 1); + + // Dimension is determinded via cmake -DDIM + Nonlinear_Elasticity::Solid solid(case_path); + solid.run(); + } + catch (std::exception &exc) + { + std::cerr << std::endl + << std::endl + << "----------------------------------------------------" + << std::endl; + std::cerr << "Exception on processing: " << std::endl + << exc.what() << std::endl + << "Aborting!" << std::endl + << "----------------------------------------------------" + << std::endl; + + return 1; + } + catch (...) + { + std::cerr << std::endl + << std::endl + << "----------------------------------------------------" + << std::endl; + std::cerr << "Unknown exception!" << std::endl + << "Aborting!" << std::endl + << "----------------------------------------------------" + << std::endl; + return 1; + } + + return 0; +} From 4edff83ab0f2a7ef64a1152716d57b0da0618d5a Mon Sep 17 00:00:00 2001 From: DavidSCN Date: Sat, 13 Feb 2021 20:34:22 +0100 Subject: [PATCH 03/18] Add cmake macros and linear elasticity header --- dealii-precice.cc | 22 ++- source/CMakeLists.txt | 4 + .../include/linear_elasticity.h | 155 ++++++++++++++++++ 3 files changed, 175 insertions(+), 6 deletions(-) create mode 100644 source/CMakeLists.txt create mode 100644 source/linear_elasticity/include/linear_elasticity.h diff --git a/dealii-precice.cc b/dealii-precice.cc index 15a08d1..e9ae8a4 100644 --- a/dealii-precice.cc +++ b/dealii-precice.cc @@ -1,7 +1,6 @@ -//#include +#include "source/linear_elasticity/include/linear_elasticity.h" #include "source/nonlinear_elasticity/include/nonlinear_elasticity.h" - int main(int argc, char **argv) { @@ -20,7 +19,7 @@ main(int argc, char **argv) const std::string adapter_info = GIT_SHORTREV == std::string("") ? "unknown" : - (GIT_SHORTREV + std::string(" on branch ") + "GIT_BRANCH"); + (GIT_SHORTREV + std::string(" on branch ") + GIT_BRANCH); const std::string dealii_info = DEAL_II_GIT_SHORTREV == std::string("") ? "unknown" : @@ -44,7 +43,7 @@ main(int argc, char **argv) std::string parameter_file; if (argc > 1) - parameter_file = argv[1]; + parameter_file = argv[2]; else parameter_file = "nonlinear_elasticity.prm"; @@ -53,9 +52,20 @@ main(int argc, char **argv) std::string case_path = std::string::npos == pos ? "" : parameter_file.substr(0, pos + 1); + std::string solver_type = argv[1]; // Dimension is determinded via cmake -DDIM - Nonlinear_Elasticity::Solid solid(case_path); - solid.run(); + if (solver_type == "-nonlinear") // nonlinear + { + Nonlinear_Elasticity::Solid solid(case_path); + solid.run(); + } + else if (solver_type == "-linear") // linear + { + Linear_Elasticity::ElastoDynamics elastic_solver(case_path); + elastic_solver.run(); + } + else + AssertThrow(false, ExcMessage("Unknown solver specified. ")) } catch (std::exception &exc) { diff --git a/source/CMakeLists.txt b/source/CMakeLists.txt new file mode 100644 index 0000000..d70fb87 --- /dev/null +++ b/source/CMakeLists.txt @@ -0,0 +1,4 @@ +CMAKE_MINIMUM_REQUIRED(VERSION 3.1.0) + +ADD_SUBDIRECTORY(nonlinear_elasticity) +ADD_SUBDIRECTORY(linear_elasticity) diff --git a/source/linear_elasticity/include/linear_elasticity.h b/source/linear_elasticity/include/linear_elasticity.h new file mode 100644 index 0000000..18c2c2e --- /dev/null +++ b/source/linear_elasticity/include/linear_elasticity.h @@ -0,0 +1,155 @@ +#ifndef LINEAR_ELASTICITY_H +#define LINEAR_ELASTICITY_H + +#include +#include +#include +#include +#include +#include + +#include +#include +#include + +#include +#include +#include +#include + +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include + +#include +#include + +#include +#include + +#include "parameter_handling.h" +#include "postprocessor.h" + +// The Linear_Elasticity case includes a linear elastic material with a one-step +// theta time integration +namespace Linear_Elasticity +{ + using namespace dealii; + + template + class ElastoDynamics + { + public: + ElastoDynamics(const std::string &case_path); + + ~ElastoDynamics(); + // As usual in dealii, the run function covers the main time loop of the + // system + void + run(); + + private: + // Create the mesh and set boundary IDs for different boundary conditions + void + make_grid(); + + // Set up the FE system and allocate data structures + void + setup_system(); + + // Compute time invariant matrices e.g. stiffness matrix and mass matrix + void + assemble_system(); + + // Assemble the Neumann contribution i.e. the coupling data obtained from + // the Fluid participant + void + assemble_rhs(); + + void + assemble_consistent_loading(); + + // Solve the linear system + void + solve(); + + // Update the displacement according to the theta scheme + void + update_displacement(); + + // Output results to vtk files + void + output_results() const; + + // Paramter class parsing all user specific input parameters + const Parameters::AllParameters parameters; + + // Boundary IDs, reserved for the respectve application + unsigned int clamped_mesh_id; + unsigned int out_of_plane_clamped_mesh_id; + const unsigned int interface_boundary_id; + + // Dealii typical objects + Triangulation triangulation; + DoFHandler dof_handler; + FESystem fe; + MappingQGeneric mapping; + const unsigned int quad_order; + + AffineConstraints hanging_node_constraints; + + // Matrices used during computations + SparsityPattern sparsity_pattern; + SparseMatrix mass_matrix; + SparseMatrix stiffness_matrix; + SparseMatrix system_matrix; + SparseMatrix stepping_matrix; + + // Time dependent variables + Vector old_velocity; + Vector velocity; + Vector old_displacement; + Vector displacement; + Vector old_stress; + Vector stress; + Vector system_rhs; + + // Body forces e.g. gravity. Values are specified in the input file + const bool body_force_enabled; + Vector body_force_vector; + + // In order to measure some timings + mutable TimerOutput timer; + + // The main adapter objects: The time class keeps track of the current time + // and time steps. The Adapter class includes all functionalities for + // coupling via preCICE. Look at the documentation of the class for more + // information. + Adapter::Time time; + Adapter::Adapter, Parameters::AllParameters> adapter; + + // Alias for all time dependent variables, which should be saved/reloaded + // in case of an implicit coupling. This vector is directly used in the + // Adapter class + std::vector *> state_variables; + // for the output directory + const std::string case_path; + }; +} // namespace Linear_Elasticity + +#endif // LINEAR_ELASTICITY_H From 7e124c65dbede71308cf025b1e4290230354f43f Mon Sep 17 00:00:00 2001 From: DavidSCN Date: Sat, 13 Feb 2021 20:36:28 +0100 Subject: [PATCH 04/18] Port linear solver to new structure --- linear_elasticity/CMakeLists.txt | 68 ---- .../include/parameter_handling.h | 356 ---------------- source/linear_elasticity/CMakeLists.txt | 18 + .../include/parameter_handling.h | 120 ++++++ .../include/postprocessor.h | 4 +- .../linear_elasticity}/linear_elasticity.cc | 383 ++++++++++-------- .../linear_elasticity}/linear_elasticity.prm | 0 7 files changed, 362 insertions(+), 587 deletions(-) delete mode 100644 linear_elasticity/CMakeLists.txt delete mode 100644 linear_elasticity/include/parameter_handling.h create mode 100644 source/linear_elasticity/CMakeLists.txt create mode 100644 source/linear_elasticity/include/parameter_handling.h rename {linear_elasticity => source/linear_elasticity}/include/postprocessor.h (98%) rename {linear_elasticity => source/linear_elasticity}/linear_elasticity.cc (77%) rename {linear_elasticity => source/linear_elasticity}/linear_elasticity.prm (100%) diff --git a/linear_elasticity/CMakeLists.txt b/linear_elasticity/CMakeLists.txt deleted file mode 100644 index c7337e6..0000000 --- a/linear_elasticity/CMakeLists.txt +++ /dev/null @@ -1,68 +0,0 @@ -## -# CMake script for the dealii-adapter: -## - -# Set the name of the project and target: -SET(TARGET "linear_elasticity") - -SET(TARGET_SRC - ${TARGET}.cc - ) - -# To make the dimension switch available via cmake -IF (NOT DEFINED DIM) - SET(DIM 2) -ENDIF() -ADD_DEFINITIONS(-DDIM=${DIM}) - -# Change default build type for tutorial cases -IF(NOT CMAKE_BUILD_TYPE) - SET(CMAKE_BUILD_TYPE "Release" CACHE STRING - "Choose the type of build, options are: Debug Release" - FORCE) -ENDIF(NOT CMAKE_BUILD_TYPE) -MESSAGE(STATUS "Build type: ${CMAKE_BUILD_TYPE}") - -# To cleanup result files -SET(CLEAN_UP_FILES - # a custom list of globs, e.g. *.log *.vtk - *.vtk -) -# Usually, you will not need to modify anything beyond this point... - -CMAKE_MINIMUM_REQUIRED(VERSION 2.8.12) - -FIND_PACKAGE(deal.II 9.2.0 QUIET - HINTS ${deal.II_DIR} ${DEAL_II_DIR} ../ ../../ $ENV{DEAL_II_DIR} - ) -IF(NOT ${deal.II_FOUND}) - MESSAGE(FATAL_ERROR "\n" - "*** Could not locate a (sufficiently recent) version of deal.II. ***\n\n" - "You may want to either pass a flag -DDEAL_II_DIR=/path/to/deal.II to cmake\n" - "or set an environment variable \"DEAL_II_DIR\" that contains this path." - ) -ENDIF() - -DEAL_II_INITIALIZE_CACHED_VARIABLES() - -# Store current source directory -SET(TMP ${CMAKE_SOURCE_DIR}) -SET(CMAKE_SOURCE_DIR "..") - -# Query git information -DEAL_II_QUERY_GIT_INFORMATION() - -# Restore source dir -SET(CMAKE_SOURCE_DIR ${TMP}) - -PROJECT(${TARGET} LANGUAGES CXX) - -DEAL_II_INVOKE_AUTOPILOT() - -FIND_PACKAGE(precice REQUIRED) -TARGET_LINK_LIBRARIES(${TARGET} precice::precice) - -SET_PROPERTY(TARGET ${TARGET} APPEND PROPERTY COMPILE_DEFINITIONS - GIT_BRANCH="${GIT_BRANCH}" - GIT_REVISION="${GIT_REVISION}" - GIT_SHORTREV="${GIT_SHORTREV}") diff --git a/linear_elasticity/include/parameter_handling.h b/linear_elasticity/include/parameter_handling.h deleted file mode 100644 index 49e4b8b..0000000 --- a/linear_elasticity/include/parameter_handling.h +++ /dev/null @@ -1,356 +0,0 @@ -#ifndef parameter_handling_h -#define parameter_handling_h - -#include - -namespace Linear_Elasticity -{ - using namespace dealii; - - /** - * This class declares all parameters, which can be specified in the parameter - * file. The subsection abut preCICE configurations is directly interlinked - * to the Adapter class. - */ - namespace Parameters - { - // TODO: Add more parameters - struct Time - { - double delta_t; - double end_time; - int output_interval; - - static void - declare_parameters(ParameterHandler &prm); - - void - parse_parameters(ParameterHandler &prm); - }; - - void - Time::declare_parameters(ParameterHandler &prm) - { - prm.enter_subsection("Time"); - { - prm.declare_entry("End time", "1", Patterns::Double(), "End time"); - - prm.declare_entry("Time step size", - "0.1", - Patterns::Double(), - "Time step size"); - - prm.declare_entry("Output interval", - "1", - Patterns::Integer(0), - "Write results every x timesteps"); - } - prm.leave_subsection(); - } - - void - Time::parse_parameters(ParameterHandler &prm) - { - prm.enter_subsection("Time"); - { - end_time = prm.get_double("End time"); - delta_t = prm.get_double("Time step size"); - output_interval = prm.get_integer("Output interval"); - } - prm.leave_subsection(); - } - - /** - * @brief Discretization: Specifies parameters for time integration by a - * theta scheme and polynomial degree of the FE system - */ - struct Discretization - { - double theta; - unsigned int poly_degree; - - static void - declare_parameters(ParameterHandler &prm); - - void - parse_parameters(ParameterHandler &prm); - }; - - void - Discretization::declare_parameters(ParameterHandler &prm) - { - prm.enter_subsection("Discretization"); - { - prm.declare_entry("theta", - "0.5", - Patterns::Double(0, 1), - "Time integration scheme"); - - prm.declare_entry("Polynomial degree", - "3", - Patterns::Integer(0), - "Polynomial degree of the FE system"); - } - prm.leave_subsection(); - } - - void - Discretization::parse_parameters(ParameterHandler &prm) - { - prm.enter_subsection("Discretization"); - { - theta = prm.get_double("theta"); - poly_degree = prm.get_integer("Polynomial degree"); - } - prm.leave_subsection(); - } - - /** - * @brief The System struct keeps material properties and body force contributions - */ - struct System - { - double mu; - double lambda; - double rho; - Tensor<1, 3, double> body_force; - - - static void - declare_parameters(ParameterHandler &prm); - - void - parse_parameters(ParameterHandler &prm); - }; - - void - System::declare_parameters(ParameterHandler &prm) - { - prm.enter_subsection("System properties"); - { - prm.declare_entry("mu", "0.5e6", Patterns::Double(), "mu"); - - prm.declare_entry("lambda", "2e6", Patterns::Double(), "lambda"); - - prm.declare_entry("rho", "1000", Patterns::Double(0.0), "density"); - - prm.declare_entry("body forces", - "0,0,0", - Patterns::List(Patterns::Double()), - "body forces x,y,z"); - } - prm.leave_subsection(); - } - - void - System::parse_parameters(ParameterHandler &prm) - { - prm.enter_subsection("System properties"); - { - mu = prm.get_double("mu"); - lambda = prm.get_double("lambda"); - rho = prm.get_double("rho"); - const std::vector body_forces_input = - Utilities::split_string_list(prm.get("body forces")); - for (uint d = 0; d < 3; ++d) - body_force[d] = Utilities::string_to_double(body_forces_input[d]); - } - prm.leave_subsection(); - } - - /** - * @brief LinearSolver: Specifies linear solver properties - */ - struct LinearSolver - { - std::string type_lin; - double tol_lin; - double max_iterations_lin; - - static void - declare_parameters(ParameterHandler &prm); - - void - parse_parameters(ParameterHandler &prm); - }; - - void - LinearSolver::declare_parameters(ParameterHandler &prm) - { - prm.enter_subsection("Linear solver"); - { - prm.declare_entry("Solver type", - "Direct", - Patterns::Selection("CG|Direct"), - "Linear solver: CG or Direct"); - - prm.declare_entry("Residual", - "1e-6", - Patterns::Double(0.0), - "Linear solver residual (scaled by residual norm)"); - - prm.declare_entry( - "Max iteration multiplier", - "1", - Patterns::Double(0.0), - "Linear solver iterations (multiples of the system matrix size)"); - } - prm.leave_subsection(); - } - - void - LinearSolver::parse_parameters(ParameterHandler &prm) - { - prm.enter_subsection("Linear solver"); - { - type_lin = prm.get("Solver type"); - tol_lin = prm.get_double("Residual"); - max_iterations_lin = prm.get_double("Max iteration multiplier"); - } - prm.leave_subsection(); - } - - /** - * @brief PreciceAdapterConfiguration: Specifies preCICE related information. - * A lot of these information need to be consistent with the - * precice-config.xml file. - */ - struct PreciceAdapterConfiguration - { - std::string scenario; - std::string config_file; - std::string participant_name; - std::string mesh_name; - std::string read_data_name; - std::string write_data_name; - double flap_location; - bool data_consistent; - - static void - declare_parameters(ParameterHandler &prm); - - void - parse_parameters(ParameterHandler &prm); - }; - - - void - PreciceAdapterConfiguration::declare_parameters(ParameterHandler &prm) - { - prm.enter_subsection("precice configuration"); - { - prm.declare_entry("Scenario", - "FSI3", - Patterns::Selection("FSI3|PF"), - "Cases: FSI3 or PF for perpendicular flap"); - prm.declare_entry("precice config-file", - "precice-config.xml", - Patterns::Anything(), - "Name of the precice configuration file"); - prm.declare_entry( - "Participant name", - "dealiisolver", - Patterns::Anything(), - "Name of the participant in the precice-config.xml file"); - prm.declare_entry( - "Mesh name", - "dealii-mesh", - Patterns::Anything(), - "Name of the coupling mesh in the precice-config.xml file"); - prm.declare_entry( - "Read data name", - "received-data", - Patterns::Anything(), - "Name of the read data in the precice-config.xml file"); - prm.declare_entry( - "Write data name", - "calculated-data", - Patterns::Anything(), - "Name of the write data in the precice-config.xml file"); - prm.declare_entry("Flap location", - "0.0", - Patterns::Double(-3, 3), - "PF x-location"); - } - prm.leave_subsection(); - } - - void - PreciceAdapterConfiguration::parse_parameters(ParameterHandler &prm) - { - prm.enter_subsection("precice configuration"); - { - scenario = prm.get("Scenario"); - config_file = prm.get("precice config-file"); - participant_name = prm.get("Participant name"); - mesh_name = prm.get("Mesh name"); - read_data_name = prm.get("Read data name"); - write_data_name = prm.get("Write data name"); - flap_location = prm.get_double("Flap location"); - } - prm.leave_subsection(); - // Look at the specific type of read data - if ((read_data_name.find("Stress") == 0)) - data_consistent = true; - else if ((read_data_name.find("Force") == 0)) - data_consistent = false; - else - AssertThrow( - false, - ExcMessage( - "Unknown read data type. Please use 'Force' or 'Stress' in the read data naming.")); - } - - - - struct AllParameters : public LinearSolver, - public Discretization, - public System, - public Time, - public PreciceAdapterConfiguration - - { - AllParameters(const std::string &input_file); - - static void - declare_parameters(ParameterHandler &prm); - - void - parse_parameters(ParameterHandler &prm); - }; - - AllParameters::AllParameters(const std::string &input_file) - { - ParameterHandler prm; - declare_parameters(prm); - prm.parse_input(input_file); - parse_parameters(prm); - - // Optional, if we want to print all parameters in the beginning of the - // simulation - // prm.print_parameters(std::cout,ParameterHandler::Text); - } - - void - AllParameters::declare_parameters(ParameterHandler &prm) - { - LinearSolver::declare_parameters(prm); - Discretization::declare_parameters(prm); - System::declare_parameters(prm); - Time::declare_parameters(prm); - PreciceAdapterConfiguration::declare_parameters(prm); - } - - void - AllParameters::parse_parameters(ParameterHandler &prm) - { - LinearSolver::parse_parameters(prm); - Discretization::parse_parameters(prm); - System::parse_parameters(prm); - Time::parse_parameters(prm); - PreciceAdapterConfiguration::parse_parameters(prm); - } - } // namespace Parameters -} // namespace Linear_Elasticity - -#endif diff --git a/source/linear_elasticity/CMakeLists.txt b/source/linear_elasticity/CMakeLists.txt new file mode 100644 index 0000000..b5c440b --- /dev/null +++ b/source/linear_elasticity/CMakeLists.txt @@ -0,0 +1,18 @@ +## +# CMake script for the dealii-adapter: +## + +# Set the name of the project and target: +SET(TARGET "linear_elasticity") + +SET(_src + ${TARGET}.cc +) + +CMAKE_MINIMUM_REQUIRED(VERSION 3.1.0) + +ADD_LIBRARY(${TARGET} OBJECT ${_src}) +FIND_PACKAGE(precice REQUIRED) +TARGET_LINK_LIBRARIES(${TARGET} precice::precice) + +DEAL_II_SETUP_TARGET(${TARGET}) diff --git a/source/linear_elasticity/include/parameter_handling.h b/source/linear_elasticity/include/parameter_handling.h new file mode 100644 index 0000000..3a22274 --- /dev/null +++ b/source/linear_elasticity/include/parameter_handling.h @@ -0,0 +1,120 @@ +#pragma once + +#include + +namespace Linear_Elasticity +{ + using namespace dealii; + + /** + * This class declares all parameters, which can be specified in the parameter + * file. The subsection abut preCICE configurations is directly interlinked + * to the Adapter class. + */ + namespace Parameters + { + // TODO: Add more parameters + struct Time + { + double delta_t; + double end_time; + int output_interval; + + static void + declare_parameters(ParameterHandler &prm); + + void + parse_parameters(ParameterHandler &prm); + }; + + /** + * @brief Discretization: Specifies parameters for time integration by a + * theta scheme and polynomial degree of the FE system + */ + struct Discretization + { + double theta; + unsigned int poly_degree; + + static void + declare_parameters(ParameterHandler &prm); + + void + parse_parameters(ParameterHandler &prm); + }; + + /** + * @brief The System struct keeps material properties and body force contributions + */ + struct System + { + double mu; + double lambda; + double rho; + Tensor<1, 3, double> body_force; + + + static void + declare_parameters(ParameterHandler &prm); + + void + parse_parameters(ParameterHandler &prm); + }; + + /** + * @brief LinearSolver: Specifies linear solver properties + */ + struct LinearSolver + { + std::string type_lin; + double tol_lin; + double max_iterations_lin; + + static void + declare_parameters(ParameterHandler &prm); + + void + parse_parameters(ParameterHandler &prm); + }; + + /** + * @brief PreciceAdapterConfiguration: Specifies preCICE related information. + * A lot of these information need to be consistent with the + * precice-config.xml file. + */ + struct PreciceAdapterConfiguration + { + std::string scenario; + std::string config_file; + std::string participant_name; + std::string mesh_name; + std::string read_data_name; + std::string write_data_name; + double flap_location; + bool data_consistent; + + static void + declare_parameters(ParameterHandler &prm); + + void + parse_parameters(ParameterHandler &prm); + }; + + + struct AllParameters : public LinearSolver, + public Discretization, + public System, + public Time, + public PreciceAdapterConfiguration + + { + AllParameters(const std::string &input_file); + + static void + declare_parameters(ParameterHandler &prm); + + void + parse_parameters(ParameterHandler &prm); + }; + } // namespace Parameters +} // namespace Linear_Elasticity diff --git a/linear_elasticity/include/postprocessor.h b/source/linear_elasticity/include/postprocessor.h similarity index 98% rename from linear_elasticity/include/postprocessor.h rename to source/linear_elasticity/include/postprocessor.h index c3abd58..a36b024 100644 --- a/linear_elasticity/include/postprocessor.h +++ b/source/linear_elasticity/include/postprocessor.h @@ -1,5 +1,4 @@ -#ifndef POSTPROCESSOR_H -#define POSTPROCESSOR_H +#pragma once #include @@ -122,4 +121,3 @@ namespace Linear_Elasticity } } // namespace Linear_Elasticity -#endif // POSTPROCESSOR_H diff --git a/linear_elasticity/linear_elasticity.cc b/source/linear_elasticity/linear_elasticity.cc similarity index 77% rename from linear_elasticity/linear_elasticity.cc rename to source/linear_elasticity/linear_elasticity.cc index 6eac8e9..6dc0ef3 100644 --- a/linear_elasticity/linear_elasticity.cc +++ b/source/linear_elasticity/linear_elasticity.cc @@ -1,3 +1,5 @@ +#include "include/linear_elasticity.h" + #include #include #include @@ -33,11 +35,12 @@ #include #include +#include +#include + #include #include -#include "../adapter/adapter.h" -#include "../adapter/time.h" #include "include/parameter_handling.h" #include "include/postprocessor.h" @@ -47,105 +50,241 @@ namespace Linear_Elasticity { using namespace dealii; - template - class ElastoDynamics + namespace Parameters { - public: - ElastoDynamics(const std::string &case_path); - ~ElastoDynamics(); - // As usual in dealii, the run function covers the main time loop of the - // system void - run(); + Time::declare_parameters(ParameterHandler &prm) + { + prm.enter_subsection("Time"); + { + prm.declare_entry("End time", "1", Patterns::Double(), "End time"); + + prm.declare_entry("Time step size", + "0.1", + Patterns::Double(), + "Time step size"); + + prm.declare_entry("Output interval", + "1", + Patterns::Integer(0), + "Write results every x timesteps"); + } + prm.leave_subsection(); + } - private: - // Create the mesh and set boundary IDs for different boundary conditions void - make_grid(); + Time::parse_parameters(ParameterHandler &prm) + { + prm.enter_subsection("Time"); + { + end_time = prm.get_double("End time"); + delta_t = prm.get_double("Time step size"); + output_interval = prm.get_integer("Output interval"); + } + prm.leave_subsection(); + } - // Set up the FE system and allocate data structures void - setup_system(); + Discretization::declare_parameters(ParameterHandler &prm) + { + prm.enter_subsection("Discretization"); + { + prm.declare_entry("theta", + "0.5", + Patterns::Double(0, 1), + "Time integration scheme"); + + prm.declare_entry("Polynomial degree", + "3", + Patterns::Integer(0), + "Polynomial degree of the FE system"); + } + prm.leave_subsection(); + } - // Compute time invariant matrices e.g. stiffness matrix and mass matrix void - assemble_system(); + Discretization::parse_parameters(ParameterHandler &prm) + { + prm.enter_subsection("Discretization"); + { + theta = prm.get_double("theta"); + poly_degree = prm.get_integer("Polynomial degree"); + } + prm.leave_subsection(); + } + + void + System::declare_parameters(ParameterHandler &prm) + { + prm.enter_subsection("System properties"); + { + prm.declare_entry("mu", "0.5e6", Patterns::Double(), "mu"); + + prm.declare_entry("lambda", "2e6", Patterns::Double(), "lambda"); + + prm.declare_entry("rho", "1000", Patterns::Double(0.0), "density"); + + prm.declare_entry("body forces", + "0,0,0", + Patterns::List(Patterns::Double()), + "body forces x,y,z"); + } + prm.leave_subsection(); + } + + void + System::parse_parameters(ParameterHandler &prm) + { + prm.enter_subsection("System properties"); + { + mu = prm.get_double("mu"); + lambda = prm.get_double("lambda"); + rho = prm.get_double("rho"); + const std::vector body_forces_input = + Utilities::split_string_list(prm.get("body forces")); + for (uint d = 0; d < 3; ++d) + body_force[d] = Utilities::string_to_double(body_forces_input[d]); + } + prm.leave_subsection(); + } + + + void + LinearSolver::declare_parameters(ParameterHandler &prm) + { + prm.enter_subsection("Linear solver"); + { + prm.declare_entry("Solver type", + "Direct", + Patterns::Selection("CG|Direct"), + "Linear solver: CG or Direct"); + + prm.declare_entry("Residual", + "1e-6", + Patterns::Double(0.0), + "Linear solver residual (scaled by residual norm)"); + + prm.declare_entry( + "Max iteration multiplier", + "1", + Patterns::Double(0.0), + "Linear solver iterations (multiples of the system matrix size)"); + } + prm.leave_subsection(); + } - // Assemble the Neumann contribution i.e. the coupling data obtained from - // the Fluid participant void - assemble_rhs(); + LinearSolver::parse_parameters(ParameterHandler &prm) + { + prm.enter_subsection("Linear solver"); + { + type_lin = prm.get("Solver type"); + tol_lin = prm.get_double("Residual"); + max_iterations_lin = prm.get_double("Max iteration multiplier"); + } + prm.leave_subsection(); + } void - assemble_consistent_loading(); + PreciceAdapterConfiguration::declare_parameters(ParameterHandler &prm) + { + prm.enter_subsection("precice configuration"); + { + prm.declare_entry("Scenario", + "FSI3", + Patterns::Selection("FSI3|PF"), + "Cases: FSI3 or PF for perpendicular flap"); + prm.declare_entry("precice config-file", + "precice-config.xml", + Patterns::Anything(), + "Name of the precice configuration file"); + prm.declare_entry( + "Participant name", + "dealiisolver", + Patterns::Anything(), + "Name of the participant in the precice-config.xml file"); + prm.declare_entry( + "Mesh name", + "dealii-mesh", + Patterns::Anything(), + "Name of the coupling mesh in the precice-config.xml file"); + prm.declare_entry( + "Read data name", + "received-data", + Patterns::Anything(), + "Name of the read data in the precice-config.xml file"); + prm.declare_entry( + "Write data name", + "calculated-data", + Patterns::Anything(), + "Name of the write data in the precice-config.xml file"); + prm.declare_entry("Flap location", + "0.0", + Patterns::Double(-3, 3), + "PF x-location"); + } + prm.leave_subsection(); + } - // Solve the linear system void - solve(); + PreciceAdapterConfiguration::parse_parameters(ParameterHandler &prm) + { + prm.enter_subsection("precice configuration"); + { + scenario = prm.get("Scenario"); + config_file = prm.get("precice config-file"); + participant_name = prm.get("Participant name"); + mesh_name = prm.get("Mesh name"); + read_data_name = prm.get("Read data name"); + write_data_name = prm.get("Write data name"); + flap_location = prm.get_double("Flap location"); + } + prm.leave_subsection(); + // Look at the specific type of read data + if ((read_data_name.find("Stress") == 0)) + data_consistent = true; + else if ((read_data_name.find("Force") == 0)) + data_consistent = false; + else + AssertThrow( + false, + ExcMessage( + "Unknown read data type. Please use 'Force' or 'Stress' in the read data naming.")); + } + + AllParameters::AllParameters(const std::string &input_file) + { + ParameterHandler prm; + declare_parameters(prm); + prm.parse_input(input_file); + parse_parameters(prm); + + // Optional, if we want to print all parameters in the beginning of the + // simulation + // prm.print_parameters(std::cout,ParameterHandler::Text); + } - // Update the displacement according to the theta scheme void - update_displacement(); + AllParameters::declare_parameters(ParameterHandler &prm) + { + LinearSolver::declare_parameters(prm); + Discretization::declare_parameters(prm); + System::declare_parameters(prm); + Time::declare_parameters(prm); + PreciceAdapterConfiguration::declare_parameters(prm); + } - // Output results to vtk files void - output_results() const; - - // Paramter class parsing all user specific input parameters - const Parameters::AllParameters parameters; - - // Boundary IDs, reserved for the respectve application - unsigned int clamped_mesh_id; - unsigned int out_of_plane_clamped_mesh_id; - const unsigned int interface_boundary_id; - - // Dealii typical objects - Triangulation triangulation; - DoFHandler dof_handler; - FESystem fe; - MappingQGeneric mapping; - const unsigned int quad_order; - - AffineConstraints hanging_node_constraints; - - // Matrices used during computations - SparsityPattern sparsity_pattern; - SparseMatrix mass_matrix; - SparseMatrix stiffness_matrix; - SparseMatrix system_matrix; - SparseMatrix stepping_matrix; - - // Time dependent variables - Vector old_velocity; - Vector velocity; - Vector old_displacement; - Vector displacement; - Vector old_stress; - Vector stress; - Vector system_rhs; - - // Body forces e.g. gravity. Values are specified in the input file - const bool body_force_enabled; - Vector body_force_vector; - - // In order to measure some timings - mutable TimerOutput timer; - - // The main adapter objects: The time class keeps track of the current time - // and time steps. The Adapter class includes all functionalities for - // coupling via preCICE. Look at the documentation of the class for more - // information. - Adapter::Time time; - Adapter::Adapter, Parameters::AllParameters> adapter; - - // Alias for all time dependent variables, which should be saved/reloaded - // in case of an implicit coupling. This vector is directly used in the - // Adapter class - std::vector *> state_variables; - // for the output directory - const std::string case_path; - }; + AllParameters::parse_parameters(ParameterHandler &prm) + { + LinearSolver::parse_parameters(prm); + Discretization::parse_parameters(prm); + System::parse_parameters(prm); + Time::parse_parameters(prm); + PreciceAdapterConfiguration::parse_parameters(prm); + } + } // namespace Parameters // Constructor @@ -807,82 +946,6 @@ namespace Linear_Elasticity // communication etc. adapter.precice.finalize(); } -} // namespace Linear_Elasticity - -int -main(int argc, char **argv) -{ - using namespace Linear_Elasticity; - using namespace dealii; - -#ifdef DEAL_II_WITH_MPI - Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv); -#endif - - try - { - // Query adapter and deal.II info - const std::string adapter_info = - GIT_SHORTREV == std::string("") ? - "unknown" : - (GIT_SHORTREV + std::string(" on branch ") + GIT_BRANCH); - const std::string dealii_info = - DEAL_II_GIT_SHORTREV == std::string("") ? - "unknown" : - (DEAL_II_GIT_SHORTREV + std::string(" on branch ") + - DEAL_II_GIT_BRANCH); - - std::cout - << "-----------------------------------------------------------------------------" - << std::endl; - std::cout << "-- . adapter revision " << adapter_info << std::endl; - std::cout << "-- . deal.II " << DEAL_II_PACKAGE_VERSION - << " (revision " << dealii_info << ")" << std::endl; - std::cout - << "-----------------------------------------------------------------------------" - << std::endl - << std::endl; - - std::string parameter_file; - if (argc > 1) - parameter_file = argv[1]; - else - parameter_file = "linear_elasticity.prm"; - - // Extract case path for the output directory - size_t pos = parameter_file.find_last_of("/"); - std::string case_path = - std::string::npos == pos ? "" : parameter_file.substr(0, pos + 1); - - ElastoDynamics elastic_solver(case_path); - elastic_solver.run(); - } - catch (std::exception &exc) - { - std::cerr << std::endl - << std::endl - << "----------------------------------------------------" - << std::endl; - std::cerr << "Exception on processing: " << std::endl - << exc.what() << std::endl - << "Aborting!" << std::endl - << "----------------------------------------------------" - << std::endl; - - return 1; - } - catch (...) - { - std::cerr << std::endl - << std::endl - << "----------------------------------------------------" - << std::endl; - std::cerr << "Unknown exception!" << std::endl - << "Aborting!" << std::endl - << "----------------------------------------------------" - << std::endl; - return 1; - } - return 0; -} + template class ElastoDynamics; +} // namespace Linear_Elasticity diff --git a/linear_elasticity/linear_elasticity.prm b/source/linear_elasticity/linear_elasticity.prm similarity index 100% rename from linear_elasticity/linear_elasticity.prm rename to source/linear_elasticity/linear_elasticity.prm From 84729672bdaccdaac1dbb1f98c61e591a6f586f8 Mon Sep 17 00:00:00 2001 From: DavidSCN Date: Sat, 13 Feb 2021 20:37:57 +0100 Subject: [PATCH 05/18] port nonlinear case to new structure --- nonlinear_elasticity/CMakeLists.txt | 68 --- .../include/parameter_handling.h | 414 ------------- source/nonlinear_elasticity/CMakeLists.txt | 18 + .../include/compressible_neo_hook_material.h | 4 +- .../include/nonlinear_elasticity.h | 338 +++++++++++ .../include/parameter_handling.h | 145 +++++ .../include/postprocessor.h | 4 +- .../nonlinear_elasticity.cc | 565 ++++++++---------- .../nonlinear_elasticity.prm | 0 9 files changed, 742 insertions(+), 814 deletions(-) delete mode 100644 nonlinear_elasticity/CMakeLists.txt delete mode 100644 nonlinear_elasticity/include/parameter_handling.h create mode 100644 source/nonlinear_elasticity/CMakeLists.txt rename {nonlinear_elasticity => source/nonlinear_elasticity}/include/compressible_neo_hook_material.h (97%) create mode 100644 source/nonlinear_elasticity/include/nonlinear_elasticity.h create mode 100644 source/nonlinear_elasticity/include/parameter_handling.h rename {nonlinear_elasticity => source/nonlinear_elasticity}/include/postprocessor.h (98%) rename {nonlinear_elasticity => source/nonlinear_elasticity}/nonlinear_elasticity.cc (80%) rename {nonlinear_elasticity => source/nonlinear_elasticity}/nonlinear_elasticity.prm (100%) diff --git a/nonlinear_elasticity/CMakeLists.txt b/nonlinear_elasticity/CMakeLists.txt deleted file mode 100644 index b72a6fc..0000000 --- a/nonlinear_elasticity/CMakeLists.txt +++ /dev/null @@ -1,68 +0,0 @@ -## -# CMake script for the dealii-adapter: -## - -# Set the name of the project and target: -SET(TARGET "nonlinear_elasticity") - -SET(TARGET_SRC - ${TARGET}.cc - ) - -# To make the dimension switch available via cmake -IF (NOT DEFINED DIM) - SET(DIM 2) -ENDIF() -ADD_DEFINITIONS(-DDIM=${DIM}) - -# Change default build type for tutorial cases -IF(NOT CMAKE_BUILD_TYPE) - SET(CMAKE_BUILD_TYPE "Release" CACHE STRING - "Choose the type of build, options are: Debug Release" - FORCE) -ENDIF(NOT CMAKE_BUILD_TYPE) -MESSAGE(STATUS "Build type: ${CMAKE_BUILD_TYPE}") - -# To cleanup result files -SET(CLEAN_UP_FILES - # a custom list of globs, e.g. *.log *.vtk - *.vtk -) -# Usually, you will not need to modify anything beyond this point... - -CMAKE_MINIMUM_REQUIRED(VERSION 2.8.12) - -FIND_PACKAGE(deal.II 9.2.0 QUIET - HINTS ${deal.II_DIR} ${DEAL_II_DIR} ../ ../../ $ENV{DEAL_II_DIR} - ) -IF(NOT ${deal.II_FOUND}) - MESSAGE(FATAL_ERROR "\n" - "*** Could not locate a (sufficiently recent) version of deal.II. ***\n\n" - "You may want to either pass a flag -DDEAL_II_DIR=/path/to/deal.II to cmake\n" - "or set an environment variable \"DEAL_II_DIR\" that contains this path." - ) -ENDIF() - -DEAL_II_INITIALIZE_CACHED_VARIABLES() - -# Store current source directory -SET(TMP ${CMAKE_SOURCE_DIR}) -SET(CMAKE_SOURCE_DIR "..") - -# Query git information -DEAL_II_QUERY_GIT_INFORMATION() - -# Restore source dir -SET(CMAKE_SOURCE_DIR ${TMP}) - -PROJECT(${TARGET} LANGUAGES CXX) - -DEAL_II_INVOKE_AUTOPILOT() - -FIND_PACKAGE(precice REQUIRED) -TARGET_LINK_LIBRARIES(${TARGET} precice::precice) - -SET_PROPERTY(TARGET ${TARGET} APPEND PROPERTY COMPILE_DEFINITIONS - GIT_BRANCH="${GIT_BRANCH}" - GIT_REVISION="${GIT_REVISION}" - GIT_SHORTREV="${GIT_SHORTREV}") diff --git a/nonlinear_elasticity/include/parameter_handling.h b/nonlinear_elasticity/include/parameter_handling.h deleted file mode 100644 index bdb5952..0000000 --- a/nonlinear_elasticity/include/parameter_handling.h +++ /dev/null @@ -1,414 +0,0 @@ -#ifndef parameter_handling_h -#define parameter_handling_h - -#include - -namespace Nonlinear_Elasticity -{ - using namespace dealii; - - /** - * This class declares all parameters, which can be specified in the parameter - * file. The subsection abut preCICE configurations is directly interlinked - * to the Adapter class. - */ - namespace Parameters - { - /** - * @brief System: Specifies system properties - */ - struct System - { - double nu; - double mu; - double rho; - Tensor<1, 3, double> body_force; - - static void - declare_parameters(ParameterHandler &prm); - - void - parse_parameters(ParameterHandler &prm); - }; - - void - System::declare_parameters(ParameterHandler &prm) - { - prm.enter_subsection("System properties"); - { - prm.declare_entry("Poisson's ratio", - "0.3", - Patterns::Double(-1.0, 0.5), - "Poisson's ratio"); - - prm.declare_entry("Shear modulus", - "0.4225e6", - Patterns::Double(), - "Shear modulus"); - - prm.declare_entry("rho", "1000", Patterns::Double(), "rho"); - - prm.declare_entry("body forces", - "0,0,0", - Patterns::List(Patterns::Double()), - "body forces x,y,z"); - } - prm.leave_subsection(); - } - - void - System::parse_parameters(ParameterHandler &prm) - { - prm.enter_subsection("System properties"); - { - nu = prm.get_double("Poisson's ratio"); - mu = prm.get_double("Shear modulus"); - rho = prm.get_double("rho"); - const std::vector body_forces_input = - Utilities::split_string_list(prm.get("body forces")); - for (uint d = 0; d < 3; ++d) - body_force[d] = Utilities::string_to_double(body_forces_input[d]); - } - prm.leave_subsection(); - } - - /** - * @brief LinearSolver: Specifies linear solver properties - */ - struct LinearSolver - { - std::string type_lin; - double tol_lin; - double max_iterations_lin; - - static void - declare_parameters(ParameterHandler &prm); - - void - parse_parameters(ParameterHandler &prm); - }; - - void - LinearSolver::declare_parameters(ParameterHandler &prm) - { - prm.enter_subsection("Linear solver"); - { - prm.declare_entry("Solver type", - "CG", - Patterns::Selection("CG|Direct"), - "Linear solver: CG or Direct"); - - prm.declare_entry("Residual", - "1e-6", - Patterns::Double(0.0), - "Linear solver residual (scaled by residual norm)"); - - prm.declare_entry( - "Max iteration multiplier", - "1", - Patterns::Double(0.0), - "Linear solver iterations (multiples of the system matrix size)"); - } - prm.leave_subsection(); - } - - void - LinearSolver::parse_parameters(ParameterHandler &prm) - { - prm.enter_subsection("Linear solver"); - { - type_lin = prm.get("Solver type"); - tol_lin = prm.get_double("Residual"); - max_iterations_lin = prm.get_double("Max iteration multiplier"); - } - prm.leave_subsection(); - } - - /** - * @brief NonlinearSolver: Specifies nonlinear solver properties - */ - struct NonlinearSolver - { - unsigned int max_iterations_NR; - double tol_f; - double tol_u; - - static void - declare_parameters(ParameterHandler &prm); - - void - parse_parameters(ParameterHandler &prm); - }; - - void - NonlinearSolver::declare_parameters(ParameterHandler &prm) - { - prm.enter_subsection("Nonlinear solver"); - { - prm.declare_entry("Max iterations Newton-Raphson", - "10", - Patterns::Integer(0), - "Number of Newton-Raphson iterations allowed"); - - prm.declare_entry("Tolerance force", - "1.0e-9", - Patterns::Double(0.0), - "Force residual tolerance"); - - prm.declare_entry("Tolerance displacement", - "1.0e-6", - Patterns::Double(0.0), - "Displacement error tolerance"); - } - prm.leave_subsection(); - } - - void - NonlinearSolver::parse_parameters(ParameterHandler &prm) - { - prm.enter_subsection("Nonlinear solver"); - { - max_iterations_NR = prm.get_integer("Max iterations Newton-Raphson"); - tol_f = prm.get_double("Tolerance force"); - tol_u = prm.get_double("Tolerance displacement"); - } - prm.leave_subsection(); - } - - /** - * @brief Time: Specifies simulation time properties - */ - struct Time - { - double delta_t; - double end_time; - int output_interval; - - static void - declare_parameters(ParameterHandler &prm); - - void - parse_parameters(ParameterHandler &prm); - }; - - void - Time::declare_parameters(ParameterHandler &prm) - { - prm.enter_subsection("Time"); - { - prm.declare_entry("End time", "1", Patterns::Double(), "End time"); - - prm.declare_entry("Time step size", - "0.1", - Patterns::Double(), - "Time step size"); - - prm.declare_entry("Output interval", - "1", - Patterns::Integer(), - "Output interval"); - } - prm.leave_subsection(); - } - - void - Time::parse_parameters(ParameterHandler &prm) - { - prm.enter_subsection("Time"); - { - end_time = prm.get_double("End time"); - delta_t = prm.get_double("Time step size"); - output_interval = prm.get_integer("Output interval"); - } - prm.leave_subsection(); - } - - /** - * @brief Discretization: Specifies parameters for time integration by an - * implicit Newmark scheme and polynomial degree of the FE system - */ - struct Discretization - { - double beta; - double gamma; - unsigned int poly_degree; - - - static void - declare_parameters(ParameterHandler &prm); - - void - parse_parameters(ParameterHandler &prm); - }; - - void - Discretization::declare_parameters(ParameterHandler &prm) - { - prm.enter_subsection("Discretization"); - { - prm.declare_entry("beta", - "0.25", - Patterns::Double(0, 0.5), - "Newmark beta"); - - prm.declare_entry("gamma", - "0.5", - Patterns::Double(0, 1), - "Newmark gamma"); - - prm.declare_entry("Polynomial degree", - "3", - Patterns::Integer(0), - "Polynomial degree of the FE system"); - } - prm.leave_subsection(); - } - - void - Discretization::parse_parameters(ParameterHandler &prm) - { - prm.enter_subsection("Discretization"); - { - beta = prm.get_double("beta"); - gamma = prm.get_double("gamma"); - poly_degree = prm.get_integer("Polynomial degree"); - } - prm.leave_subsection(); - } - - /** - * @brief PreciceAdapterConfiguration: Specifies preCICE related information. - * A lot of these information need to be consistent with the - * precice-config.xml file. - */ - struct PreciceAdapterConfiguration - { - std::string scenario; - std::string config_file; - std::string participant_name; - std::string mesh_name; - std::string read_data_name; - std::string write_data_name; - double flap_location; - - static void - declare_parameters(ParameterHandler &prm); - - void - parse_parameters(ParameterHandler &prm); - }; - - - void - PreciceAdapterConfiguration::declare_parameters(ParameterHandler &prm) - { - prm.enter_subsection("precice configuration"); - { - prm.declare_entry("Scenario", - "FSI3", - Patterns::Selection("FSI3|PF"), - "Cases: FSI3 or PF for perpendicular flap"); - prm.declare_entry("precice config-file", - "precice-config.xml", - Patterns::Anything(), - "Name of the precice configuration file"); - prm.declare_entry( - "Participant name", - "dealiisolver", - Patterns::Anything(), - "Name of the participant in the precice-config.xml file"); - prm.declare_entry( - "Mesh name", - "dealii-mesh", - Patterns::Anything(), - "Name of the coupling mesh in the precice-config.xml file"); - prm.declare_entry( - "Read data name", - "received-data", - Patterns::Anything(), - "Name of the read data in the precice-config.xml file"); - prm.declare_entry( - "Write data name", - "calculated-data", - Patterns::Anything(), - "Name of the write data in the precice-config.xml file"); - prm.declare_entry("Flap location", - "0.0", - Patterns::Double(-3, 3), - "PF x-location"); - } - prm.leave_subsection(); - } - - void - PreciceAdapterConfiguration::parse_parameters(ParameterHandler &prm) - { - prm.enter_subsection("precice configuration"); - { - scenario = prm.get("Scenario"); - config_file = prm.get("precice config-file"); - participant_name = prm.get("Participant name"); - mesh_name = prm.get("Mesh name"); - read_data_name = prm.get("Read data name"); - write_data_name = prm.get("Write data name"); - flap_location = prm.get_double("Flap location"); - } - prm.leave_subsection(); - } - - - - struct AllParameters : public System, - public LinearSolver, - public NonlinearSolver, - public Time, - public Discretization, - public PreciceAdapterConfiguration - - { - AllParameters(const std::string &input_file); - - static void - declare_parameters(ParameterHandler &prm); - - void - parse_parameters(ParameterHandler &prm); - }; - - AllParameters::AllParameters(const std::string &input_file) - { - ParameterHandler prm; - declare_parameters(prm); - prm.parse_input(input_file); - parse_parameters(prm); - - // Optional, if we want to print all parameters in the beginning of the - // simulation - // prm.print_parameters(std::cout,ParameterHandler::Text); - } - - void - AllParameters::declare_parameters(ParameterHandler &prm) - { - System::declare_parameters(prm); - LinearSolver::declare_parameters(prm); - NonlinearSolver::declare_parameters(prm); - Time::declare_parameters(prm); - Discretization::declare_parameters(prm); - PreciceAdapterConfiguration::declare_parameters(prm); - } - - void - AllParameters::parse_parameters(ParameterHandler &prm) - { - System::parse_parameters(prm); - LinearSolver::parse_parameters(prm); - NonlinearSolver::parse_parameters(prm); - Time::parse_parameters(prm); - Discretization::parse_parameters(prm); - PreciceAdapterConfiguration::parse_parameters(prm); - } - } // namespace Parameters -} // namespace Nonlinear_Elasticity - -#endif diff --git a/source/nonlinear_elasticity/CMakeLists.txt b/source/nonlinear_elasticity/CMakeLists.txt new file mode 100644 index 0000000..3c8a2a5 --- /dev/null +++ b/source/nonlinear_elasticity/CMakeLists.txt @@ -0,0 +1,18 @@ +## +# CMake script for the dealii-adapter: +## + +# Set the name of the project and target: +SET(TARGET "nonlinear_elasticity") + +SET(_src + ${TARGET}.cc +) + +CMAKE_MINIMUM_REQUIRED(VERSION 2.8.12) + +ADD_LIBRARY(${TARGET} OBJECT ${_src}) +FIND_PACKAGE(precice REQUIRED) +TARGET_LINK_LIBRARIES(${TARGET} precice::precice) + +DEAL_II_SETUP_TARGET(${TARGET}) diff --git a/nonlinear_elasticity/include/compressible_neo_hook_material.h b/source/nonlinear_elasticity/include/compressible_neo_hook_material.h similarity index 97% rename from nonlinear_elasticity/include/compressible_neo_hook_material.h rename to source/nonlinear_elasticity/include/compressible_neo_hook_material.h index ca841aa..75ca8f3 100644 --- a/nonlinear_elasticity/include/compressible_neo_hook_material.h +++ b/source/nonlinear_elasticity/include/compressible_neo_hook_material.h @@ -1,5 +1,4 @@ -#ifndef COMPRESSIBLE_NEO_HOOK_MATERIAL_H -#define COMPRESSIBLE_NEO_HOOK_MATERIAL_H +#pragma once #include @@ -139,4 +138,3 @@ namespace Nonlinear_Elasticity } }; } // namespace Nonlinear_Elasticity -#endif // COMPRESSIBLE_NEO_HOOK_MATERIAL_H diff --git a/source/nonlinear_elasticity/include/nonlinear_elasticity.h b/source/nonlinear_elasticity/include/nonlinear_elasticity.h new file mode 100644 index 0000000..e68f240 --- /dev/null +++ b/source/nonlinear_elasticity/include/nonlinear_elasticity.h @@ -0,0 +1,338 @@ +#ifndef NONLINEAR_ELASTICITY_H +#define NONLINEAR_ELASTICITY_H + +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include + +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include + +#include +#include +#include + +#include +#include + +#include +#include + +#include "compressible_neo_hook_material.h" +#include "parameter_handling.h" +#include "postprocessor.h" + +namespace Nonlinear_Elasticity +{ + using namespace dealii; + + // PointHistory class offers a method for storing data at the quadrature + // points. Here each quadrature point holds a pointer to a material + // description + template + class PointHistory + { + public: + PointHistory() + {} + + virtual ~PointHistory() + {} + + void + setup_lqp(const Parameters::AllParameters ¶meters) + { + material.reset( + new Material_Compressible_Neo_Hook_One_Field( + parameters.mu, parameters.nu, parameters.rho)); + } + + // Strain energy + NumberType + get_Psi(const NumberType & det_F, + const SymmetricTensor<2, dim, NumberType> &b_bar) const + { + return material->get_Psi(det_F, b_bar); + } + // Kirchhoff stress + SymmetricTensor<2, dim, NumberType> + get_tau(const NumberType & det_F, + const SymmetricTensor<2, dim, NumberType> &b_bar) const + { + return material->get_tau(det_F, b_bar); + } + // Tangent + SymmetricTensor<4, dim, NumberType> + get_Jc(const NumberType & det_F, + const SymmetricTensor<2, dim, NumberType> &b_bar) const + { + return material->get_Jc(det_F, b_bar); + } + // Density + NumberType + get_rho() const + { + return material->get_rho(); + } + + private: + std::shared_ptr> + material; + }; + + + // Forward declarations for classes that will perform assembly of the + // linearized system + template + struct Assembler_Base; + template + struct Assembler; + + // The Solid class is the central class in that it represents the problem at + // hand. It follows the usual scheme in that all it really has is a + // constructor, destructor and a run() function that dispatches all the work + // to private functions of this class: + template + class Solid + { + public: + Solid(const std::string &case_path); + + virtual ~Solid(); + + void + run(); + + private: + // Generates grid and sets boundary IDs, which are needed for different BCs + void + make_grid(); + // Set up the finite element system to be solved + void + system_setup(); + + // Several functions to assemble the system and right hand side matrices + // using multithreading. Each of them comes as a wrapper function, one that + // is executed to do the work in the WorkStream model on one cell, and one + // that copies the work done on this one cell into the global object that + // represents it + void + assemble_system(const BlockVector &solution_delta, + const BlockVector &acceleration, + const BlockVector &external_stress); + + // We use a separate data structure to perform the assembly. It needs access + // to some low-level data, so we simply befriend the class instead of + // creating a complex interface to provide access as necessary. + friend struct Assembler_Base; + friend struct Assembler; + + // Apply Dirichlet boundary conditions on the displacement field + void + make_constraints(const int &it_nr); + + // Create and update the quadrature points. Here, no data needs to be copied + // into a global object, so the copy_local_to_global function is empty: + void + setup_qph(); + + // Solve for the displacement using a Newton-Raphson method. + void + solve_nonlinear_timestep(BlockVector &solution_delta); + + std::pair + solve_linear_system(BlockVector &newton_update); + + // Solution retrieval + BlockVector + get_total_solution(const BlockVector &solution_delta) const; + + // Update fnuctions for time dependent variables according to Newmarks + // scheme + void + update_acceleration(BlockVector &displacement_delta); + + void + update_velocity(BlockVector &displacement_delta); + + void + update_old_variables(); + + // Post-processing and writing data to file : + void + output_results() const; + + const Parameters::AllParameters parameters; + + double vol_reference; + double vol_current; + + Triangulation triangulation; + + CellDataStorage::cell_iterator, + PointHistory> + quadrature_point_history; + + const unsigned int degree; + const FESystem fe; + DoFHandler dof_handler_ref; + const unsigned int dofs_per_cell; + const FEValuesExtractors::Vector u_fe; + + // Description of how the block-system is arranged. There is just 1 block, + // that contains a vector DOF u. This is a legacy of the original work + // (step-44) + static const unsigned int n_blocks = 1; + static const unsigned int n_components = dim; + static const unsigned int first_u_component = 0; + + enum + { + u_dof = 0 + }; + + std::vector dofs_per_block; + + const QGauss qf_cell; + const QGauss qf_face; + const unsigned int n_q_points; + const unsigned int n_q_points_f; + // Interface ID, which is later assigned to the mesh region for coupling + // It is chosen arbotrarily + const unsigned int boundary_interface_id; + + // Newmark parameters + // Coefficients, which are needed for time dependencies + const double alpha_1 = + 1. / (parameters.beta * std::pow(parameters.delta_t, 2)); + const double alpha_2 = 1. / (parameters.beta * parameters.delta_t); + const double alpha_3 = (1 - (2 * parameters.beta)) / (2 * parameters.beta); + const double alpha_4 = + parameters.gamma / (parameters.beta * parameters.delta_t); + const double alpha_5 = 1 - (parameters.gamma / parameters.beta); + const double alpha_6 = + (1 - (parameters.gamma / (2 * parameters.beta))) * parameters.delta_t; + + // Read body force from parameter file + const Tensor<1, 3, double> body_force = parameters.body_force; + + // Clamped boundary ID to be used consistently + const unsigned int clamped_boundary_id = 1; + const unsigned int out_of_plane_clamped_mesh_id = 8; + + // ..and store the directory, in order to output the result files there + const std::string case_path; + + AffineConstraints constraints; + BlockSparsityPattern sparsity_pattern; + BlockSparseMatrix tangent_matrix; + BlockVector system_rhs; + BlockVector total_displacement; + BlockVector total_displacement_old; + BlockVector velocity; + BlockVector velocity_old; + BlockVector acceleration; + BlockVector acceleration_old; + + // Alias to collect all time dependent variables in a single vector + // This is directly passed to the Adapter routine in order to + // store these variables for implicit couplings. + std::vector *> state_variables; + + // Global vector, which keeps all contributions of the Fluid participant + // i.e. stress data for assembly. This vector is filled properly in the + // Adapter + BlockVector external_stress; + + // In order to measure some timings + mutable TimerOutput timer; + + // The main adapter objects: The time class keeps track of the current time + // and time steps. The Adapter class includes all functionalities for + // coupling via preCICE. Look at the documentation of the class for more + // information. + Adapter::Time time; + Adapter::Adapter, Parameters::AllParameters> + adapter; + + // Then define a number of variables to store norms and update norms and + // normalisation factors. + struct Errors + { + Errors() + : u(1.0) + {} + + void + reset() + { + u = 1.0; + } + void + normalise(const Errors &val) + { + if (val.u != 0.0) + u /= val.u; + } + + double u; + }; + + Errors error_residual, error_residual_0, error_residual_norm, error_update, + error_update_0, error_update_norm; + + // Methods to calculate erros + void + get_error_residual(Errors &error_residual); + + void + get_error_update(const BlockVector &newton_update, + Errors & error_update); + + // Print information to screen during simulation + static void + print_conv_header(); + + void + print_conv_footer(); + }; + +} // namespace Nonlinear_Elasticity + +#endif // NONLINEAR_ELASTICITY_H diff --git a/source/nonlinear_elasticity/include/parameter_handling.h b/source/nonlinear_elasticity/include/parameter_handling.h new file mode 100644 index 0000000..7c62131 --- /dev/null +++ b/source/nonlinear_elasticity/include/parameter_handling.h @@ -0,0 +1,145 @@ +#pragma once + +#include + +namespace Nonlinear_Elasticity +{ + using namespace dealii; + + /** + * This class declares all parameters, which can be specified in the parameter + * file. The subsection abut preCICE configurations is directly interlinked + * to the Adapter class. + */ + namespace Parameters + { + /** + * @brief System: Specifies system properties + */ + struct System + { + double nu; + double mu; + double rho; + Tensor<1, 3, double> body_force; + + static void + declare_parameters(ParameterHandler &prm); + + void + parse_parameters(ParameterHandler &prm); + }; + + + /** + * @brief LinearSolver: Specifies linear solver properties + */ + struct LinearSolver + { + std::string type_lin; + double tol_lin; + double max_iterations_lin; + + static void + declare_parameters(ParameterHandler &prm); + + void + parse_parameters(ParameterHandler &prm); + }; + + + /** + * @brief NonlinearSolver: Specifies nonlinear solver properties + */ + struct NonlinearSolver + { + unsigned int max_iterations_NR; + double tol_f; + double tol_u; + + static void + declare_parameters(ParameterHandler &prm); + + void + parse_parameters(ParameterHandler &prm); + }; + + + /** + * @brief Time: Specifies simulation time properties + */ + struct Time + { + double delta_t; + double end_time; + int output_interval; + + static void + declare_parameters(ParameterHandler &prm); + + void + parse_parameters(ParameterHandler &prm); + }; + + /** + * @brief Discretization: Specifies parameters for time integration by an + * implicit Newmark scheme and polynomial degree of the FE system + */ + struct Discretization + { + double beta; + double gamma; + unsigned int poly_degree; + + + static void + declare_parameters(ParameterHandler &prm); + + void + parse_parameters(ParameterHandler &prm); + }; + + + /** + * @brief PreciceAdapterConfiguration: Specifies preCICE related information. + * A lot of these information need to be consistent with the + * precice-config.xml file. + */ + struct PreciceAdapterConfiguration + { + std::string scenario; + std::string config_file; + std::string participant_name; + std::string mesh_name; + std::string read_data_name; + std::string write_data_name; + double flap_location; + + static void + declare_parameters(ParameterHandler &prm); + + void + parse_parameters(ParameterHandler &prm); + }; + + + struct AllParameters : public System, + public LinearSolver, + public NonlinearSolver, + public Time, + public Discretization, + public PreciceAdapterConfiguration + + { + AllParameters(const std::string &input_file); + + static void + declare_parameters(ParameterHandler &prm); + + void + parse_parameters(ParameterHandler &prm); + }; + + + } // namespace Parameters +} // namespace Nonlinear_Elasticity diff --git a/nonlinear_elasticity/include/postprocessor.h b/source/nonlinear_elasticity/include/postprocessor.h similarity index 98% rename from nonlinear_elasticity/include/postprocessor.h rename to source/nonlinear_elasticity/include/postprocessor.h index 631b20c..d4acf42 100644 --- a/nonlinear_elasticity/include/postprocessor.h +++ b/source/nonlinear_elasticity/include/postprocessor.h @@ -1,5 +1,4 @@ -#ifndef POSTPROCESSOR_H -#define POSTPROCESSOR_H +#pragma once #include @@ -122,4 +121,3 @@ namespace Nonlinear_Elasticity } } // namespace Nonlinear_Elasticity -#endif // POSTPROCESSOR_H diff --git a/nonlinear_elasticity/nonlinear_elasticity.cc b/source/nonlinear_elasticity/nonlinear_elasticity.cc similarity index 80% rename from nonlinear_elasticity/nonlinear_elasticity.cc rename to source/nonlinear_elasticity/nonlinear_elasticity.cc index 4a20e5c..bc85aa3 100644 --- a/nonlinear_elasticity/nonlinear_elasticity.cc +++ b/source/nonlinear_elasticity/nonlinear_elasticity.cc @@ -1,3 +1,5 @@ +#include "include/nonlinear_elasticity.h" + #include #include #include @@ -48,9 +50,6 @@ #include #include -#include "../adapter/adapter.h" -#include "../adapter/time.h" -#include "include/compressible_neo_hook_material.h" #include "include/parameter_handling.h" #include "include/postprocessor.h" @@ -58,278 +57,277 @@ namespace Nonlinear_Elasticity { using namespace dealii; - - // PointHistory class offers a method for storing data at the quadrature - // points. Here each quadrature point holds a pointer to a material - // description - template - class PointHistory + namespace Parameters { - public: - PointHistory() - {} - - virtual ~PointHistory() - {} - - void - setup_lqp(const Parameters::AllParameters ¶meters) + AllParameters::AllParameters(const std::string &input_file) { - material.reset( - new Material_Compressible_Neo_Hook_One_Field( - parameters.mu, parameters.nu, parameters.rho)); + ParameterHandler prm; + declare_parameters(prm); + prm.parse_input(input_file); + parse_parameters(prm); + + // Optional, if we want to print all parameters in the beginning of the + // simulation + // prm.print_parameters(std::cout,ParameterHandler::Text); } - // Strain energy - NumberType - get_Psi(const NumberType & det_F, - const SymmetricTensor<2, dim, NumberType> &b_bar) const - { - return material->get_Psi(det_F, b_bar); - } - // Kirchhoff stress - SymmetricTensor<2, dim, NumberType> - get_tau(const NumberType & det_F, - const SymmetricTensor<2, dim, NumberType> &b_bar) const - { - return material->get_tau(det_F, b_bar); - } - // Tangent - SymmetricTensor<4, dim, NumberType> - get_Jc(const NumberType & det_F, - const SymmetricTensor<2, dim, NumberType> &b_bar) const - { - return material->get_Jc(det_F, b_bar); - } - // Density - NumberType - get_rho() const + void + Discretization::declare_parameters(ParameterHandler &prm) { - return material->get_rho(); + prm.enter_subsection("Discretization"); + { + prm.declare_entry("beta", + "0.25", + Patterns::Double(0, 0.5), + "Newmark beta"); + + prm.declare_entry("gamma", + "0.5", + Patterns::Double(0, 1), + "Newmark gamma"); + + prm.declare_entry("Polynomial degree", + "3", + Patterns::Integer(0), + "Polynomial degree of the FE system"); + } + prm.leave_subsection(); } - private: - std::shared_ptr> - material; - }; - - - // Forward declarations for classes that will perform assembly of the - // linearized system - template - struct Assembler_Base; - template - struct Assembler; - - // The Solid class is the central class in that it represents the problem at - // hand. It follows the usual scheme in that all it really has is a - // constructor, destructor and a run() function that dispatches all the work - // to private functions of this class: - template - class Solid - { - public: - Solid(const std::string &case_path); - - virtual ~Solid(); - - void - run(); - - private: - // Generates grid and sets boundary IDs, which are needed for different BCs void - make_grid(); - // Set up the finite element system to be solved - void - system_setup(); - - // Several functions to assemble the system and right hand side matrices - // using multithreading. Each of them comes as a wrapper function, one that - // is executed to do the work in the WorkStream model on one cell, and one - // that copies the work done on this one cell into the global object that - // represents it - void - assemble_system(const BlockVector &solution_delta, - const BlockVector &acceleration, - const BlockVector &external_stress); - - // We use a separate data structure to perform the assembly. It needs access - // to some low-level data, so we simply befriend the class instead of - // creating a complex interface to provide access as necessary. - friend struct Assembler_Base; - friend struct Assembler; + Discretization::parse_parameters(ParameterHandler &prm) + { + prm.enter_subsection("Discretization"); + { + beta = prm.get_double("beta"); + gamma = prm.get_double("gamma"); + poly_degree = prm.get_integer("Polynomial degree"); + } + prm.leave_subsection(); + } - // Apply Dirichlet boundary conditions on the displacement field void - make_constraints(const int &it_nr); + PreciceAdapterConfiguration::declare_parameters(ParameterHandler &prm) + { + prm.enter_subsection("precice configuration"); + { + prm.declare_entry("Scenario", + "FSI3", + Patterns::Selection("FSI3|PF"), + "Cases: FSI3 or PF for perpendicular flap"); + prm.declare_entry("precice config-file", + "precice-config.xml", + Patterns::Anything(), + "Name of the precice configuration file"); + prm.declare_entry( + "Participant name", + "dealiisolver", + Patterns::Anything(), + "Name of the participant in the precice-config.xml file"); + prm.declare_entry( + "Mesh name", + "dealii-mesh", + Patterns::Anything(), + "Name of the coupling mesh in the precice-config.xml file"); + prm.declare_entry( + "Read data name", + "received-data", + Patterns::Anything(), + "Name of the read data in the precice-config.xml file"); + prm.declare_entry( + "Write data name", + "calculated-data", + Patterns::Anything(), + "Name of the write data in the precice-config.xml file"); + prm.declare_entry("Flap location", + "0.0", + Patterns::Double(-3, 3), + "PF x-location"); + } + prm.leave_subsection(); + } - // Create and update the quadrature points. Here, no data needs to be copied - // into a global object, so the copy_local_to_global function is empty: void - setup_qph(); + PreciceAdapterConfiguration::parse_parameters(ParameterHandler &prm) + { + prm.enter_subsection("precice configuration"); + { + scenario = prm.get("Scenario"); + config_file = prm.get("precice config-file"); + participant_name = prm.get("Participant name"); + mesh_name = prm.get("Mesh name"); + read_data_name = prm.get("Read data name"); + write_data_name = prm.get("Write data name"); + flap_location = prm.get_double("Flap location"); + } + prm.leave_subsection(); + } - // Solve for the displacement using a Newton-Raphson method. void - solve_nonlinear_timestep(BlockVector &solution_delta); - - std::pair - solve_linear_system(BlockVector &newton_update); + Time::declare_parameters(ParameterHandler &prm) + { + prm.enter_subsection("Time"); + { + prm.declare_entry("End time", "1", Patterns::Double(), "End time"); - // Solution retrieval - BlockVector - get_total_solution(const BlockVector &solution_delta) const; + prm.declare_entry("Time step size", + "0.1", + Patterns::Double(), + "Time step size"); - // Update fnuctions for time dependent variables according to Newmarks - // scheme - void - update_acceleration(BlockVector &displacement_delta); + prm.declare_entry("Output interval", + "1", + Patterns::Integer(), + "Output interval"); + } + prm.leave_subsection(); + } void - update_velocity(BlockVector &displacement_delta); + Time::parse_parameters(ParameterHandler &prm) + { + prm.enter_subsection("Time"); + { + end_time = prm.get_double("End time"); + delta_t = prm.get_double("Time step size"); + output_interval = prm.get_integer("Output interval"); + } + prm.leave_subsection(); + } void - update_old_variables(); + NonlinearSolver::declare_parameters(ParameterHandler &prm) + { + prm.enter_subsection("Nonlinear solver"); + { + prm.declare_entry("Max iterations Newton-Raphson", + "10", + Patterns::Integer(0), + "Number of Newton-Raphson iterations allowed"); + + prm.declare_entry("Tolerance force", + "1.0e-9", + Patterns::Double(0.0), + "Force residual tolerance"); + + prm.declare_entry("Tolerance displacement", + "1.0e-6", + Patterns::Double(0.0), + "Displacement error tolerance"); + } + prm.leave_subsection(); + } - // Post-processing and writing data to file : void - output_results() const; - - const Parameters::AllParameters parameters; - - double vol_reference; - double vol_current; - - Triangulation triangulation; - - CellDataStorage::cell_iterator, - PointHistory> - quadrature_point_history; - - const unsigned int degree; - const FESystem fe; - DoFHandler dof_handler_ref; - const unsigned int dofs_per_cell; - const FEValuesExtractors::Vector u_fe; - - // Description of how the block-system is arranged. There is just 1 block, - // that contains a vector DOF u. This is a legacy of the original work - // (step-44) - static const unsigned int n_blocks = 1; - static const unsigned int n_components = dim; - static const unsigned int first_u_component = 0; - - enum + NonlinearSolver::parse_parameters(ParameterHandler &prm) { - u_dof = 0 - }; + prm.enter_subsection("Nonlinear solver"); + { + max_iterations_NR = prm.get_integer("Max iterations Newton-Raphson"); + tol_f = prm.get_double("Tolerance force"); + tol_u = prm.get_double("Tolerance displacement"); + } + prm.leave_subsection(); + } - std::vector dofs_per_block; - - const QGauss qf_cell; - const QGauss qf_face; - const unsigned int n_q_points; - const unsigned int n_q_points_f; - // Interface ID, which is later assigned to the mesh region for coupling - // It is chosen arbotrarily - const unsigned int boundary_interface_id; - - // Newmark parameters - // Coefficients, which are needed for time dependencies - const double alpha_1 = - 1. / (parameters.beta * std::pow(parameters.delta_t, 2)); - const double alpha_2 = 1. / (parameters.beta * parameters.delta_t); - const double alpha_3 = (1 - (2 * parameters.beta)) / (2 * parameters.beta); - const double alpha_4 = - parameters.gamma / (parameters.beta * parameters.delta_t); - const double alpha_5 = 1 - (parameters.gamma / parameters.beta); - const double alpha_6 = - (1 - (parameters.gamma / (2 * parameters.beta))) * parameters.delta_t; - - // Read body force from parameter file - const Tensor<1, 3, double> body_force = parameters.body_force; - - // Clamped boundary ID to be used consistently - const unsigned int clamped_boundary_id = 1; - const unsigned int out_of_plane_clamped_mesh_id = 8; - - // ..and store the directory, in order to output the result files there - const std::string case_path; - - AffineConstraints constraints; - BlockSparsityPattern sparsity_pattern; - BlockSparseMatrix tangent_matrix; - BlockVector system_rhs; - BlockVector total_displacement; - BlockVector total_displacement_old; - BlockVector velocity; - BlockVector velocity_old; - BlockVector acceleration; - BlockVector acceleration_old; - - // Alias to collect all time dependent variables in a single vector - // This is directly passed to the Adapter routine in order to - // store these variables for implicit couplings. - std::vector *> state_variables; - - // Global vector, which keeps all contributions of the Fluid participant - // i.e. stress data for assembly. This vector is filled properly in the - // Adapter - BlockVector external_stress; - - // In order to measure some timings - mutable TimerOutput timer; - - // The main adapter objects: The time class keeps track of the current time - // and time steps. The Adapter class includes all functionalities for - // coupling via preCICE. Look at the documentation of the class for more - // information. - Adapter::Time time; - Adapter::Adapter, Parameters::AllParameters> - adapter; - - // Then define a number of variables to store norms and update norms and - // normalisation factors. - struct Errors + void + LinearSolver::declare_parameters(ParameterHandler &prm) { - Errors() - : u(1.0) - {} - - void - reset() + prm.enter_subsection("Linear solver"); { - u = 1.0; + prm.declare_entry("Solver type", + "CG", + Patterns::Selection("CG|Direct"), + "Linear solver: CG or Direct"); + + prm.declare_entry("Residual", + "1e-6", + Patterns::Double(0.0), + "Linear solver residual (scaled by residual norm)"); + + prm.declare_entry( + "Max iteration multiplier", + "1", + Patterns::Double(0.0), + "Linear solver iterations (multiples of the system matrix size)"); } - void - normalise(const Errors &val) + prm.leave_subsection(); + } + void + System::declare_parameters(ParameterHandler &prm) + { + prm.enter_subsection("System properties"); { - if (val.u != 0.0) - u /= val.u; + prm.declare_entry("Poisson's ratio", + "0.3", + Patterns::Double(-1.0, 0.5), + "Poisson's ratio"); + + prm.declare_entry("Shear modulus", + "0.4225e6", + Patterns::Double(), + "Shear modulus"); + + prm.declare_entry("rho", "1000", Patterns::Double(), "rho"); + + prm.declare_entry("body forces", + "0,0,0", + Patterns::List(Patterns::Double()), + "body forces x,y,z"); } + prm.leave_subsection(); + } - double u; - }; - - Errors error_residual, error_residual_0, error_residual_norm, error_update, - error_update_0, error_update_norm; - - // Methods to calculate erros void - get_error_residual(Errors &error_residual); + System::parse_parameters(ParameterHandler &prm) + { + prm.enter_subsection("System properties"); + { + nu = prm.get_double("Poisson's ratio"); + mu = prm.get_double("Shear modulus"); + rho = prm.get_double("rho"); + const std::vector body_forces_input = + Utilities::split_string_list(prm.get("body forces")); + for (uint d = 0; d < 3; ++d) + body_force[d] = Utilities::string_to_double(body_forces_input[d]); + } + prm.leave_subsection(); + } void - get_error_update(const BlockVector &newton_update, - Errors & error_update); - - // Print information to screen during simulation - static void - print_conv_header(); + LinearSolver::parse_parameters(ParameterHandler &prm) + { + prm.enter_subsection("Linear solver"); + { + type_lin = prm.get("Solver type"); + tol_lin = prm.get_double("Residual"); + max_iterations_lin = prm.get_double("Max iteration multiplier"); + } + prm.leave_subsection(); + } + } // namespace Parameters - void - print_conv_footer(); - }; + void + Parameters::AllParameters::declare_parameters(ParameterHandler &prm) + { + System::declare_parameters(prm); + LinearSolver::declare_parameters(prm); + NonlinearSolver::declare_parameters(prm); + Time::declare_parameters(prm); + Discretization::declare_parameters(prm); + PreciceAdapterConfiguration::declare_parameters(prm); + } + void + Parameters::AllParameters::parse_parameters(ParameterHandler &prm) + { + System::parse_parameters(prm); + LinearSolver::parse_parameters(prm); + NonlinearSolver::parse_parameters(prm); + Time::parse_parameters(prm); + Discretization::parse_parameters(prm); + PreciceAdapterConfiguration::parse_parameters(prm); + } // Constructor initializes member variables and reads the parameter file template @@ -1499,90 +1497,5 @@ namespace Nonlinear_Elasticity timer.leave_subsection("Output results"); } + template class Solid; } // namespace Nonlinear_Elasticity - -int -main(int argc, char **argv) -{ - using namespace Nonlinear_Elasticity; - using namespace dealii; - -#ifdef DEAL_II_WITH_MPI - Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv); -#endif - - try - { - deallog.depth_console(0); - static const unsigned int n_threads = MultithreadInfo::n_threads(); - - // Query adapter and deal.II info - const std::string adapter_info = - GIT_SHORTREV == std::string("") ? - "unknown" : - (GIT_SHORTREV + std::string(" on branch ") + GIT_BRANCH); - const std::string dealii_info = - DEAL_II_GIT_SHORTREV == std::string("") ? - "unknown" : - (DEAL_II_GIT_SHORTREV + std::string(" on branch ") + - DEAL_II_GIT_BRANCH); - - std::cout - << "-----------------------------------------------------------------------------" - << std::endl - << "-- . running with " << n_threads << " thread" - << (n_threads == 1 ? "" : "s") << std::endl; - - std::cout << "-- . adapter revision " << adapter_info << std::endl; - std::cout << "-- . deal.II " << DEAL_II_PACKAGE_VERSION - << " (revision " << dealii_info << ")" << std::endl; - std::cout - << "-----------------------------------------------------------------------------" - << std::endl - << std::endl; - - - std::string parameter_file; - if (argc > 1) - parameter_file = argv[1]; - else - parameter_file = "nonlinear_elasticity.prm"; - - // Extract case path for the output directory - size_t pos = parameter_file.find_last_of("/"); - std::string case_path = - std::string::npos == pos ? "" : parameter_file.substr(0, pos + 1); - - // Dimension is determinded via cmake -DDIM - Solid solid(case_path); - solid.run(); - } - catch (std::exception &exc) - { - std::cerr << std::endl - << std::endl - << "----------------------------------------------------" - << std::endl; - std::cerr << "Exception on processing: " << std::endl - << exc.what() << std::endl - << "Aborting!" << std::endl - << "----------------------------------------------------" - << std::endl; - - return 1; - } - catch (...) - { - std::cerr << std::endl - << std::endl - << "----------------------------------------------------" - << std::endl; - std::cerr << "Unknown exception!" << std::endl - << "Aborting!" << std::endl - << "----------------------------------------------------" - << std::endl; - return 1; - } - - return 0; -} diff --git a/nonlinear_elasticity/nonlinear_elasticity.prm b/source/nonlinear_elasticity/nonlinear_elasticity.prm similarity index 100% rename from nonlinear_elasticity/nonlinear_elasticity.prm rename to source/nonlinear_elasticity/nonlinear_elasticity.prm From f4a83051977ef00596aaf773cfa211666093f2c3 Mon Sep 17 00:00:00 2001 From: DavidSCN Date: Sat, 13 Feb 2021 20:45:13 +0100 Subject: [PATCH 06/18] Fix building test --- .github/workflows/building.yml | 4 ---- 1 file changed, 4 deletions(-) diff --git a/.github/workflows/building.yml b/.github/workflows/building.yml index 166e60d..9af3f96 100644 --- a/.github/workflows/building.yml +++ b/.github/workflows/building.yml @@ -22,10 +22,6 @@ jobs: cd dealii-adapter && \ git fetch origin ${{ github.ref }} && \ git checkout FETCH_HEAD && \ - cd linear_elasticity && \ - cmake . && \ - make && \ - cd ../nonlinear_elasticity && \ cmake . && \ make"; From bf653bd2d8cb77f918b244374be8d27856d118f7 Mon Sep 17 00:00:00 2001 From: DavidSCN Date: Wed, 3 Mar 2021 08:49:54 +0100 Subject: [PATCH 07/18] Move and unify parameter header file into include path --- include/adapter/parameters.h | 133 ++++++++++++++++ .../include/parameter_handling.h | 120 --------------- .../include/parameter_handling.h | 145 ------------------ 3 files changed, 133 insertions(+), 265 deletions(-) create mode 100644 include/adapter/parameters.h delete mode 100644 source/linear_elasticity/include/parameter_handling.h delete mode 100644 source/nonlinear_elasticity/include/parameter_handling.h diff --git a/include/adapter/parameters.h b/include/adapter/parameters.h new file mode 100644 index 0000000..bdad65d --- /dev/null +++ b/include/adapter/parameters.h @@ -0,0 +1,133 @@ +#ifndef PARAMETERS_H +#define PARAMETERS_H + +#include + +/** + * This class declares all parameters, which can be specified in the parameter + * file. The subsection abut preCICE configurations is directly interlinked + * to the Adapter class. + */ +namespace Parameters +{ + using namespace dealii; + /** + * @brief Time: Specifies simulation time properties + */ + struct Time + { + double end_time; + double delta_t; + int output_interval; + + static void + declare_parameters(ParameterHandler &prm); + + void + parse_parameters(ParameterHandler &prm); + }; + + /** + * @brief The System struct keeps material properties and body force contributions + */ + struct System + { + double nu; + double mu; + double lambda; + double rho; + Tensor<1, 3, double> body_force; + + + static void + declare_parameters(ParameterHandler &prm); + + void + parse_parameters(ParameterHandler &prm); + }; + + + /** + * @brief LinearSolver: Specifies linear solver properties + */ + struct Solver + { + std::string type_lin; + double tol_lin; + double max_iterations_lin; + unsigned int max_iterations_NR; + double tol_f; + double tol_u; + + static void + declare_parameters(ParameterHandler &prm); + + void + parse_parameters(ParameterHandler &prm); + }; + + + + /** + * @brief Discretization: Specifies parameters for time integration by a + * theta scheme and polynomial degree of the FE system + */ + struct Discretization + { + unsigned int poly_degree; + // For the linear elastic model (theta-scheme) + double theta; + // For the nonlinear elastic model (Newmark) + double beta; + double gamma ; + + static void + declare_parameters(ParameterHandler &prm); + + void + parse_parameters(ParameterHandler &prm); + }; + + + /** + * @brief PreciceAdapterConfiguration: Specifies preCICE related information. + * A lot of these information need to be consistent with the + * precice-config.xml file. + */ + struct PreciceAdapterConfiguration + { + std::string scenario; + std::string config_file; + std::string participant_name; + std::string mesh_name; + std::string read_data_name; + std::string write_data_name; + double flap_location; + bool data_consistent; + + static void + declare_parameters(ParameterHandler &prm); + + void + parse_parameters(ParameterHandler &prm); + }; + + + struct AllParameters : public Solver, + public Discretization, + public System, + public Time, + public PreciceAdapterConfiguration + + { + AllParameters(const std::string &input_file); + + static void + declare_parameters(ParameterHandler &prm); + + void + parse_parameters(ParameterHandler &prm); + }; +} // namespace Parameters + +#endif // PARAMETERS_H diff --git a/source/linear_elasticity/include/parameter_handling.h b/source/linear_elasticity/include/parameter_handling.h deleted file mode 100644 index 3a22274..0000000 --- a/source/linear_elasticity/include/parameter_handling.h +++ /dev/null @@ -1,120 +0,0 @@ -#pragma once - -#include - -namespace Linear_Elasticity -{ - using namespace dealii; - - /** - * This class declares all parameters, which can be specified in the parameter - * file. The subsection abut preCICE configurations is directly interlinked - * to the Adapter class. - */ - namespace Parameters - { - // TODO: Add more parameters - struct Time - { - double delta_t; - double end_time; - int output_interval; - - static void - declare_parameters(ParameterHandler &prm); - - void - parse_parameters(ParameterHandler &prm); - }; - - /** - * @brief Discretization: Specifies parameters for time integration by a - * theta scheme and polynomial degree of the FE system - */ - struct Discretization - { - double theta; - unsigned int poly_degree; - - static void - declare_parameters(ParameterHandler &prm); - - void - parse_parameters(ParameterHandler &prm); - }; - - /** - * @brief The System struct keeps material properties and body force contributions - */ - struct System - { - double mu; - double lambda; - double rho; - Tensor<1, 3, double> body_force; - - - static void - declare_parameters(ParameterHandler &prm); - - void - parse_parameters(ParameterHandler &prm); - }; - - /** - * @brief LinearSolver: Specifies linear solver properties - */ - struct LinearSolver - { - std::string type_lin; - double tol_lin; - double max_iterations_lin; - - static void - declare_parameters(ParameterHandler &prm); - - void - parse_parameters(ParameterHandler &prm); - }; - - /** - * @brief PreciceAdapterConfiguration: Specifies preCICE related information. - * A lot of these information need to be consistent with the - * precice-config.xml file. - */ - struct PreciceAdapterConfiguration - { - std::string scenario; - std::string config_file; - std::string participant_name; - std::string mesh_name; - std::string read_data_name; - std::string write_data_name; - double flap_location; - bool data_consistent; - - static void - declare_parameters(ParameterHandler &prm); - - void - parse_parameters(ParameterHandler &prm); - }; - - - struct AllParameters : public LinearSolver, - public Discretization, - public System, - public Time, - public PreciceAdapterConfiguration - - { - AllParameters(const std::string &input_file); - - static void - declare_parameters(ParameterHandler &prm); - - void - parse_parameters(ParameterHandler &prm); - }; - } // namespace Parameters -} // namespace Linear_Elasticity diff --git a/source/nonlinear_elasticity/include/parameter_handling.h b/source/nonlinear_elasticity/include/parameter_handling.h deleted file mode 100644 index 7c62131..0000000 --- a/source/nonlinear_elasticity/include/parameter_handling.h +++ /dev/null @@ -1,145 +0,0 @@ -#pragma once - -#include - -namespace Nonlinear_Elasticity -{ - using namespace dealii; - - /** - * This class declares all parameters, which can be specified in the parameter - * file. The subsection abut preCICE configurations is directly interlinked - * to the Adapter class. - */ - namespace Parameters - { - /** - * @brief System: Specifies system properties - */ - struct System - { - double nu; - double mu; - double rho; - Tensor<1, 3, double> body_force; - - static void - declare_parameters(ParameterHandler &prm); - - void - parse_parameters(ParameterHandler &prm); - }; - - - /** - * @brief LinearSolver: Specifies linear solver properties - */ - struct LinearSolver - { - std::string type_lin; - double tol_lin; - double max_iterations_lin; - - static void - declare_parameters(ParameterHandler &prm); - - void - parse_parameters(ParameterHandler &prm); - }; - - - /** - * @brief NonlinearSolver: Specifies nonlinear solver properties - */ - struct NonlinearSolver - { - unsigned int max_iterations_NR; - double tol_f; - double tol_u; - - static void - declare_parameters(ParameterHandler &prm); - - void - parse_parameters(ParameterHandler &prm); - }; - - - /** - * @brief Time: Specifies simulation time properties - */ - struct Time - { - double delta_t; - double end_time; - int output_interval; - - static void - declare_parameters(ParameterHandler &prm); - - void - parse_parameters(ParameterHandler &prm); - }; - - /** - * @brief Discretization: Specifies parameters for time integration by an - * implicit Newmark scheme and polynomial degree of the FE system - */ - struct Discretization - { - double beta; - double gamma; - unsigned int poly_degree; - - - static void - declare_parameters(ParameterHandler &prm); - - void - parse_parameters(ParameterHandler &prm); - }; - - - /** - * @brief PreciceAdapterConfiguration: Specifies preCICE related information. - * A lot of these information need to be consistent with the - * precice-config.xml file. - */ - struct PreciceAdapterConfiguration - { - std::string scenario; - std::string config_file; - std::string participant_name; - std::string mesh_name; - std::string read_data_name; - std::string write_data_name; - double flap_location; - - static void - declare_parameters(ParameterHandler &prm); - - void - parse_parameters(ParameterHandler &prm); - }; - - - struct AllParameters : public System, - public LinearSolver, - public NonlinearSolver, - public Time, - public Discretization, - public PreciceAdapterConfiguration - - { - AllParameters(const std::string &input_file); - - static void - declare_parameters(ParameterHandler &prm); - - void - parse_parameters(ParameterHandler &prm); - }; - - - } // namespace Parameters -} // namespace Nonlinear_Elasticity From c3fe0064ddfbfaa4ee6cbd57ca409892044f7414 Mon Sep 17 00:00:00 2001 From: DavidSCN Date: Wed, 3 Mar 2021 08:51:56 +0100 Subject: [PATCH 08/18] Move and unify parameter sources into include path --- include/adapter/parameters.cc | 277 ++++++++++++++++++ source/linear_elasticity/linear_elasticity.cc | 237 --------------- .../nonlinear_elasticity.cc | 272 ----------------- 3 files changed, 277 insertions(+), 509 deletions(-) create mode 100644 include/adapter/parameters.cc diff --git a/include/adapter/parameters.cc b/include/adapter/parameters.cc new file mode 100644 index 0000000..7a6a4aa --- /dev/null +++ b/include/adapter/parameters.cc @@ -0,0 +1,277 @@ +#include + +namespace Parameters +{ + void + Time::declare_parameters(ParameterHandler &prm) + { + prm.enter_subsection("Time"); + { + prm.declare_entry("End time", "1", Patterns::Double(), "End time"); + + prm.declare_entry("Time step size", + "0.1", + Patterns::Double(), + "Time step size"); + + prm.declare_entry("Output interval", + "1", + Patterns::Integer(0), + "Write results every x timesteps"); + } + prm.leave_subsection(); + } + + void + Time::parse_parameters(ParameterHandler &prm) + { + prm.enter_subsection("Time"); + { + end_time = prm.get_double("End time"); + delta_t = prm.get_double("Time step size"); + output_interval = prm.get_integer("Output interval"); + } + prm.leave_subsection(); + } + + void + System::declare_parameters(ParameterHandler &prm) + { + prm.enter_subsection("System properties"); + { + prm.declare_entry("Shear modulus", + "0.4225e6", + Patterns::Double(), + "Shear modulus"); + + prm.declare_entry("Poisson's ratio", + "0.3", + Patterns::Double(-1.0, 0.5), + "Poisson's ratio"); + + prm.declare_entry("rho", "1000", Patterns::Double(0.0), "density"); + + prm.declare_entry("body forces", + "0,0,0", + Patterns::List(Patterns::Double()), + "body forces x,y,z"); + } + prm.leave_subsection(); + } + + void + System::parse_parameters(ParameterHandler &prm) + { + prm.enter_subsection("System properties"); + { + mu = prm.get_double("Shear modulus"); + nu = prm.get_double("Poisson's ratio"); + rho = prm.get_double("rho"); + const std::vector body_forces_input = + Utilities::split_string_list(prm.get("body forces")); + for (uint d = 0; d < 3; ++d) + body_force[d] = Utilities::string_to_double(body_forces_input[d]); + } + prm.leave_subsection(); + + lambda = 2 * mu * nu / (1 - 2 * nu); + } + + + void + Solver::declare_parameters(ParameterHandler &prm) + { + prm.enter_subsection("Linear solver"); + { + prm.declare_entry("Solver type", + "Direct", + Patterns::Selection("CG|Direct"), + "Linear solver: CG or Direct"); + + prm.declare_entry("Residual", + "1e-6", + Patterns::Double(0.0), + "Linear solver residual (scaled by residual norm)"); + + prm.declare_entry( + "Max iteration multiplier", + "1", + Patterns::Double(0.0), + "Linear solver iterations (multiples of the system matrix size)"); + + prm.declare_entry("Max iterations Newton-Raphson", + "10", + Patterns::Integer(0), + "Number of Newton-Raphson iterations allowed"); + + prm.declare_entry("Tolerance force", + "1.0e-9", + Patterns::Double(0.0), + "Force residual tolerance"); + + prm.declare_entry("Tolerance displacement", + "1.0e-6", + Patterns::Double(0.0), + "Displacement error tolerance"); + } + prm.leave_subsection(); + } + + + void + Solver::parse_parameters(ParameterHandler &prm) + { + prm.enter_subsection("Solver"); + { + type_lin = prm.get("Solver type"); + tol_lin = prm.get_double("Residual"); + max_iterations_lin = prm.get_double("Max iteration multiplier"); + + max_iterations_NR = prm.get_integer("Max iterations Newton-Raphson"); + tol_f = prm.get_double("Tolerance force"); + tol_u = prm.get_double("Tolerance displacement"); + } + prm.leave_subsection(); + } + + + void + Discretization::declare_parameters(ParameterHandler &prm) + { + prm.enter_subsection("Discretization"); + { + prm.declare_entry("Polynomial degree", + "3", + Patterns::Integer(0), + "Polynomial degree of the FE system"); + + prm.declare_entry("theta", + "0.5", + Patterns::Double(0, 1), + "Time integration scheme"); + + prm.declare_entry("beta", + "0.25", + Patterns::Double(0, 0.5), + "Newmark beta"); + + prm.declare_entry("gamma", + "0.5", + Patterns::Double(0, 1), + "Newmark gamma"); + } + prm.leave_subsection(); + } + + void + Discretization::parse_parameters(ParameterHandler &prm) + { + prm.enter_subsection("Discretization"); + { + poly_degree = prm.get_integer("Polynomial degree"); + theta = prm.get_double("theta"); + beta = prm.get_double("beta"); + gamma = prm.get_double("gamma"); + } + prm.leave_subsection(); + } + + + void + PreciceAdapterConfiguration::declare_parameters(ParameterHandler &prm) + { + prm.enter_subsection("precice configuration"); + { + prm.declare_entry("Scenario", + "FSI3", + Patterns::Selection("FSI3|PF"), + "Cases: FSI3 or PF for perpendicular flap"); + prm.declare_entry("precice config-file", + "precice-config.xml", + Patterns::Anything(), + "Name of the precice configuration file"); + prm.declare_entry( + "Participant name", + "dealiisolver", + Patterns::Anything(), + "Name of the participant in the precice-config.xml file"); + prm.declare_entry( + "Mesh name", + "dealii-mesh", + Patterns::Anything(), + "Name of the coupling mesh in the precice-config.xml file"); + prm.declare_entry("Read data name", + "received-data", + Patterns::Anything(), + "Name of the read data in the precice-config.xml file"); + prm.declare_entry( + "Write data name", + "calculated-data", + Patterns::Anything(), + "Name of the write data in the precice-config.xml file"); + prm.declare_entry("Flap location", + "0.0", + Patterns::Double(-3, 3), + "PF x-location"); + } + prm.leave_subsection(); + } + + void + PreciceAdapterConfiguration::parse_parameters(ParameterHandler &prm) + { + prm.enter_subsection("precice configuration"); + { + scenario = prm.get("Scenario"); + config_file = prm.get("precice config-file"); + participant_name = prm.get("Participant name"); + mesh_name = prm.get("Mesh name"); + read_data_name = prm.get("Read data name"); + write_data_name = prm.get("Write data name"); + flap_location = prm.get_double("Flap location"); + } + prm.leave_subsection(); + // Look at the specific type of read data + if ((read_data_name.find("Stress") == 0)) + data_consistent = true; + else if ((read_data_name.find("Force") == 0)) + data_consistent = false; + else + AssertThrow( + false, + ExcMessage( + "Unknown read data type. Please use 'Force' or 'Stress' in the read data naming.")); + } + + AllParameters::AllParameters(const std::string &input_file) + { + ParameterHandler prm; + declare_parameters(prm); + prm.parse_input(input_file); + parse_parameters(prm); + + // Optional, if we want to print all parameters in the beginning of the + // simulation + // prm.print_parameters(std::cout,ParameterHandler::Text); + } + + void + AllParameters::declare_parameters(ParameterHandler &prm) + { + Solver::declare_parameters(prm); + Discretization::declare_parameters(prm); + System::declare_parameters(prm); + Time::declare_parameters(prm); + PreciceAdapterConfiguration::declare_parameters(prm); + } + + void + AllParameters::parse_parameters(ParameterHandler &prm) + { + Solver::parse_parameters(prm); + Discretization::parse_parameters(prm); + System::parse_parameters(prm); + Time::parse_parameters(prm); + PreciceAdapterConfiguration::parse_parameters(prm); + } +} // namespace Parameters diff --git a/source/linear_elasticity/linear_elasticity.cc b/source/linear_elasticity/linear_elasticity.cc index 6dc0ef3..a69584d 100644 --- a/source/linear_elasticity/linear_elasticity.cc +++ b/source/linear_elasticity/linear_elasticity.cc @@ -50,243 +50,6 @@ namespace Linear_Elasticity { using namespace dealii; - namespace Parameters - { - void - Time::declare_parameters(ParameterHandler &prm) - { - prm.enter_subsection("Time"); - { - prm.declare_entry("End time", "1", Patterns::Double(), "End time"); - - prm.declare_entry("Time step size", - "0.1", - Patterns::Double(), - "Time step size"); - - prm.declare_entry("Output interval", - "1", - Patterns::Integer(0), - "Write results every x timesteps"); - } - prm.leave_subsection(); - } - - void - Time::parse_parameters(ParameterHandler &prm) - { - prm.enter_subsection("Time"); - { - end_time = prm.get_double("End time"); - delta_t = prm.get_double("Time step size"); - output_interval = prm.get_integer("Output interval"); - } - prm.leave_subsection(); - } - - void - Discretization::declare_parameters(ParameterHandler &prm) - { - prm.enter_subsection("Discretization"); - { - prm.declare_entry("theta", - "0.5", - Patterns::Double(0, 1), - "Time integration scheme"); - - prm.declare_entry("Polynomial degree", - "3", - Patterns::Integer(0), - "Polynomial degree of the FE system"); - } - prm.leave_subsection(); - } - - void - Discretization::parse_parameters(ParameterHandler &prm) - { - prm.enter_subsection("Discretization"); - { - theta = prm.get_double("theta"); - poly_degree = prm.get_integer("Polynomial degree"); - } - prm.leave_subsection(); - } - - void - System::declare_parameters(ParameterHandler &prm) - { - prm.enter_subsection("System properties"); - { - prm.declare_entry("mu", "0.5e6", Patterns::Double(), "mu"); - - prm.declare_entry("lambda", "2e6", Patterns::Double(), "lambda"); - - prm.declare_entry("rho", "1000", Patterns::Double(0.0), "density"); - - prm.declare_entry("body forces", - "0,0,0", - Patterns::List(Patterns::Double()), - "body forces x,y,z"); - } - prm.leave_subsection(); - } - - void - System::parse_parameters(ParameterHandler &prm) - { - prm.enter_subsection("System properties"); - { - mu = prm.get_double("mu"); - lambda = prm.get_double("lambda"); - rho = prm.get_double("rho"); - const std::vector body_forces_input = - Utilities::split_string_list(prm.get("body forces")); - for (uint d = 0; d < 3; ++d) - body_force[d] = Utilities::string_to_double(body_forces_input[d]); - } - prm.leave_subsection(); - } - - - void - LinearSolver::declare_parameters(ParameterHandler &prm) - { - prm.enter_subsection("Linear solver"); - { - prm.declare_entry("Solver type", - "Direct", - Patterns::Selection("CG|Direct"), - "Linear solver: CG or Direct"); - - prm.declare_entry("Residual", - "1e-6", - Patterns::Double(0.0), - "Linear solver residual (scaled by residual norm)"); - - prm.declare_entry( - "Max iteration multiplier", - "1", - Patterns::Double(0.0), - "Linear solver iterations (multiples of the system matrix size)"); - } - prm.leave_subsection(); - } - - void - LinearSolver::parse_parameters(ParameterHandler &prm) - { - prm.enter_subsection("Linear solver"); - { - type_lin = prm.get("Solver type"); - tol_lin = prm.get_double("Residual"); - max_iterations_lin = prm.get_double("Max iteration multiplier"); - } - prm.leave_subsection(); - } - - void - PreciceAdapterConfiguration::declare_parameters(ParameterHandler &prm) - { - prm.enter_subsection("precice configuration"); - { - prm.declare_entry("Scenario", - "FSI3", - Patterns::Selection("FSI3|PF"), - "Cases: FSI3 or PF for perpendicular flap"); - prm.declare_entry("precice config-file", - "precice-config.xml", - Patterns::Anything(), - "Name of the precice configuration file"); - prm.declare_entry( - "Participant name", - "dealiisolver", - Patterns::Anything(), - "Name of the participant in the precice-config.xml file"); - prm.declare_entry( - "Mesh name", - "dealii-mesh", - Patterns::Anything(), - "Name of the coupling mesh in the precice-config.xml file"); - prm.declare_entry( - "Read data name", - "received-data", - Patterns::Anything(), - "Name of the read data in the precice-config.xml file"); - prm.declare_entry( - "Write data name", - "calculated-data", - Patterns::Anything(), - "Name of the write data in the precice-config.xml file"); - prm.declare_entry("Flap location", - "0.0", - Patterns::Double(-3, 3), - "PF x-location"); - } - prm.leave_subsection(); - } - - void - PreciceAdapterConfiguration::parse_parameters(ParameterHandler &prm) - { - prm.enter_subsection("precice configuration"); - { - scenario = prm.get("Scenario"); - config_file = prm.get("precice config-file"); - participant_name = prm.get("Participant name"); - mesh_name = prm.get("Mesh name"); - read_data_name = prm.get("Read data name"); - write_data_name = prm.get("Write data name"); - flap_location = prm.get_double("Flap location"); - } - prm.leave_subsection(); - // Look at the specific type of read data - if ((read_data_name.find("Stress") == 0)) - data_consistent = true; - else if ((read_data_name.find("Force") == 0)) - data_consistent = false; - else - AssertThrow( - false, - ExcMessage( - "Unknown read data type. Please use 'Force' or 'Stress' in the read data naming.")); - } - - AllParameters::AllParameters(const std::string &input_file) - { - ParameterHandler prm; - declare_parameters(prm); - prm.parse_input(input_file); - parse_parameters(prm); - - // Optional, if we want to print all parameters in the beginning of the - // simulation - // prm.print_parameters(std::cout,ParameterHandler::Text); - } - - void - AllParameters::declare_parameters(ParameterHandler &prm) - { - LinearSolver::declare_parameters(prm); - Discretization::declare_parameters(prm); - System::declare_parameters(prm); - Time::declare_parameters(prm); - PreciceAdapterConfiguration::declare_parameters(prm); - } - - void - AllParameters::parse_parameters(ParameterHandler &prm) - { - LinearSolver::parse_parameters(prm); - Discretization::parse_parameters(prm); - System::parse_parameters(prm); - Time::parse_parameters(prm); - PreciceAdapterConfiguration::parse_parameters(prm); - } - - } // namespace Parameters - - // Constructor template ElastoDynamics::ElastoDynamics(const std::string &case_path) diff --git a/source/nonlinear_elasticity/nonlinear_elasticity.cc b/source/nonlinear_elasticity/nonlinear_elasticity.cc index bc85aa3..2425d3b 100644 --- a/source/nonlinear_elasticity/nonlinear_elasticity.cc +++ b/source/nonlinear_elasticity/nonlinear_elasticity.cc @@ -57,278 +57,6 @@ namespace Nonlinear_Elasticity { using namespace dealii; - namespace Parameters - { - AllParameters::AllParameters(const std::string &input_file) - { - ParameterHandler prm; - declare_parameters(prm); - prm.parse_input(input_file); - parse_parameters(prm); - - // Optional, if we want to print all parameters in the beginning of the - // simulation - // prm.print_parameters(std::cout,ParameterHandler::Text); - } - - void - Discretization::declare_parameters(ParameterHandler &prm) - { - prm.enter_subsection("Discretization"); - { - prm.declare_entry("beta", - "0.25", - Patterns::Double(0, 0.5), - "Newmark beta"); - - prm.declare_entry("gamma", - "0.5", - Patterns::Double(0, 1), - "Newmark gamma"); - - prm.declare_entry("Polynomial degree", - "3", - Patterns::Integer(0), - "Polynomial degree of the FE system"); - } - prm.leave_subsection(); - } - - void - Discretization::parse_parameters(ParameterHandler &prm) - { - prm.enter_subsection("Discretization"); - { - beta = prm.get_double("beta"); - gamma = prm.get_double("gamma"); - poly_degree = prm.get_integer("Polynomial degree"); - } - prm.leave_subsection(); - } - - void - PreciceAdapterConfiguration::declare_parameters(ParameterHandler &prm) - { - prm.enter_subsection("precice configuration"); - { - prm.declare_entry("Scenario", - "FSI3", - Patterns::Selection("FSI3|PF"), - "Cases: FSI3 or PF for perpendicular flap"); - prm.declare_entry("precice config-file", - "precice-config.xml", - Patterns::Anything(), - "Name of the precice configuration file"); - prm.declare_entry( - "Participant name", - "dealiisolver", - Patterns::Anything(), - "Name of the participant in the precice-config.xml file"); - prm.declare_entry( - "Mesh name", - "dealii-mesh", - Patterns::Anything(), - "Name of the coupling mesh in the precice-config.xml file"); - prm.declare_entry( - "Read data name", - "received-data", - Patterns::Anything(), - "Name of the read data in the precice-config.xml file"); - prm.declare_entry( - "Write data name", - "calculated-data", - Patterns::Anything(), - "Name of the write data in the precice-config.xml file"); - prm.declare_entry("Flap location", - "0.0", - Patterns::Double(-3, 3), - "PF x-location"); - } - prm.leave_subsection(); - } - - void - PreciceAdapterConfiguration::parse_parameters(ParameterHandler &prm) - { - prm.enter_subsection("precice configuration"); - { - scenario = prm.get("Scenario"); - config_file = prm.get("precice config-file"); - participant_name = prm.get("Participant name"); - mesh_name = prm.get("Mesh name"); - read_data_name = prm.get("Read data name"); - write_data_name = prm.get("Write data name"); - flap_location = prm.get_double("Flap location"); - } - prm.leave_subsection(); - } - - void - Time::declare_parameters(ParameterHandler &prm) - { - prm.enter_subsection("Time"); - { - prm.declare_entry("End time", "1", Patterns::Double(), "End time"); - - prm.declare_entry("Time step size", - "0.1", - Patterns::Double(), - "Time step size"); - - prm.declare_entry("Output interval", - "1", - Patterns::Integer(), - "Output interval"); - } - prm.leave_subsection(); - } - - void - Time::parse_parameters(ParameterHandler &prm) - { - prm.enter_subsection("Time"); - { - end_time = prm.get_double("End time"); - delta_t = prm.get_double("Time step size"); - output_interval = prm.get_integer("Output interval"); - } - prm.leave_subsection(); - } - - void - NonlinearSolver::declare_parameters(ParameterHandler &prm) - { - prm.enter_subsection("Nonlinear solver"); - { - prm.declare_entry("Max iterations Newton-Raphson", - "10", - Patterns::Integer(0), - "Number of Newton-Raphson iterations allowed"); - - prm.declare_entry("Tolerance force", - "1.0e-9", - Patterns::Double(0.0), - "Force residual tolerance"); - - prm.declare_entry("Tolerance displacement", - "1.0e-6", - Patterns::Double(0.0), - "Displacement error tolerance"); - } - prm.leave_subsection(); - } - - void - NonlinearSolver::parse_parameters(ParameterHandler &prm) - { - prm.enter_subsection("Nonlinear solver"); - { - max_iterations_NR = prm.get_integer("Max iterations Newton-Raphson"); - tol_f = prm.get_double("Tolerance force"); - tol_u = prm.get_double("Tolerance displacement"); - } - prm.leave_subsection(); - } - - void - LinearSolver::declare_parameters(ParameterHandler &prm) - { - prm.enter_subsection("Linear solver"); - { - prm.declare_entry("Solver type", - "CG", - Patterns::Selection("CG|Direct"), - "Linear solver: CG or Direct"); - - prm.declare_entry("Residual", - "1e-6", - Patterns::Double(0.0), - "Linear solver residual (scaled by residual norm)"); - - prm.declare_entry( - "Max iteration multiplier", - "1", - Patterns::Double(0.0), - "Linear solver iterations (multiples of the system matrix size)"); - } - prm.leave_subsection(); - } - void - System::declare_parameters(ParameterHandler &prm) - { - prm.enter_subsection("System properties"); - { - prm.declare_entry("Poisson's ratio", - "0.3", - Patterns::Double(-1.0, 0.5), - "Poisson's ratio"); - - prm.declare_entry("Shear modulus", - "0.4225e6", - Patterns::Double(), - "Shear modulus"); - - prm.declare_entry("rho", "1000", Patterns::Double(), "rho"); - - prm.declare_entry("body forces", - "0,0,0", - Patterns::List(Patterns::Double()), - "body forces x,y,z"); - } - prm.leave_subsection(); - } - - void - System::parse_parameters(ParameterHandler &prm) - { - prm.enter_subsection("System properties"); - { - nu = prm.get_double("Poisson's ratio"); - mu = prm.get_double("Shear modulus"); - rho = prm.get_double("rho"); - const std::vector body_forces_input = - Utilities::split_string_list(prm.get("body forces")); - for (uint d = 0; d < 3; ++d) - body_force[d] = Utilities::string_to_double(body_forces_input[d]); - } - prm.leave_subsection(); - } - - void - LinearSolver::parse_parameters(ParameterHandler &prm) - { - prm.enter_subsection("Linear solver"); - { - type_lin = prm.get("Solver type"); - tol_lin = prm.get_double("Residual"); - max_iterations_lin = prm.get_double("Max iteration multiplier"); - } - prm.leave_subsection(); - } - } // namespace Parameters - - void - Parameters::AllParameters::declare_parameters(ParameterHandler &prm) - { - System::declare_parameters(prm); - LinearSolver::declare_parameters(prm); - NonlinearSolver::declare_parameters(prm); - Time::declare_parameters(prm); - Discretization::declare_parameters(prm); - PreciceAdapterConfiguration::declare_parameters(prm); - } - - void - Parameters::AllParameters::parse_parameters(ParameterHandler &prm) - { - System::parse_parameters(prm); - LinearSolver::parse_parameters(prm); - NonlinearSolver::parse_parameters(prm); - Time::parse_parameters(prm); - Discretization::parse_parameters(prm); - PreciceAdapterConfiguration::parse_parameters(prm); - } - // Constructor initializes member variables and reads the parameter file template Solid::Solid(const std::string &case_path) From 4057bbe34218b69caeef9302383f5f67fb711be1 Mon Sep 17 00:00:00 2001 From: DavidSCN Date: Wed, 3 Mar 2021 08:52:40 +0100 Subject: [PATCH 09/18] Configure CMake for new source file and adjust includes --- CMakeLists.txt | 2 +- source/linear_elasticity/include/linear_elasticity.h | 2 +- source/linear_elasticity/linear_elasticity.cc | 2 +- source/nonlinear_elasticity/include/nonlinear_elasticity.h | 2 +- source/nonlinear_elasticity/nonlinear_elasticity.cc | 3 ++- 5 files changed, 6 insertions(+), 5 deletions(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index 97e278a..fd8f940 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -19,7 +19,7 @@ ADD_DEFINITIONS(-DDIM=${DIM}) # Set the target and the target source SET( TARGET "dealii-precice" ) -SET( TARGET_SRC ${TARGET}.cc ) +SET( TARGET_SRC ${TARGET}.cc include/adapter/parameters.cc) DEAL_II_INITIALIZE_CACHED_VARIABLES() diff --git a/source/linear_elasticity/include/linear_elasticity.h b/source/linear_elasticity/include/linear_elasticity.h index 18c2c2e..78ac2b0 100644 --- a/source/linear_elasticity/include/linear_elasticity.h +++ b/source/linear_elasticity/include/linear_elasticity.h @@ -37,12 +37,12 @@ #include #include +#include #include #include #include -#include "parameter_handling.h" #include "postprocessor.h" // The Linear_Elasticity case includes a linear elastic material with a one-step diff --git a/source/linear_elasticity/linear_elasticity.cc b/source/linear_elasticity/linear_elasticity.cc index a69584d..0120e1d 100644 --- a/source/linear_elasticity/linear_elasticity.cc +++ b/source/linear_elasticity/linear_elasticity.cc @@ -36,12 +36,12 @@ #include #include +#include #include #include #include -#include "include/parameter_handling.h" #include "include/postprocessor.h" // The Linear_Elasticity case includes a linear elastic material with a one-step diff --git a/source/nonlinear_elasticity/include/nonlinear_elasticity.h b/source/nonlinear_elasticity/include/nonlinear_elasticity.h index e68f240..a4a1270 100644 --- a/source/nonlinear_elasticity/include/nonlinear_elasticity.h +++ b/source/nonlinear_elasticity/include/nonlinear_elasticity.h @@ -49,13 +49,13 @@ #include #include +#include #include #include #include #include "compressible_neo_hook_material.h" -#include "parameter_handling.h" #include "postprocessor.h" namespace Nonlinear_Elasticity diff --git a/source/nonlinear_elasticity/nonlinear_elasticity.cc b/source/nonlinear_elasticity/nonlinear_elasticity.cc index 2425d3b..19ea161 100644 --- a/source/nonlinear_elasticity/nonlinear_elasticity.cc +++ b/source/nonlinear_elasticity/nonlinear_elasticity.cc @@ -47,10 +47,11 @@ #include #include +#include + #include #include -#include "include/parameter_handling.h" #include "include/postprocessor.h" namespace Nonlinear_Elasticity From 8f8369c7caf3a5dfe28ac56a4290e9180e392a13 Mon Sep 17 00:00:00 2001 From: DavidSCN Date: Wed, 3 Mar 2021 09:57:43 +0100 Subject: [PATCH 10/18] Simplify parameter handling using add_parameter --- include/adapter/parameters.cc | 316 +++++++++++++--------------------- include/adapter/parameters.h | 82 ++++----- 2 files changed, 146 insertions(+), 252 deletions(-) diff --git a/include/adapter/parameters.cc b/include/adapter/parameters.cc index 7a6a4aa..97b24fa 100644 --- a/include/adapter/parameters.cc +++ b/include/adapter/parameters.cc @@ -3,74 +3,47 @@ namespace Parameters { void - Time::declare_parameters(ParameterHandler &prm) + Time::add_output_parameters(ParameterHandler &prm) { prm.enter_subsection("Time"); { - prm.declare_entry("End time", "1", Patterns::Double(), "End time"); + prm.add_parameter("End time", end_time, "End time", Patterns::Double()); - prm.declare_entry("Time step size", - "0.1", - Patterns::Double(), - "Time step size"); + prm.add_parameter("Time step size", + delta_t, + "Time step size", + Patterns::Double()); - prm.declare_entry("Output interval", - "1", - Patterns::Integer(0), - "Write results every x timesteps"); + prm.add_parameter("Output interval", + output_interval, + "Write results every x timesteps", + Patterns::Integer(0)); } prm.leave_subsection(); } - void - Time::parse_parameters(ParameterHandler &prm) - { - prm.enter_subsection("Time"); - { - end_time = prm.get_double("End time"); - delta_t = prm.get_double("Time step size"); - output_interval = prm.get_integer("Output interval"); - } - prm.leave_subsection(); - } - - void - System::declare_parameters(ParameterHandler &prm) - { - prm.enter_subsection("System properties"); - { - prm.declare_entry("Shear modulus", - "0.4225e6", - Patterns::Double(), - "Shear modulus"); - - prm.declare_entry("Poisson's ratio", - "0.3", - Patterns::Double(-1.0, 0.5), - "Poisson's ratio"); - - prm.declare_entry("rho", "1000", Patterns::Double(0.0), "density"); - - prm.declare_entry("body forces", - "0,0,0", - Patterns::List(Patterns::Double()), - "body forces x,y,z"); - } - prm.leave_subsection(); - } void - System::parse_parameters(ParameterHandler &prm) + System::add_output_parameters(ParameterHandler &prm) { prm.enter_subsection("System properties"); { - mu = prm.get_double("Shear modulus"); - nu = prm.get_double("Poisson's ratio"); - rho = prm.get_double("rho"); - const std::vector body_forces_input = - Utilities::split_string_list(prm.get("body forces")); - for (uint d = 0; d < 3; ++d) - body_force[d] = Utilities::string_to_double(body_forces_input[d]); + prm.add_parameter("Shear modulus", + mu, + "Shear modulus", + Patterns::Double()); + + prm.add_parameter("Poisson's ratio", + nu, + "Poisson's ratio", + Patterns::Double(-1.0, 0.5)); + + prm.add_parameter("rho", rho, "density", Patterns::Double(0.0)); + + prm.add_parameter("body forces", + body_force, + "body forces x,y,z", + Patterns::List(Patterns::Double())); } prm.leave_subsection(); @@ -79,158 +52,115 @@ namespace Parameters void - Solver::declare_parameters(ParameterHandler &prm) + Solver::add_output_parameters(ParameterHandler &prm) { - prm.enter_subsection("Linear solver"); + prm.enter_subsection("Solver"); { - prm.declare_entry("Solver type", - "Direct", - Patterns::Selection("CG|Direct"), - "Linear solver: CG or Direct"); + prm.add_parameter("Solver type", + type_lin, + "Linear solver: CG or Direct", + Patterns::Selection("CG|Direct")); - prm.declare_entry("Residual", - "1e-6", - Patterns::Double(0.0), - "Linear solver residual (scaled by residual norm)"); + prm.add_parameter("Residual", + tol_lin, + "Linear solver residual (scaled by residual norm)", + Patterns::Double(0.0)); - prm.declare_entry( + prm.add_parameter( "Max iteration multiplier", - "1", - Patterns::Double(0.0), - "Linear solver iterations (multiples of the system matrix size)"); - - prm.declare_entry("Max iterations Newton-Raphson", - "10", - Patterns::Integer(0), - "Number of Newton-Raphson iterations allowed"); - - prm.declare_entry("Tolerance force", - "1.0e-9", - Patterns::Double(0.0), - "Force residual tolerance"); - - prm.declare_entry("Tolerance displacement", - "1.0e-6", - Patterns::Double(0.0), - "Displacement error tolerance"); - } - prm.leave_subsection(); - } - - - void - Solver::parse_parameters(ParameterHandler &prm) - { - prm.enter_subsection("Solver"); - { - type_lin = prm.get("Solver type"); - tol_lin = prm.get_double("Residual"); - max_iterations_lin = prm.get_double("Max iteration multiplier"); - - max_iterations_NR = prm.get_integer("Max iterations Newton-Raphson"); - tol_f = prm.get_double("Tolerance force"); - tol_u = prm.get_double("Tolerance displacement"); + max_iterations_lin, + "Linear solver iterations (multiples of the system matrix size)", + Patterns::Double(0.0)); + + prm.add_parameter("Max iterations Newton-Raphson", + max_iterations_NR, + "Number of Newton-Raphson iterations allowed", + Patterns::Integer(0)); + + prm.add_parameter("Tolerance force", + tol_f, + "Force residual tolerance", + Patterns::Double(0.0)); + + prm.add_parameter("Tolerance displacement", + tol_u, + "Displacement error tolerance", + Patterns::Double(0.0)); } prm.leave_subsection(); } void - Discretization::declare_parameters(ParameterHandler &prm) + Discretization::add_output_parameters(ParameterHandler &prm) { prm.enter_subsection("Discretization"); { - prm.declare_entry("Polynomial degree", - "3", - Patterns::Integer(0), - "Polynomial degree of the FE system"); - - prm.declare_entry("theta", - "0.5", - Patterns::Double(0, 1), - "Time integration scheme"); - - prm.declare_entry("beta", - "0.25", - Patterns::Double(0, 0.5), - "Newmark beta"); - - prm.declare_entry("gamma", - "0.5", - Patterns::Double(0, 1), - "Newmark gamma"); - } - prm.leave_subsection(); - } - - void - Discretization::parse_parameters(ParameterHandler &prm) - { - prm.enter_subsection("Discretization"); - { - poly_degree = prm.get_integer("Polynomial degree"); - theta = prm.get_double("theta"); - beta = prm.get_double("beta"); - gamma = prm.get_double("gamma"); + prm.add_parameter("Polynomial degree", + poly_degree, + "Polynomial degree of the FE system", + Patterns::Integer(0)); + + prm.add_parameter("theta", + theta, + "Time integration scheme", + Patterns::Double(0, 1)); + + prm.add_parameter("beta", beta, "Newmark beta", Patterns::Double(0, 0.5)); + + prm.add_parameter("gamma", + gamma, + "Newmark gamma", + Patterns::Double(0, 1)); } prm.leave_subsection(); } void - PreciceAdapterConfiguration::declare_parameters(ParameterHandler &prm) + PreciceAdapterConfiguration::add_output_parameters(ParameterHandler &prm) { prm.enter_subsection("precice configuration"); { - prm.declare_entry("Scenario", - "FSI3", - Patterns::Selection("FSI3|PF"), - "Cases: FSI3 or PF for perpendicular flap"); - prm.declare_entry("precice config-file", - "precice-config.xml", - Patterns::Anything(), - "Name of the precice configuration file"); - prm.declare_entry( + prm.add_parameter("Scenario", + scenario, + "Cases: FSI3 or PF for perpendicular flap", + Patterns::Selection("FSI3|PF")); + + prm.add_parameter("precice config-file", + config_file, + "Name of the precice configuration file", + Patterns::Anything()); + + prm.add_parameter( "Participant name", - "dealiisolver", - Patterns::Anything(), - "Name of the participant in the precice-config.xml file"); - prm.declare_entry( + participant_name, + "Name of the participant in the precice-config.xml file", + Patterns::Anything()); + + prm.add_parameter( "Mesh name", - "dealii-mesh", - Patterns::Anything(), - "Name of the coupling mesh in the precice-config.xml file"); - prm.declare_entry("Read data name", - "received-data", - Patterns::Anything(), - "Name of the read data in the precice-config.xml file"); - prm.declare_entry( - "Write data name", - "calculated-data", - Patterns::Anything(), - "Name of the write data in the precice-config.xml file"); - prm.declare_entry("Flap location", - "0.0", - Patterns::Double(-3, 3), - "PF x-location"); + mesh_name, + "Name of the coupling mesh in the precice-config.xml file", + Patterns::Anything()); + + prm.add_parameter("Read data name", + read_data_name, + "Name of the read data in the precice-config.xml file", + Patterns::Anything()); + + prm.add_parameter("Write data name", + write_data_name, + "Name of the write data in the precice-config.xml file", + Patterns::Anything()); + + prm.add_parameter("Flap location", + flap_location, + "PF x-location", + Patterns::Double(-3, 3)); } prm.leave_subsection(); - } - void - PreciceAdapterConfiguration::parse_parameters(ParameterHandler &prm) - { - prm.enter_subsection("precice configuration"); - { - scenario = prm.get("Scenario"); - config_file = prm.get("precice config-file"); - participant_name = prm.get("Participant name"); - mesh_name = prm.get("Mesh name"); - read_data_name = prm.get("Read data name"); - write_data_name = prm.get("Write data name"); - flap_location = prm.get_double("Flap location"); - } - prm.leave_subsection(); // Look at the specific type of read data if ((read_data_name.find("Stress") == 0)) data_consistent = true; @@ -243,35 +173,21 @@ namespace Parameters "Unknown read data type. Please use 'Force' or 'Stress' in the read data naming.")); } + AllParameters::AllParameters(const std::string &input_file) { ParameterHandler prm; - declare_parameters(prm); + + Solver::add_output_parameters(prm); + Discretization::add_output_parameters(prm); + System::add_output_parameters(prm); + Time::add_output_parameters(prm); + PreciceAdapterConfiguration::add_output_parameters(prm); + prm.parse_input(input_file); - parse_parameters(prm); // Optional, if we want to print all parameters in the beginning of the // simulation // prm.print_parameters(std::cout,ParameterHandler::Text); } - - void - AllParameters::declare_parameters(ParameterHandler &prm) - { - Solver::declare_parameters(prm); - Discretization::declare_parameters(prm); - System::declare_parameters(prm); - Time::declare_parameters(prm); - PreciceAdapterConfiguration::declare_parameters(prm); - } - - void - AllParameters::parse_parameters(ParameterHandler &prm) - { - Solver::parse_parameters(prm); - Discretization::parse_parameters(prm); - System::parse_parameters(prm); - Time::parse_parameters(prm); - PreciceAdapterConfiguration::parse_parameters(prm); - } } // namespace Parameters diff --git a/include/adapter/parameters.h b/include/adapter/parameters.h index bdad65d..862c148 100644 --- a/include/adapter/parameters.h +++ b/include/adapter/parameters.h @@ -16,15 +16,12 @@ namespace Parameters */ struct Time { - double end_time; - double delta_t; - int output_interval; - - static void - declare_parameters(ParameterHandler &prm); + double end_time = 1; + double delta_t = 0.1; + int output_interval = 1; void - parse_parameters(ParameterHandler &prm); + add_output_parameters(ParameterHandler &prm); }; /** @@ -32,18 +29,14 @@ namespace Parameters */ struct System { - double nu; - double mu; - double lambda; - double rho; + double nu = 0.3; + double mu = 1538462; + double lambda = 2307692; + double rho = 1000; Tensor<1, 3, double> body_force; - - static void - declare_parameters(ParameterHandler &prm); - void - parse_parameters(ParameterHandler &prm); + add_output_parameters(ParameterHandler &prm); }; @@ -52,18 +45,15 @@ namespace Parameters */ struct Solver { - std::string type_lin; - double tol_lin; - double max_iterations_lin; - unsigned int max_iterations_NR; - double tol_f; - double tol_u; - - static void - declare_parameters(ParameterHandler &prm); + std::string type_lin = "Direct"; + double tol_lin = 1e-6; + double max_iterations_lin = 1; + unsigned int max_iterations_NR = 10; + double tol_f = 1e-9; + double tol_u = 1e-6; void - parse_parameters(ParameterHandler &prm); + add_output_parameters(ParameterHandler &prm); }; @@ -74,18 +64,15 @@ namespace Parameters */ struct Discretization { - unsigned int poly_degree; + unsigned int poly_degree = 3; // For the linear elastic model (theta-scheme) - double theta; + double theta = 0.5; // For the nonlinear elastic model (Newmark) - double beta; - double gamma ; - - static void - declare_parameters(ParameterHandler &prm); + double beta = 0.25; + double gamma = 0.5; void - parse_parameters(ParameterHandler &prm); + add_output_parameters(ParameterHandler &prm); }; @@ -96,20 +83,17 @@ namespace Parameters */ struct PreciceAdapterConfiguration { - std::string scenario; - std::string config_file; - std::string participant_name; - std::string mesh_name; - std::string read_data_name; - std::string write_data_name; - double flap_location; - bool data_consistent; - - static void - declare_parameters(ParameterHandler &prm); + std::string scenario = "FSI3"; + std::string config_file = "precice-config.xml"; + std::string participant_name = "dealiisolver"; + std::string mesh_name = "dealii-mesh"; + std::string read_data_name = "Stress"; + std::string write_data_name = "Displacement"; + double flap_location = 0.0; + bool data_consistent = true; void - parse_parameters(ParameterHandler &prm); + add_output_parameters(ParameterHandler &prm); }; @@ -121,12 +105,6 @@ namespace Parameters { AllParameters(const std::string &input_file); - - static void - declare_parameters(ParameterHandler &prm); - - void - parse_parameters(ParameterHandler &prm); }; } // namespace Parameters From bb7ed54b5101508d011338168aad05a3b7952f9c Mon Sep 17 00:00:00 2001 From: DavidSCN Date: Wed, 3 Mar 2021 10:52:49 +0100 Subject: [PATCH 11/18] Remove specific parameter files and add unified one --- dealii-precice.cc | 2 +- .../linear_elasticity.prm => parameters.prm | 29 +++++++++++-------- 2 files changed, 18 insertions(+), 13 deletions(-) rename source/linear_elasticity/linear_elasticity.prm => parameters.prm (76%) diff --git a/dealii-precice.cc b/dealii-precice.cc index e9ae8a4..fb45460 100644 --- a/dealii-precice.cc +++ b/dealii-precice.cc @@ -45,7 +45,7 @@ main(int argc, char **argv) if (argc > 1) parameter_file = argv[2]; else - parameter_file = "nonlinear_elasticity.prm"; + parameter_file = "parameters.prm"; // Extract case path for the output directory size_t pos = parameter_file.find_last_of("/"); diff --git a/source/linear_elasticity/linear_elasticity.prm b/parameters.prm similarity index 76% rename from source/linear_elasticity/linear_elasticity.prm rename to parameters.prm index 62edcbe..9cf3410 100644 --- a/source/linear_elasticity/linear_elasticity.prm +++ b/parameters.prm @@ -14,29 +14,25 @@ subsection Time end subsection Discretization - # Time integration scheme - # 0 = forward, 1 = backward - set theta = 0.5 - # Polynomial degree of the FE system set Polynomial degree = 3 end subsection System properties - # mu (shear modulus) - set mu = 0.5e6 + # Poisson's ratio + set Poisson's ratio = 0.4 - # lambda - set lambda = 2e6 + # Shear modulus + set Shear modulus = 0.5e6 - # density - set rho = 1000 + # Density + set rho = 1000 - # body forces x,y,z + # Body forces x,y,z set body forces = 0.0,0.0,0.0 end -subsection Linear solver +subsection Solver # Linear solver iterations (multiples of the system matrix size) set Max iteration multiplier = 1 @@ -45,6 +41,15 @@ subsection Linear solver # Linear solver: CG or Direct set Solver type = Direct + + # Number of Newton-Raphson iterations allowed + set Max iterations Newton-Raphson = 10 + + # Displacement error tolerance + set Tolerance displacement = 1.0e-6 + + # Force residual tolerance + set Tolerance force = 1.0e-9 end subsection precice configuration From d691fa22a5efa6c4812733d1e3308fdf9fb045c1 Mon Sep 17 00:00:00 2001 From: DavidSCN Date: Wed, 3 Mar 2021 12:02:47 +0100 Subject: [PATCH 12/18] Add parameter option for model selection --- dealii-precice.cc | 24 +++++++++++++----------- include/adapter/parameters.cc | 5 +++++ include/adapter/parameters.h | 1 + parameters.prm | 3 +++ 4 files changed, 22 insertions(+), 11 deletions(-) diff --git a/dealii-precice.cc b/dealii-precice.cc index fb45460..7749303 100644 --- a/dealii-precice.cc +++ b/dealii-precice.cc @@ -1,3 +1,5 @@ +#include + #include "source/linear_elasticity/include/linear_elasticity.h" #include "source/nonlinear_elasticity/include/nonlinear_elasticity.h" @@ -40,32 +42,32 @@ main(int argc, char **argv) << std::endl << std::endl; - - std::string parameter_file; - if (argc > 1) - parameter_file = argv[2]; - else - parameter_file = "parameters.prm"; + // Store the name of the parameter file + const std::string parameter_file = argc > 1 ? argv[1] : "parameters.prm"; // Extract case path for the output directory size_t pos = parameter_file.find_last_of("/"); std::string case_path = std::string::npos == pos ? "" : parameter_file.substr(0, pos + 1); - std::string solver_type = argv[1]; - // Dimension is determinded via cmake -DDIM - if (solver_type == "-nonlinear") // nonlinear + // Query solver type from the parameter file + ParameterHandler prm; + Parameters::Solver solver; + solver.add_output_parameters(prm); + prm.parse_input(case_path + parameter_file, "", true); + + if (solver.model == "neo-Hook") // nonlinear { Nonlinear_Elasticity::Solid solid(case_path); solid.run(); } - else if (solver_type == "-linear") // linear + else if (solver.model == "linear") // linear { Linear_Elasticity::ElastoDynamics elastic_solver(case_path); elastic_solver.run(); } else - AssertThrow(false, ExcMessage("Unknown solver specified. ")) + AssertThrow(false, ExcNotImplemented()) } catch (std::exception &exc) { diff --git a/include/adapter/parameters.cc b/include/adapter/parameters.cc index 97b24fa..088591d 100644 --- a/include/adapter/parameters.cc +++ b/include/adapter/parameters.cc @@ -56,6 +56,11 @@ namespace Parameters { prm.enter_subsection("Solver"); { + prm.add_parameter("Model", + model, + "Solid model to be used: linear or neo-Hook", + Patterns::Selection("linear|neo-Hook")); + prm.add_parameter("Solver type", type_lin, "Linear solver: CG or Direct", diff --git a/include/adapter/parameters.h b/include/adapter/parameters.h index 862c148..fc6c6f8 100644 --- a/include/adapter/parameters.h +++ b/include/adapter/parameters.h @@ -45,6 +45,7 @@ namespace Parameters */ struct Solver { + std::string model = "linear"; std::string type_lin = "Direct"; double tol_lin = 1e-6; double max_iterations_lin = 1; diff --git a/parameters.prm b/parameters.prm index 9cf3410..bcaf5e6 100644 --- a/parameters.prm +++ b/parameters.prm @@ -33,6 +33,9 @@ subsection System properties end subsection Solver + # Structural model to be used: linear or neo-Hook + set Model = linear + # Linear solver iterations (multiples of the system matrix size) set Max iteration multiplier = 1 From 8997e9fbca474d9b1740a278b62290947310d13f Mon Sep 17 00:00:00 2001 From: DavidSCN Date: Wed, 3 Mar 2021 12:11:33 +0100 Subject: [PATCH 13/18] Add an entry for the parameter file location --- dealii-precice.cc | 7 ++++--- include/adapter/parameters.cc | 2 +- source/linear_elasticity/include/linear_elasticity.h | 3 ++- source/linear_elasticity/linear_elasticity.cc | 5 +++-- source/nonlinear_elasticity/include/nonlinear_elasticity.h | 2 +- source/nonlinear_elasticity/nonlinear_elasticity.cc | 6 +++--- 6 files changed, 14 insertions(+), 11 deletions(-) diff --git a/dealii-precice.cc b/dealii-precice.cc index 7749303..d646f7e 100644 --- a/dealii-precice.cc +++ b/dealii-precice.cc @@ -54,16 +54,17 @@ main(int argc, char **argv) ParameterHandler prm; Parameters::Solver solver; solver.add_output_parameters(prm); - prm.parse_input(case_path + parameter_file, "", true); + prm.parse_input(parameter_file, "", true); if (solver.model == "neo-Hook") // nonlinear { - Nonlinear_Elasticity::Solid solid(case_path); + Nonlinear_Elasticity::Solid solid(case_path, parameter_file); solid.run(); } else if (solver.model == "linear") // linear { - Linear_Elasticity::ElastoDynamics elastic_solver(case_path); + Linear_Elasticity::ElastoDynamics elastic_solver(case_path, + parameter_file); elastic_solver.run(); } else diff --git a/include/adapter/parameters.cc b/include/adapter/parameters.cc index 088591d..37b52f3 100644 --- a/include/adapter/parameters.cc +++ b/include/adapter/parameters.cc @@ -58,7 +58,7 @@ namespace Parameters { prm.add_parameter("Model", model, - "Solid model to be used: linear or neo-Hook", + "Structural model to be used: linear or neo-Hook", Patterns::Selection("linear|neo-Hook")); prm.add_parameter("Solver type", diff --git a/source/linear_elasticity/include/linear_elasticity.h b/source/linear_elasticity/include/linear_elasticity.h index 78ac2b0..ae270fa 100644 --- a/source/linear_elasticity/include/linear_elasticity.h +++ b/source/linear_elasticity/include/linear_elasticity.h @@ -55,7 +55,8 @@ namespace Linear_Elasticity class ElastoDynamics { public: - ElastoDynamics(const std::string &case_path); + ElastoDynamics(const std::string &case_path, + const std::string ¶meter_file); ~ElastoDynamics(); // As usual in dealii, the run function covers the main time loop of the diff --git a/source/linear_elasticity/linear_elasticity.cc b/source/linear_elasticity/linear_elasticity.cc index 0120e1d..cca4a1a 100644 --- a/source/linear_elasticity/linear_elasticity.cc +++ b/source/linear_elasticity/linear_elasticity.cc @@ -52,8 +52,9 @@ namespace Linear_Elasticity // Constructor template - ElastoDynamics::ElastoDynamics(const std::string &case_path) - : parameters(case_path + "linear_elasticity.prm") + ElastoDynamics::ElastoDynamics(const std::string &case_path, + const std::string ¶meter_file) + : parameters(parameter_file) , interface_boundary_id(6) , dof_handler(triangulation) , fe(FE_Q(parameters.poly_degree), dim) diff --git a/source/nonlinear_elasticity/include/nonlinear_elasticity.h b/source/nonlinear_elasticity/include/nonlinear_elasticity.h index a4a1270..b354c8f 100644 --- a/source/nonlinear_elasticity/include/nonlinear_elasticity.h +++ b/source/nonlinear_elasticity/include/nonlinear_elasticity.h @@ -132,7 +132,7 @@ namespace Nonlinear_Elasticity class Solid { public: - Solid(const std::string &case_path); + Solid(const std::string &case_path, const std::string ¶meter_file); virtual ~Solid(); diff --git a/source/nonlinear_elasticity/nonlinear_elasticity.cc b/source/nonlinear_elasticity/nonlinear_elasticity.cc index 19ea161..46abc3b 100644 --- a/source/nonlinear_elasticity/nonlinear_elasticity.cc +++ b/source/nonlinear_elasticity/nonlinear_elasticity.cc @@ -60,9 +60,9 @@ namespace Nonlinear_Elasticity // Constructor initializes member variables and reads the parameter file template - Solid::Solid(const std::string &case_path) - : parameters( - Parameters::AllParameters(case_path + "nonlinear_elasticity.prm")) + Solid::Solid(const std::string &case_path, + const std::string ¶meter_file) + : parameters(Parameters::AllParameters(parameter_file)) , vol_reference(0.0) , vol_current(0.0) , triangulation(Triangulation::maximum_smoothing) From 5476c4925dcc1cff28ff5ca389ffa05951018f32 Mon Sep 17 00:00:00 2001 From: DavidSCN Date: Wed, 3 Mar 2021 12:15:56 +0100 Subject: [PATCH 14/18] Fix Hook Hooke typo --- dealii-precice.cc | 2 +- include/adapter/parameters.cc | 4 ++-- parameters.prm | 2 +- 3 files changed, 4 insertions(+), 4 deletions(-) diff --git a/dealii-precice.cc b/dealii-precice.cc index d646f7e..5cf6e63 100644 --- a/dealii-precice.cc +++ b/dealii-precice.cc @@ -56,7 +56,7 @@ main(int argc, char **argv) solver.add_output_parameters(prm); prm.parse_input(parameter_file, "", true); - if (solver.model == "neo-Hook") // nonlinear + if (solver.model == "neo-Hooke") // nonlinear { Nonlinear_Elasticity::Solid solid(case_path, parameter_file); solid.run(); diff --git a/include/adapter/parameters.cc b/include/adapter/parameters.cc index 37b52f3..6a1c06e 100644 --- a/include/adapter/parameters.cc +++ b/include/adapter/parameters.cc @@ -58,8 +58,8 @@ namespace Parameters { prm.add_parameter("Model", model, - "Structural model to be used: linear or neo-Hook", - Patterns::Selection("linear|neo-Hook")); + "Structural model to be used: linear or neo-Hooke", + Patterns::Selection("linear|neo-Hooke")); prm.add_parameter("Solver type", type_lin, diff --git a/parameters.prm b/parameters.prm index bcaf5e6..1396b38 100644 --- a/parameters.prm +++ b/parameters.prm @@ -33,7 +33,7 @@ subsection System properties end subsection Solver - # Structural model to be used: linear or neo-Hook + # Structural model to be used: linear or neo-Hooke set Model = linear # Linear solver iterations (multiples of the system matrix size) From 81363a6649c1f8fb95c5a093725c14497ce58d91 Mon Sep 17 00:00:00 2001 From: DavidSCN Date: Sun, 21 Mar 2021 14:16:30 +0100 Subject: [PATCH 15/18] Apply suggestions from code review --- CMakeLists.txt | 2 +- README.md | 6 +-- dealii-precice.cc => elasticity.cc | 2 +- include/adapter/parameters.cc | 46 +++++++++++-------- parameters.prm | 21 +++++---- source/linear_elasticity/linear_elasticity.cc | 2 +- 6 files changed, 44 insertions(+), 35 deletions(-) rename dealii-precice.cc => elasticity.cc (98%) diff --git a/CMakeLists.txt b/CMakeLists.txt index fd8f940..0cbbdbd 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -18,7 +18,7 @@ ENDIF() ADD_DEFINITIONS(-DDIM=${DIM}) # Set the target and the target source -SET( TARGET "dealii-precice" ) +SET( TARGET "elasticity" ) SET( TARGET_SRC ${TARGET}.cc include/adapter/parameters.cc) DEAL_II_INITIALIZE_CACHED_VARIABLES() diff --git a/README.md b/README.md index a559ae7..64c320b 100644 --- a/README.md +++ b/README.md @@ -6,10 +6,10 @@ Coupled structural solvers written with the C++ finite element library deal.II: -- `linear_elasticity` contains a linear-elastic solver based on the [step-8 tutorial program](https://www.dealii.org/9.2.0/doxygen/deal.II/step_8.html) of deal.II -- `nonlinear_elasticity` contains a nonlinear elastic solver, which builds on previous work of Jean-Paul Pelteret and Andrew McBride in their deal.II code gallery program '[Quasi-Static Finite-Strain Compressible Elasticity](https://www.dealii.org/9.0.0/doxygen/deal.II/code_gallery_Quasi_static_Finite_strain_Compressible_Elasticity.html).' This solver supports shared-memory parallelization. +- `source/linear_elasticity` contains a linear-elastic solver based on the [step-8 tutorial program](https://www.dealii.org/9.2.0/doxygen/deal.II/step_8.html) of deal.II +- `source/nonlinear_elasticity` contains a nonlinear elastic solver, which builds on previous work of Jean-Paul Pelteret and Andrew McBride in their deal.II code gallery program '[Quasi-Static Finite-Strain Compressible Elasticity](https://www.dealii.org/9.0.0/doxygen/deal.II/code_gallery_Quasi_static_Finite_strain_Compressible_Elasticity.html).' This solver supports shared-memory parallelization. -Applied coupling functionalities have been separated and can be found in the `adapter` directory. +Applied coupling functionalities have been separated and can be found in the `include/adapter` directory. ## Start here Our [wiki](https://github.com/precice/dealii-adapter/wiki) will help you start. If you are missing something, [let us know](https://www.precice.org/resources/#contact). diff --git a/dealii-precice.cc b/elasticity.cc similarity index 98% rename from dealii-precice.cc rename to elasticity.cc index 5cf6e63..006cdef 100644 --- a/dealii-precice.cc +++ b/elasticity.cc @@ -56,7 +56,7 @@ main(int argc, char **argv) solver.add_output_parameters(prm); prm.parse_input(parameter_file, "", true); - if (solver.model == "neo-Hooke") // nonlinear + if (solver.model == "neo-Hookean") // nonlinear { Nonlinear_Elasticity::Solid solid(case_path, parameter_file); solid.run(); diff --git a/include/adapter/parameters.cc b/include/adapter/parameters.cc index 6a1c06e..bea72b0 100644 --- a/include/adapter/parameters.cc +++ b/include/adapter/parameters.cc @@ -58,39 +58,43 @@ namespace Parameters { prm.add_parameter("Model", model, - "Structural model to be used: linear or neo-Hooke", - Patterns::Selection("linear|neo-Hooke")); + "Structural model to be used: linear or neo-Hookean", + Patterns::Selection("linear|neo-Hookean")); prm.add_parameter("Solver type", type_lin, "Linear solver: CG or Direct", Patterns::Selection("CG|Direct")); - prm.add_parameter("Residual", - tol_lin, - "Linear solver residual (scaled by residual norm)", - Patterns::Double(0.0)); + prm.add_parameter( + "Residual", + tol_lin, + "CG solver residual (multiplied by residual norm, ignored if Model == linear)", + Patterns::Double(0.0)); prm.add_parameter( "Max iteration multiplier", max_iterations_lin, - "Linear solver iterations (multiples of the system matrix size)", + "Max CG solver iterations (multiples of the system matrix size)", Patterns::Double(0.0)); - prm.add_parameter("Max iterations Newton-Raphson", - max_iterations_NR, - "Number of Newton-Raphson iterations allowed", - Patterns::Integer(0)); + prm.add_parameter( + "Max iterations Newton-Raphson", + max_iterations_NR, + "Number of Newton-Raphson iterations allowed (ignored if Model == linear)", + Patterns::Integer(0)); - prm.add_parameter("Tolerance force", - tol_f, - "Force residual tolerance", - Patterns::Double(0.0)); + prm.add_parameter( + "Tolerance force", + tol_f, + "Force residual tolerance for non-linear iteration (ignored if Model == linear)", + Patterns::Double(0.0)); - prm.add_parameter("Tolerance displacement", - tol_u, - "Displacement error tolerance", - Patterns::Double(0.0)); + prm.add_parameter( + "Tolerance displacement", + tol_u, + "Displacement error tolerance for non-linear iteration (ignored if Model == linear)", + Patterns::Double(0.0)); } prm.leave_subsection(); } @@ -191,6 +195,10 @@ namespace Parameters prm.parse_input(input_file); + AssertThrow((data_consistent && model == "neo-Hookean") || + model == "linear", + ExcNotImplemented()); + // Optional, if we want to print all parameters in the beginning of the // simulation // prm.print_parameters(std::cout,ParameterHandler::Text); diff --git a/parameters.prm b/parameters.prm index 1396b38..3485920 100644 --- a/parameters.prm +++ b/parameters.prm @@ -33,25 +33,26 @@ subsection System properties end subsection Solver - # Structural model to be used: linear or neo-Hooke + # Structural model to be used: linear or neo-Hookean set Model = linear - # Linear solver iterations (multiples of the system matrix size) + # Linear solver: CG or Direct + set Solver type = Direct + + # Max CG solver iterations (multiples of the system matrix size) + # In 2-d, this value is best set at 2. In 3-d, a value of 1 work fine. set Max iteration multiplier = 1 - # Linear solver residual (scaled by residual norm) + # Absolute CG solver residual (multiplied by residual norm, ignored if Model == linear) set Residual = 1e-6 - # Linear solver: CG or Direct - set Solver type = Direct - - # Number of Newton-Raphson iterations allowed + # Number of Newton-Raphson iterations allowed (ignored if Model == linear) set Max iterations Newton-Raphson = 10 - # Displacement error tolerance + # Relative displacement error tolerance for non-linear iteration (ignored if Model == linear) set Tolerance displacement = 1.0e-6 - # Force residual tolerance + # Relative force residual tolerance for non-linear iteration (ignored if Model == linear) set Tolerance force = 1.0e-9 end @@ -66,7 +67,7 @@ subsection precice configuration set Participant name = Solid # Name of the coupling mesh in the precice-config.xml file - set Mesh name = Solid_mesh + set Mesh name = Solid-Mesh # Name of the read data in the precice-config.xml file set Read data name = Stress diff --git a/source/linear_elasticity/linear_elasticity.cc b/source/linear_elasticity/linear_elasticity.cc index cca4a1a..1c7dfc8 100644 --- a/source/linear_elasticity/linear_elasticity.cc +++ b/source/linear_elasticity/linear_elasticity.cc @@ -541,7 +541,7 @@ namespace Linear_Elasticity const int solver_its = system_matrix.m() * parameters.max_iterations_lin; - const double tol_sol = parameters.tol_lin * system_rhs.l2_norm(); + const double tol_sol = 1.e-10; SolverControl solver_control(solver_its, tol_sol); GrowingVectorMemory<> GVM; From 005535ac7962c8352ff61e5ab81cb310ba08bd85 Mon Sep 17 00:00:00 2001 From: DavidSCN Date: Sun, 21 Mar 2021 14:47:22 +0100 Subject: [PATCH 16/18] Fix 'Force' data check for non-linear model --- include/adapter/parameters.cc | 24 ++++++++----------- .../nonlinear_elasticity.cc | 8 ++++++- 2 files changed, 17 insertions(+), 15 deletions(-) diff --git a/include/adapter/parameters.cc b/include/adapter/parameters.cc index bea72b0..dc73797 100644 --- a/include/adapter/parameters.cc +++ b/include/adapter/parameters.cc @@ -169,17 +169,6 @@ namespace Parameters Patterns::Double(-3, 3)); } prm.leave_subsection(); - - // Look at the specific type of read data - if ((read_data_name.find("Stress") == 0)) - data_consistent = true; - else if ((read_data_name.find("Force") == 0)) - data_consistent = false; - else - AssertThrow( - false, - ExcMessage( - "Unknown read data type. Please use 'Force' or 'Stress' in the read data naming.")); } @@ -195,9 +184,16 @@ namespace Parameters prm.parse_input(input_file); - AssertThrow((data_consistent && model == "neo-Hookean") || - model == "linear", - ExcNotImplemented()); + // Look at the specific type of read data + if ((read_data_name.find("Stress") == 0)) + data_consistent = true; + else if ((read_data_name.find("Force") == 0)) + data_consistent = false; + else + AssertThrow( + false, + ExcMessage( + "Unknown read data type. Please use 'Force' or 'Stress' in the read data naming.")); // Optional, if we want to print all parameters in the beginning of the // simulation diff --git a/source/nonlinear_elasticity/nonlinear_elasticity.cc b/source/nonlinear_elasticity/nonlinear_elasticity.cc index 46abc3b..5821989 100644 --- a/source/nonlinear_elasticity/nonlinear_elasticity.cc +++ b/source/nonlinear_elasticity/nonlinear_elasticity.cc @@ -82,7 +82,13 @@ namespace Nonlinear_Elasticity , timer(std::cout, TimerOutput::summary, TimerOutput::wall_times) , time(parameters.end_time, parameters.delta_t) , adapter(parameters, boundary_interface_id) - {} + { + AssertThrow( + parameters.data_consistent, + ExcMessage( + "The neo-Hookean solid doesn't support 'Force' data reading. Please switch to 'Stress' " + "data on the Fluid side or use the linear model of the solid solver")); + } // Destructor clears the DoFHandler template From 54893cc98a1a9528a9b83077db1140426ebeaaa8 Mon Sep 17 00:00:00 2001 From: DavidSCN Date: Mon, 22 Mar 2021 08:30:39 +0100 Subject: [PATCH 17/18] Enfore compatibility to tutorial parameter file --- parameters.prm | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/parameters.prm b/parameters.prm index 3485920..c10c444 100644 --- a/parameters.prm +++ b/parameters.prm @@ -40,7 +40,7 @@ subsection Solver set Solver type = Direct # Max CG solver iterations (multiples of the system matrix size) - # In 2-d, this value is best set at 2. In 3-d, a value of 1 work fine. + # In 2D, this value is best set at 2. In 3D, a value of 1 works fine. set Max iteration multiplier = 1 # Absolute CG solver residual (multiplied by residual norm, ignored if Model == linear) From 0edace1f91e193fac2cd58b6c8005d164f39632b Mon Sep 17 00:00:00 2001 From: DavidSCN Date: Tue, 23 Mar 2021 15:29:10 +0100 Subject: [PATCH 18/18] Specify linear residual in parameter file more in detail --- parameters.prm | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/parameters.prm b/parameters.prm index c10c444..d767d85 100644 --- a/parameters.prm +++ b/parameters.prm @@ -43,7 +43,8 @@ subsection Solver # In 2D, this value is best set at 2. In 3D, a value of 1 works fine. set Max iteration multiplier = 1 - # Absolute CG solver residual (multiplied by residual norm, ignored if Model == linear) + # Relative drop criterion for CG solver residual (multiplied by residual norm, ignored if Model == linear or solver == direct) + # Hard-coded to absolute criterion of 1e-10 for the linear model. set Residual = 1e-6 # Number of Newton-Raphson iterations allowed (ignored if Model == linear)