Commit 734a9868 authored by Dmitry Yu. Naumov's avatar Dmitry Yu. Naumov

Merge branch 'ChemicalBC' into 'master'

[CL] Solution dependent Dirichlet boundary condition

See merge request !3085
parents f2d152d1 b514b165
Pipeline #1448 passed with stages
in 118 minutes and 8 seconds
......@@ -62,6 +62,13 @@ public:
// A hook added for solution dependent dirichlet
}
virtual void postTimestep(const double /*t*/,
std::vector<GlobalVector*> const& /*x*/,
int const /*process_id*/)
{
// A hook added for solution dependent dirichlet
}
virtual ~BoundaryCondition() = default;
};
......
......@@ -63,6 +63,15 @@ public:
}
}
void postTimestep(const double t, std::vector<GlobalVector*> const& x,
int const process_id)
{
for (auto const& bc_ptr : _boundary_conditions)
{
bc_ptr->postTimestep(t, x, process_id);
}
}
private:
mutable std::vector<NumLib::IndexValueVector<GlobalIndexType>> _dirichlet_bcs;
std::vector<std::unique_ptr<BoundaryCondition>> _boundary_conditions;
......
......@@ -21,6 +21,7 @@
#include "PhaseFieldIrreversibleDamageOracleBoundaryCondition.h"
#include "PrimaryVariableConstraintDirichletBoundaryCondition.h"
#include "RobinBoundaryCondition.h"
#include "SolutionDependentDirichletBoundaryCondition.h"
#include "VariableDependentNeumannBoundaryCondition.h"
#include "BaseLib/TimeInterval.h"
......@@ -114,6 +115,12 @@ std::unique_ptr<BoundaryCondition> createBoundaryCondition(
config.config, config.boundary_mesh, dof_table, variable_id,
*config.component_id, parameters);
}
if (type == "SolutionDependentDirichlet")
{
return ProcessLib::createSolutionDependentDirichletBoundaryCondition(
config.config, config.boundary_mesh, dof_table, variable_id,
*config.component_id, parameters);
}
if (type == "HCNonAdvectiveFreeComponentFlowBoundary")
{
return createHCNonAdvectiveFreeComponentFlowBoundaryCondition(
......
/**
* \file
* \copyright
* Copyright (c) 2012-2020, 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 "SolutionDependentDirichletBoundaryCondition.h"
#include <memory>
#include <vector>
#include "BaseLib/ConfigTree.h"
#include "BaseLib/Logging.h"
#include "DirichletBoundaryConditionAuxiliaryFunctions.h"
#include "MeshLib/Mesh.h"
#include "NumLib/DOF/LocalToGlobalIndexMap.h"
#include "ParameterLib/Utils.h"
namespace ProcessLib
{
SolutionDependentDirichletBoundaryCondition::
SolutionDependentDirichletBoundaryCondition(
ParameterLib::Parameter<double> const& parameter,
MeshLib::Mesh const& bc_mesh,
NumLib::LocalToGlobalIndexMap const& dof_table_bulk,
int const variable_id, int const component_id)
: _bc_mesh(bc_mesh), _variable_id(variable_id), _component_id(component_id)
{
checkParametersOfDirichletBoundaryCondition(_bc_mesh, dof_table_bulk,
_variable_id, _component_id);
std::vector<MeshLib::Node*> const& bc_nodes = bc_mesh.getNodes();
MeshLib::MeshSubset bc_mesh_subset(_bc_mesh, bc_nodes);
// Create local DOF table from the BC mesh subset for the given variable
// and component id.
_dof_table_boundary.reset(dof_table_bulk.deriveBoundaryConstrainedMap(
variable_id, {component_id}, std::move(bc_mesh_subset)));
std::string const property_name = "solution_dependent_bc";
if (bc_mesh.getProperties().existsPropertyVector<double>(property_name))
{
OGS_FATAL(
"Found mesh property '{:s}' which is the built-in property of the "
"class SolutionDependentDirichletBoundaryCondition.",
property_name);
}
auto& solution_dependent_bc = *MeshLib::getOrCreateMeshProperty<double>(
const_cast<MeshLib::Mesh&>(bc_mesh), property_name,
MeshLib::MeshItemType::Node, 1);
solution_dependent_bc.resize(bc_mesh.getNumberOfNodes());
ParameterLib::SpatialPosition pos;
auto const& nodes = bc_mesh.getNodes();
for (std::size_t i = 0; i < solution_dependent_bc.size(); ++i)
{
auto const id = nodes[i]->getID();
pos.setNodeID(id);
solution_dependent_bc[i] = parameter(0, pos)[0];
}
_parameter = std::make_unique<ParameterLib::MeshNodeParameter<double>>(
property_name, bc_mesh, solution_dependent_bc);
}
void SolutionDependentDirichletBoundaryCondition::getEssentialBCValues(
double const t, GlobalVector const& x,
NumLib::IndexValueVector<GlobalIndexType>& bc_values) const
{
getEssentialBCValuesLocal(*_parameter, _bc_mesh, _bc_mesh.getNodes(),
*_dof_table_boundary, _variable_id, _component_id,
t, x, bc_values);
}
void SolutionDependentDirichletBoundaryCondition::postTimestep(
double const /*t*/, std::vector<GlobalVector*> const& x,
int const process_id)
{
auto& solution_dependent_bc =
*const_cast<MeshLib::Properties&>(_bc_mesh.getProperties())
.getPropertyVector<double>("solution_dependent_bc");
auto const& nodes = _bc_mesh.getNodes();
for (std::size_t i = 0; i < solution_dependent_bc.size(); ++i)
{
auto const id = nodes[i]->getID();
auto const global_index = _dof_table_boundary->getGlobalIndex(
{_bc_mesh.getID(), MeshLib::MeshItemType::Node, id}, _variable_id,
_component_id);
assert(global_index >= 0);
solution_dependent_bc[i] = x[process_id]->get(global_index);
}
}
std::unique_ptr<SolutionDependentDirichletBoundaryCondition>
createSolutionDependentDirichletBoundaryCondition(
BaseLib::ConfigTree const& config, MeshLib::Mesh const& bc_mesh,
NumLib::LocalToGlobalIndexMap const& dof_table_bulk, int const variable_id,
int const component_id,
std::vector<std::unique_ptr<ParameterLib::ParameterBase>> const& parameters)
{
DBUG(
"Constructing SolutionDependentDirichletBoundaryCondition from "
"config.");
//! \ogs_file_param{prj__process_variables__process_variable__boundary_conditions__boundary_condition__type}
config.checkConfigParameter("type", "SolutionDependentDirichlet");
auto& initial_value_parameter = ParameterLib::findParameter<double>(
//! \ogs_file_param{prj__process_variables__process_variable__boundary_conditions__boundary_condition__SolutionDependentDirichlet__initial_value_parameter}
config.getConfigParameter<std::string>("initial_value_parameter"),
parameters, 1, &bc_mesh);
// In case of partitioned mesh the boundary could be empty, i.e. there is no
// boundary condition.
#ifdef USE_PETSC
// This can be extracted to createBoundaryCondition() but then the config
// parameters are not read and will cause an error.
// TODO (naumov): Add a function to ConfigTree for skipping the tags of the
// subtree and move the code up in createBoundaryCondition().
if (bc_mesh.getDimension() == 0 && bc_mesh.getNumberOfNodes() == 0 &&
bc_mesh.getNumberOfElements() == 0)
{
return nullptr;
}
#endif // USE_PETSC
return std::make_unique<SolutionDependentDirichletBoundaryCondition>(
initial_value_parameter, bc_mesh, dof_table_bulk, variable_id,
component_id);
}
} // namespace ProcessLib
/**
* \file
* \copyright
* Copyright (c) 2012-2020, OpenGeoSys Community (http://www.opengeosys.org)
* Distributed under a Modified BSD License.
* See accompanying file LICENSE.txt or
* http://www.opengeosys.org/project/license
*
*/
#pragma once
#include "BoundaryCondition.h"
#include "ParameterLib/MeshNodeParameter.h"
namespace ProcessLib
{
/// The SolutionDependentDirichletBoundaryCondition belongs to the
/// Dirichlet-type boundary condition.
///
/// This class is a special category of Dirichlet boundary condition,
/// applied in the situation where the value assigned for the boundary condition
/// is dependent on the process solution of the last time step. This particular
/// boundary condition is widely used in the reactive transport problems and has
/// the potential to be used in other processes.
class SolutionDependentDirichletBoundaryCondition final
: public BoundaryCondition
{
public:
SolutionDependentDirichletBoundaryCondition(
ParameterLib::Parameter<double> const& parameter,
MeshLib::Mesh const& bc_mesh,
NumLib::LocalToGlobalIndexMap const& dof_table_bulk,
int const variable_id, int const component_id);
void getEssentialBCValues(
double const t, GlobalVector const& x,
NumLib::IndexValueVector<GlobalIndexType>& bc_values) const override;
/// Renchao: The original idea to place the update of the boundary condition
/// value at the preTimestep stage. The update could be achieved within the
/// class function "Process::preTimestep". Whereas, I find it not doable to
/// implement in this way. The class function "Process::preTimestep" is
/// called in a row by the functions "preTimestepForAllProcesses" and
/// "TimeLoop::outputSolutions". These two functions are called when
/// initializing and subsequently looping over the "TimeLoop". Actually, it
/// is not intended to make the boundary condition value updated in the
/// first loop. Instead, the update is intended to start from the second
/// loop. For these two reasons, I think it more proper to do the
/// implementation at the postTimestep stage.
void postTimestep(double const /*t*/,
std::vector<GlobalVector*> const& x,
int const process_id) override;
private:
MeshLib::Mesh const& _bc_mesh;
int const _variable_id;
int const _component_id;
std::unique_ptr<NumLib::LocalToGlobalIndexMap const> _dof_table_boundary;
std::unique_ptr<ParameterLib::MeshNodeParameter<double>> _parameter;
};
std::unique_ptr<SolutionDependentDirichletBoundaryCondition>
createSolutionDependentDirichletBoundaryCondition(
BaseLib::ConfigTree const& config, MeshLib::Mesh const& bc_mesh,
NumLib::LocalToGlobalIndexMap const& dof_table_bulk, int const variable_id,
int const component_id,
std::vector<std::unique_ptr<ParameterLib::ParameterBase>> const&
parameters);
} // namespace ProcessLib
......@@ -384,6 +384,8 @@ void Process::postTimestep(std::vector<GlobalVector*> const& x, const double t,
for (auto* const solution : x)
MathLib::LinAlg::setLocalAccessibleVector(*solution);
postTimestepConcreteProcess(x, t, delta_t, process_id);
_boundary_conditions[process_id].postTimestep(t, x, process_id);
}
void Process::postNonLinearSolver(GlobalVector const& x, const double t,
......
......@@ -438,8 +438,8 @@
<boundary_conditions>
<boundary_condition>
<mesh>1d_isofrac_upstream</mesh>
<type>Dirichlet</type>
<parameter>c_Synthetica</parameter>
<type>SolutionDependentDirichlet</type>
<initial_value_parameter>c_Synthetica</initial_value_parameter>
</boundary_condition>
</boundary_conditions>
</process_variable>
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment