diff --git a/Applications/ApplicationsLib/ProjectData.cpp b/Applications/ApplicationsLib/ProjectData.cpp index e5176f05673f6c9be2ba6e8ded82e8192ed08473..5787f506a2c01a08749586bd195635b1025efc34 100644 --- a/Applications/ApplicationsLib/ProjectData.cpp +++ b/Applications/ApplicationsLib/ProjectData.cpp @@ -33,6 +33,7 @@ #include "ProcessLib/GroundwaterFlow/CreateGroundwaterFlowProcess.h" #include "ProcessLib/TES/CreateTESProcess.h" +#include "ProcessLib/HeatConduction/CreateHeatConductionProcess.h" namespace detail { @@ -263,6 +264,11 @@ void ProjectData::parseProcesses(BaseLib::ConfigTree const& processes_config) process = ProcessLib::TES::createTESProcess( *_mesh_vec[0], _process_variables, _parameters, process_config); } + else if (type == "HEAT_CONDUCTION") + { + process = ProcessLib::HeatConduction::createHeatConductionProcess( + *_mesh_vec[0], _process_variables, _parameters, process_config); + } else { OGS_FATAL("Unknown process type: %s", type.c_str()); diff --git a/Applications/CLI/Tests.cmake b/Applications/CLI/Tests.cmake index 2c218b4c432b345f8c3f5e987711ba7c0ec5468e..4af4d051e1f6e3a379ca254ef8ece9dd91c15f6f 100644 --- a/Applications/CLI/Tests.cmake +++ b/Applications/CLI/Tests.cmake @@ -217,6 +217,32 @@ if(NOT OGS_USE_MPI) tes_zeolite_discharge_large_pcs_0_ts_28_t_1_000000.vtu tes_zeolite_discharge_large_pcs_0_ts_28_t_1.000000.vtu solid_density fct_solid_density ) + AddTest( + NAME 1D_HeatConduction_dirichlet + PATH Parabolic/T/1D_dirichlet + EXECUTABLE ogs + EXECUTABLE_ARGS line_60_heat.prj + WRAPPER time + TESTER vtkdiff + ABSTOL 1e-5 RELTOL 1e-5 + DIFF_DATA + temperature_analytical.vtu line_60_heat_pcs_0_ts_65_t_5078125.000000.vtu Temperature_Analytical_2months temperature + temperature_analytical.vtu line_60_heat_pcs_0_ts_405_t_31640625.000000.vtu Temperature_Analytical_1year temperature + ) + + AddTest( + NAME 1D_HeatConduction_neumann + PATH Parabolic/T/1D_neumann + EXECUTABLE ogs + EXECUTABLE_ARGS line_60_heat.prj + WRAPPER time + TESTER vtkdiff + ABSTOL 1e-4 RELTOL 1e-4 + DIFF_DATA + temperature_analytical.vtu line_60_heat_pcs_0_ts_65_t_5078125.000000.vtu Temperature_Analytical_2months temperature + temperature_analytical.vtu line_60_heat_pcs_0_ts_405_t_31640625.000000.vtu Temperature_Analytical_1year temperature + ) + else() # MPI groundwater flow tests AddTest( diff --git a/ProcessLib/CMakeLists.txt b/ProcessLib/CMakeLists.txt index d3883b0d7784c4fda37015eb1b3dd52308cd35f5..ebc1893a65ac7372204d865674626f414f80f24e 100644 --- a/ProcessLib/CMakeLists.txt +++ b/ProcessLib/CMakeLists.txt @@ -15,6 +15,9 @@ APPEND_SOURCE_FILES(SOURCES GroundwaterFlow) add_subdirectory(TES) APPEND_SOURCE_FILES(SOURCES TES) +add_subdirectory(HeatConduction) +APPEND_SOURCE_FILES(SOURCES HeatConduction) + APPEND_SOURCE_FILES(SOURCES Utils) add_library(ProcessLib ${SOURCES}) diff --git a/ProcessLib/HeatConduction/CMakeLists.txt b/ProcessLib/HeatConduction/CMakeLists.txt new file mode 100644 index 0000000000000000000000000000000000000000..e69de29bb2d1d6434b8b29ae775ad8c2e48c5391 diff --git a/ProcessLib/HeatConduction/CreateHeatConductionProcess.cpp b/ProcessLib/HeatConduction/CreateHeatConductionProcess.cpp new file mode 100644 index 0000000000000000000000000000000000000000..8241a0b33a1f0acdbdc0e39dec8aabf237cca0ef --- /dev/null +++ b/ProcessLib/HeatConduction/CreateHeatConductionProcess.cpp @@ -0,0 +1,80 @@ +/** + * \copyright + * Copyright (c) 2012-2016, OpenGeoSys Community (http://www.opengeosys.org) + * Distributed under a Modified BSD License. + * See accompanying file LICENSE.txt or + * http://www.opengeosys.org/project/license + * + */ + +#include "CreateHeatConductionProcess.h" + +#include "ProcessLib/Utils/ParseSecondaryVariables.h" +#include "ProcessLib/Utils/ProcessUtils.h" +#include "HeatConductionProcess.h" +#include "HeatConductionProcessData.h" + +namespace ProcessLib +{ +namespace HeatConduction +{ +std::unique_ptr<Process> createHeatConductionProcess( + MeshLib::Mesh& mesh, + std::vector<ProcessVariable> const& variables, + std::vector<std::unique_ptr<ParameterBase>> const& parameters, + BaseLib::ConfigTree const& config) +{ + //! \ogs_file_param{process__type} + config.checkConfigParameter("type", "HEAT_CONDUCTION"); + + DBUG("Create HeatConductionProcess."); + + // Process variable. + auto process_variables = findProcessVariables( + variables, config, + {//! \ogs_file_param_special{process__HEAT_CONDUCTION__process_variables__process_variable} + "process_variable"}); + + // thermal conductivity parameter. + auto& thermal_conductivity = findParameter<double>( + config, + //! \ogs_file_param_special{process__HEAT_CONDUCTION__thermal_conductivity} + "thermal_conductivity", parameters, 1); + + DBUG("Use \'%s\' as thermal conductivity parameter.", + thermal_conductivity.name.c_str()); + + // heat capacity parameter. + auto& heat_capacity = findParameter<double>( + config, + //! \ogs_file_param_special{process__HEAT_CONDUCTION__heat_capacity} + "heat_capacity", parameters, 1); + + DBUG("Use \'%s\' as heat capacity parameter.", heat_capacity.name.c_str()); + + // density parameter. + auto& density = findParameter<double>( + config, + //! \ogs_file_param_special{process__HEAT_CONDUCTION__density} + "density", parameters, 1); + + DBUG("Use \'%s\' as density parameter.", density.name.c_str()); + + HeatConductionProcessData process_data{thermal_conductivity, heat_capacity, + density}; + + SecondaryVariableCollection secondary_variables; + + NumLib::NamedFunctionCaller named_function_caller( + {"HeatConduction_temperature"}); + + ProcessLib::parseSecondaryVariables(config, secondary_variables, + named_function_caller); + + return std::unique_ptr<Process>{new HeatConductionProcess{ + mesh, parameters, std::move(process_variables), std::move(process_data), + std::move(secondary_variables), std::move(named_function_caller)}}; +} + +} // namespace HeatConduction +} // namespace ProcessLib diff --git a/ProcessLib/HeatConduction/CreateHeatConductionProcess.h b/ProcessLib/HeatConduction/CreateHeatConductionProcess.h new file mode 100644 index 0000000000000000000000000000000000000000..43979abaa7733645de5e06d4f70d0c0ef45ec0dc --- /dev/null +++ b/ProcessLib/HeatConduction/CreateHeatConductionProcess.h @@ -0,0 +1,29 @@ +/** + * \copyright + * Copyright (c) 2012-2016, OpenGeoSys Community (http://www.opengeosys.org) + * Distributed under a Modified BSD License. + * See accompanying file LICENSE.txt or + * http://www.opengeosys.org/project/license + * + */ + +#ifndef PROCESS_LIB_CREATE_HEATCONDUCTIONPROCESS_H_ +#define PROCESS_LIB_CREATE_HEATCONDUCTIONPROCESS_H_ + +#include <memory> +#include "ProcessLib/Process.h" + +namespace ProcessLib +{ +namespace HeatConduction +{ +std::unique_ptr<Process> createHeatConductionProcess( + MeshLib::Mesh& mesh, + std::vector<ProcessVariable> const& variables, + std::vector<std::unique_ptr<ParameterBase>> const& parameters, + BaseLib::ConfigTree const& config); + +} // namespace HeatConduction +} // namespace ProcessLib + +#endif // PROCESS_LIB_CREATE_HEATCONDUCTIONPROCESS_H_ diff --git a/ProcessLib/HeatConduction/HeatConductionFEM.h b/ProcessLib/HeatConduction/HeatConductionFEM.h new file mode 100644 index 0000000000000000000000000000000000000000..5289d59e9685dcdd9ee244f9d3213037c278eadc --- /dev/null +++ b/ProcessLib/HeatConduction/HeatConductionFEM.h @@ -0,0 +1,175 @@ +/** + * \copyright + * Copyright (c) 2012-2016, OpenGeoSys Community (http://www.opengeosys.org) + * Distributed under a Modified BSD License. + * See accompanying file LICENSE.txt or + * http://www.opengeosys.org/project/license + * + */ + +#ifndef PROCESS_LIB_HEATCONDUCTION_FEM_H_ +#define PROCESS_LIB_HEATCONDUCTION_FEM_H_ + +#include <vector> + +#include "NumLib/Extrapolation/ExtrapolatableElement.h" +#include "NumLib/Fem/FiniteElement/TemplateIsoparametric.h" +#include "NumLib/Fem/ShapeMatrixPolicy.h" +#include "ProcessLib/LocalAssemblerInterface.h" +#include "ProcessLib/LocalAssemblerTraits.h" +#include "ProcessLib/Parameter/Parameter.h" +#include "ProcessLib/Utils/InitShapeMatrices.h" +#include "HeatConductionProcessData.h" + +namespace ProcessLib +{ +namespace HeatConduction +{ +const unsigned NUM_NODAL_DOF = 1; + +class HeatConductionLocalAssemblerInterface + : public ProcessLib::LocalAssemblerInterface, + public NumLib::ExtrapolatableElement +{ +public: + virtual std::vector<double> const& getIntPtHeatFluxX( + std::vector<double>& /*cache*/) const = 0; + + virtual std::vector<double> const& getIntPtHeatFluxY( + std::vector<double>& /*cache*/) const = 0; + + virtual std::vector<double> const& getIntPtHeatFluxZ( + std::vector<double>& /*cache*/) const = 0; +}; + +template <typename ShapeFunction, typename IntegrationMethod, + unsigned GlobalDim> +class LocalAssemblerData : public HeatConductionLocalAssemblerInterface +{ + using ShapeMatricesType = ShapeMatrixPolicyType<ShapeFunction, GlobalDim>; + using ShapeMatrices = typename ShapeMatricesType::ShapeMatrices; + + using LocalAssemblerTraits = ProcessLib::LocalAssemblerTraits< + ShapeMatricesType, ShapeFunction::NPOINTS, NUM_NODAL_DOF, GlobalDim>; + + using NodalMatrixType = typename LocalAssemblerTraits::LocalMatrix; + using NodalVectorType = typename LocalAssemblerTraits::LocalVector; + using GlobalDimVectorType = typename ShapeMatricesType::GlobalDimVectorType; + +public: + /// The thermal_conductivity factor is directly integrated into the local + /// element matrix. + LocalAssemblerData(MeshLib::Element const& element, + std::size_t const local_matrix_size, + unsigned const integration_order, + HeatConductionProcessData const& process_data) + : _element(element), + _process_data(process_data), + _localK(local_matrix_size, + local_matrix_size), // TODO narrowing conversion + _localM(local_matrix_size, local_matrix_size), + _localRhs(local_matrix_size), + _integration_method(integration_order), + _shape_matrices(initShapeMatrices<ShapeFunction, ShapeMatricesType, + IntegrationMethod, GlobalDim>( + element, _integration_method)) + { + // This assertion is valid only if all nodal d.o.f. use the same shape + // matrices. + assert(local_matrix_size == ShapeFunction::NPOINTS * NUM_NODAL_DOF); + } + + void assembleConcrete( + double const t, std::vector<double> const& local_x, + NumLib::LocalToGlobalIndexMap::RowColumnIndices const& indices, + GlobalMatrix& M, GlobalMatrix& K, GlobalVector& b) override + { + _localK.setZero(); + _localM.setZero(); + _localRhs.setZero(); + + unsigned const n_integration_points = + _integration_method.getNumberOfPoints(); + + SpatialPosition pos; + pos.setElementID(_element.getID()); + + for (unsigned ip = 0; ip < n_integration_points; ip++) + { + pos.setIntegrationPoint(ip); + auto const& sm = _shape_matrices[ip]; + auto const& wp = _integration_method.getWeightedPoint(ip); + auto const k = _process_data.thermal_conductivity(t, pos)[0]; + auto const heat_capacity = _process_data.heat_capacity(t, pos)[0]; + auto const density = _process_data.density(t, pos)[0]; + + _localK.noalias() += + sm.dNdx.transpose() * k * sm.dNdx * sm.detJ * wp.getWeight(); + _localM.noalias() += sm.N.transpose() * density * heat_capacity * + sm.N * sm.detJ * wp.getWeight(); + // heat flux only computed for output. + GlobalDimVectorType const heat_flux = + -k * sm.dNdx * Eigen::Map<const NodalVectorType>( + local_x.data(), ShapeFunction::NPOINTS); + + for (unsigned d = 0; d < GlobalDim; ++d) + { + _heat_fluxes[d][ip] = heat_flux[d]; + } + } + + K.add(indices, _localK); + M.add(indices, _localM); + b.add(indices.rows, _localRhs); + } + + Eigen::Map<const Eigen::RowVectorXd> getShapeMatrix( + const unsigned integration_point) const override + { + auto const& N = _shape_matrices[integration_point].N; + + // assumes N is stored contiguously in memory + return Eigen::Map<const Eigen::RowVectorXd>(N.data(), N.size()); + } + + std::vector<double> const& getIntPtHeatFluxX( + std::vector<double>& /*cache*/) const override + { + assert(_heat_fluxes.size() > 0); + return _heat_fluxes[0]; + } + + std::vector<double> const& getIntPtHeatFluxY( + std::vector<double>& /*cache*/) const override + { + assert(_heat_fluxes.size() > 1); + return _heat_fluxes[1]; + } + + std::vector<double> const& getIntPtHeatFluxZ( + std::vector<double>& /*cache*/) const override + { + assert(_heat_fluxes.size() > 2); + return _heat_fluxes[2]; + } + +private: + MeshLib::Element const& _element; + HeatConductionProcessData const& _process_data; + + NodalMatrixType _localK; + NodalMatrixType _localM; + NodalVectorType _localRhs; + + IntegrationMethod const _integration_method; + std::vector<ShapeMatrices> _shape_matrices; + + std::vector<std::vector<double>> _heat_fluxes = + std::vector<std::vector<double>>( + GlobalDim, std::vector<double>(ShapeFunction::NPOINTS)); +}; + +} // namespace HeatConduction +} // namespace ProcessLib + +#endif // PROCESS_LIB_HEATCONDUCTION_FEM_H_ diff --git a/ProcessLib/HeatConduction/HeatConductionProcess.cpp b/ProcessLib/HeatConduction/HeatConductionProcess.cpp new file mode 100644 index 0000000000000000000000000000000000000000..1152cd0db1a24e562bc7a5244b09b93b607d1e2d --- /dev/null +++ b/ProcessLib/HeatConduction/HeatConductionProcess.cpp @@ -0,0 +1,81 @@ +/** + * \copyright + * Copyright (c) 2012-2016, OpenGeoSys Community (http://www.opengeosys.org) + * Distributed under a Modified BSD License. + * See accompanying file LICENSE.txt or + * http://www.opengeosys.org/project/license + * + */ + +#include "HeatConductionProcess.h" + +#include <cassert> + +#include "ProcessLib/Utils/CreateLocalAssemblers.h" + +namespace ProcessLib +{ +namespace HeatConduction +{ +HeatConductionProcess::HeatConductionProcess( + MeshLib::Mesh& mesh, + std::vector<std::unique_ptr<ParameterBase>> const& parameters, + std::vector<std::reference_wrapper<ProcessVariable>>&& process_variables, + HeatConductionProcessData&& process_data, + SecondaryVariableCollection&& secondary_variables, + NumLib::NamedFunctionCaller&& named_function_caller) + : Process(mesh, parameters, std::move(process_variables), + std::move(secondary_variables), std::move(named_function_caller)), + _process_data(std::move(process_data)) +{ +} + +void HeatConductionProcess::initializeConcreteProcess( + NumLib::LocalToGlobalIndexMap const& dof_table, + MeshLib::Mesh const& mesh, + unsigned const integration_order) +{ + ProcessLib::createLocalAssemblers<LocalAssemblerData>( + mesh.getDimension(), mesh.getElements(), dof_table, integration_order, + _local_assemblers, _process_data); + + _secondary_variables.addSecondaryVariable( + "heat_flux_x", 1, + makeExtrapolator( + getExtrapolator(), _local_assemblers, + &HeatConductionLocalAssemblerInterface::getIntPtHeatFluxX)); + + if (mesh.getDimension() > 1) + { + _secondary_variables.addSecondaryVariable( + "heat_flux_y", 1, + makeExtrapolator( + getExtrapolator(), _local_assemblers, + &HeatConductionLocalAssemblerInterface::getIntPtHeatFluxY)); + } + if (mesh.getDimension() > 2) + { + _secondary_variables.addSecondaryVariable( + "heat_flux_z", 1, + makeExtrapolator( + getExtrapolator(), _local_assemblers, + &HeatConductionLocalAssemblerInterface::getIntPtHeatFluxZ)); + } +} + +void HeatConductionProcess::assembleConcreteProcess(const double t, + GlobalVector const& x, + GlobalMatrix& M, + GlobalMatrix& K, + GlobalVector& b) +{ + DBUG("Assemble HeatConductionProcess."); + + // Call global assembler for each local assembly item. + GlobalExecutor::executeMemberOnDereferenced( + &HeatConductionLocalAssemblerInterface::assemble, _local_assemblers, + *_local_to_global_index_map, t, x, M, K, b); +} + +} // namespace HeatConduction +} // namespace ProcessLib diff --git a/ProcessLib/HeatConduction/HeatConductionProcess.h b/ProcessLib/HeatConduction/HeatConductionProcess.h new file mode 100644 index 0000000000000000000000000000000000000000..280eb5023344a25f60c2cf8e2ba2f79400df710b --- /dev/null +++ b/ProcessLib/HeatConduction/HeatConductionProcess.h @@ -0,0 +1,60 @@ +/** + * \copyright + * Copyright (c) 2012-2016, OpenGeoSys Community (http://www.opengeosys.org) + * Distributed under a Modified BSD License. + * See accompanying file LICENSE.txt or + * http://www.opengeosys.org/project/license + * + */ + +#ifndef PROCESS_LIB_HEATCONDUCTIONPROCESS_H_ +#define PROCESS_LIB_HEATCONDUCTIONPROCESS_H_ + +#include "NumLib/Extrapolation/LocalLinearLeastSquaresExtrapolator.h" +#include "ProcessLib/Process.h" +#include "HeatConductionFEM.h" +#include "HeatConductionProcessData.h" + +namespace ProcessLib +{ +namespace HeatConduction +{ +class HeatConductionProcess final : public Process +{ + +public: + HeatConductionProcess( + MeshLib::Mesh& mesh, + std::vector<std::unique_ptr<ParameterBase>> const& parameters, + std::vector<std::reference_wrapper<ProcessVariable>>&& + process_variables, + HeatConductionProcessData&& process_data, + SecondaryVariableCollection&& secondary_variables, + NumLib::NamedFunctionCaller&& named_function_caller); + + //! \name ODESystem interface + //! @{ + + bool isLinear() const override { return true; } + //! @} + +private: + void initializeConcreteProcess( + NumLib::LocalToGlobalIndexMap const& dof_table, + MeshLib::Mesh const& mesh, + unsigned const integration_order) override; + + void assembleConcreteProcess(const double t, GlobalVector const& x, + GlobalMatrix& M, GlobalMatrix& K, + GlobalVector& b) override; + + HeatConductionProcessData _process_data; + + std::vector<std::unique_ptr<HeatConductionLocalAssemblerInterface>> + _local_assemblers; +}; + +} // namespace HeatConduction +} // namespace ProcessLib + +#endif // PROCESS_LIB_HEATCONDUCTIONPROCESS_H_ diff --git a/ProcessLib/HeatConduction/HeatConductionProcessData.h b/ProcessLib/HeatConduction/HeatConductionProcessData.h new file mode 100644 index 0000000000000000000000000000000000000000..b95012d1ab941ae502c6a6491bb7690a529ab4f2 --- /dev/null +++ b/ProcessLib/HeatConduction/HeatConductionProcessData.h @@ -0,0 +1,61 @@ +/** + * \copyright + * Copyright (c) 2012-2016, OpenGeoSys Community (http://www.opengeosys.org) + * Distributed under a Modified BSD License. + * See accompanying file LICENSE.txt or + * http://www.opengeosys.org/project/license + * + */ + +#ifndef PROCESSLIB_HEATCONDUCTION_HEATCONDUCTIONPROCESSDATA_H +#define PROCESSLIB_HEATCONDUCTION_HEATCONDUCTIONPROCESSDATA_H + +namespace MeshLib +{ +class Element; +} + +namespace ProcessLib +{ +template <typename T> +struct Parameter; + +namespace HeatConduction +{ + +struct HeatConductionProcessData +{ + HeatConductionProcessData(Parameter<double> const& thermal_conductivity_, + Parameter<double> const& heat_capacity_, + Parameter<double> const& density_) + : thermal_conductivity(thermal_conductivity_), + heat_capacity(heat_capacity_), + density(density_) + { + } + + HeatConductionProcessData(HeatConductionProcessData&& other) + : thermal_conductivity(other.thermal_conductivity), + heat_capacity(other.heat_capacity), + density(other.density) + { + } + + //! Copies are forbidden. + HeatConductionProcessData(HeatConductionProcessData const&) = delete; + + //! Assignments are not needed. + void operator=(HeatConductionProcessData const&) = delete; + + //! Assignments are not needed. + void operator=(HeatConductionProcessData&&) = delete; + + Parameter<double> const& thermal_conductivity; + Parameter<double> const& heat_capacity; + Parameter<double> const& density; +}; + +} // namespace HeatConduction +} // namespace ProcessLib + +#endif // PROCESSLIB_HEATCONDUCTION_HEATCONDUCTIONPROCESSDATA_H diff --git a/Tests/Data b/Tests/Data index 7c791bb3c554551fa381f17c8f7a04ff11e99039..117219b7e7a6ae337257bb50c49730cdc818b26b 160000 --- a/Tests/Data +++ b/Tests/Data @@ -1 +1 @@ -Subproject commit 7c791bb3c554551fa381f17c8f7a04ff11e99039 +Subproject commit 117219b7e7a6ae337257bb50c49730cdc818b26b