From 9b47442c8481b5e66a31dbd094825e58823e257c Mon Sep 17 00:00:00 2001
From: Christoph Lehmann <christoph.lehmann@ufz.de>
Date: Fri, 19 Aug 2016 16:51:12 +0200
Subject: [PATCH] [PL] made Process Jacobian-assembly-ready
---
ProcessLib/LocalAssemblerInterface.cpp | 40 ++++++-----------------
ProcessLib/LocalAssemblerInterface.h | 38 +++++++++++-----------
ProcessLib/Process.cpp | 45 ++++++--------------------
ProcessLib/Process.h | 22 +++++++++----
4 files changed, 53 insertions(+), 92 deletions(-)
diff --git a/ProcessLib/LocalAssemblerInterface.cpp b/ProcessLib/LocalAssemblerInterface.cpp
index e4567e3b962..5e572b76b05 100644
--- a/ProcessLib/LocalAssemblerInterface.cpp
+++ b/ProcessLib/LocalAssemblerInterface.cpp
@@ -8,43 +8,21 @@
*/
#include "LocalAssemblerInterface.h"
+#include <cassert>
#include "NumLib/DOF/DOFTableUtil.h"
namespace ProcessLib
{
-void LocalAssemblerInterface::assemble(
- const std::size_t mesh_item_id,
- const NumLib::LocalToGlobalIndexMap& dof_table, const double t,
- const GlobalVector& x, GlobalMatrix& M, GlobalMatrix& K, GlobalVector& b)
-{
- auto const indices = NumLib::getIndices(mesh_item_id, dof_table);
- auto const local_x = x.get(indices);
- auto const r_c_indices =
- NumLib::LocalToGlobalIndexMap::RowColumnIndices(indices, indices);
-
- assembleConcrete(t, local_x, r_c_indices, M, K, b);
-}
-
-void LocalAssemblerInterface::assembleJacobian(
- const std::size_t mesh_item_id,
- const NumLib::LocalToGlobalIndexMap& dof_table, const double t,
- const GlobalVector& x, GlobalMatrix& Jac)
-{
- auto const indices = NumLib::getIndices(mesh_item_id, dof_table);
- auto const local_x = x.get(indices);
- auto const r_c_indices =
- NumLib::LocalToGlobalIndexMap::RowColumnIndices(indices, indices);
-
- assembleJacobianConcrete(t, local_x, r_c_indices, Jac);
-}
-
-void LocalAssemblerInterface::assembleJacobianConcrete(
- double const /*t*/, std::vector<double> const& /*local_x*/,
- NumLib::LocalToGlobalIndexMap::RowColumnIndices const& /*indices*/,
- GlobalMatrix& /*Jac*/)
+void LocalAssemblerInterface::assembleWithJacobian(
+ double const /*t*/, std::vector<double> const& /*local_x*/,
+ std::vector<double> const& /*local_xdot*/, const double /*dxdot_dx*/,
+ const double /*dx_dx*/, std::vector<double>& /*local_M_data*/,
+ std::vector<double>& /*local_K_data*/,
+ std::vector<double>& /*local_b_data*/,
+ std::vector<double>& /*local_Jac_data*/)
{
OGS_FATAL(
- "The assembleJacobian() function is not implemented in the local "
+ "The assembleWithJacobian() function is not implemented in the local "
"assembler.");
}
diff --git a/ProcessLib/LocalAssemblerInterface.h b/ProcessLib/LocalAssemblerInterface.h
index 574b6444724..355bb5d8f1c 100644
--- a/ProcessLib/LocalAssemblerInterface.h
+++ b/ProcessLib/LocalAssemblerInterface.h
@@ -10,9 +10,13 @@
#ifndef PROCESSLIB_LOCALASSEMBLERINTERFACE_H
#define PROCESSLIB_LOCALASSEMBLERINTERFACE_H
-#include "NumLib/DOF/LocalToGlobalIndexMap.h"
#include "NumLib/NumericsConfig.h"
+namespace NumLib
+{
+class LocalToGlobalIndexMap;
+} // NumLib
+
namespace ProcessLib
{
@@ -26,15 +30,19 @@ class LocalAssemblerInterface
public:
virtual ~LocalAssemblerInterface() = default;
- void assemble(std::size_t const mesh_item_id,
- NumLib::LocalToGlobalIndexMap const& dof_table,
- double const t, GlobalVector const& x, GlobalMatrix& M,
- GlobalMatrix& K, GlobalVector& b);
+ virtual void assemble(
+ double const t, std::vector<double> const& local_x,
+ std::vector<double>& local_M_data, std::vector<double>& local_K_data,
+ std::vector<double>& local_b_data) = 0;
- void assembleJacobian(std::size_t const mesh_item_id,
- NumLib::LocalToGlobalIndexMap const& dof_table,
- double const t, GlobalVector const& x,
- GlobalMatrix& Jac);
+ virtual void assembleWithJacobian(double const t,
+ std::vector<double> const& local_x,
+ std::vector<double> const& local_xdot,
+ const double dxdot_dx, const double dx_dx,
+ std::vector<double>& local_M_data,
+ std::vector<double>& local_K_data,
+ std::vector<double>& local_b_data,
+ std::vector<double>& local_Jac_data);
virtual void preTimestep(std::size_t const mesh_item_id,
NumLib::LocalToGlobalIndexMap const& dof_table,
@@ -45,17 +53,7 @@ public:
NumLib::LocalToGlobalIndexMap const& dof_table,
GlobalVector const& x);
-protected:
- virtual void assembleConcrete(
- double const t, std::vector<double> const& local_x,
- NumLib::LocalToGlobalIndexMap::RowColumnIndices const& indices,
- GlobalMatrix& M, GlobalMatrix& K, GlobalVector& b) = 0;
-
- virtual void assembleJacobianConcrete(
- double const t, std::vector<double> const& local_x,
- NumLib::LocalToGlobalIndexMap::RowColumnIndices const& indices,
- GlobalMatrix& Jac);
-
+private:
virtual void preTimestepConcrete(std::vector<double> const& /*local_x*/,
double const /*t*/, double const /*dt*/)
{
diff --git a/ProcessLib/Process.cpp b/ProcessLib/Process.cpp
index 7c21ef4da2e..407b6e12fb2 100644
--- a/ProcessLib/Process.cpp
+++ b/ProcessLib/Process.cpp
@@ -20,6 +20,7 @@ namespace ProcessLib
{
Process::Process(
MeshLib::Mesh& mesh,
+ std::unique_ptr<ProcessLib::AbstractJacobianAssembler>&& jacobian_assembler,
std::vector<std::unique_ptr<ParameterBase>> const& parameters,
std::vector<std::reference_wrapper<ProcessVariable>>&& process_variables,
SecondaryVariableCollection&& secondary_variables,
@@ -27,6 +28,7 @@ Process::Process(
: _mesh(mesh),
_secondary_variables(std::move(secondary_variables)),
_named_function_caller(std::move(named_function_caller)),
+ _global_assembler(std::move(jacobian_assembler)),
_process_variables(std::move(process_variables)),
_boundary_conditions(parameters)
{
@@ -115,43 +117,16 @@ void Process::assemble(const double t, GlobalVector const& x, GlobalMatrix& M,
_boundary_conditions.applyNaturalBC(t, x, K, b);
}
-void Process::assembleJacobian(const double t, GlobalVector const& x,
- GlobalVector const& xdot, const double dxdot_dx,
- GlobalMatrix const& M, const double dx_dx,
- GlobalMatrix const& K, GlobalMatrix& Jac)
+void Process::assembleWithJacobian(const double t, GlobalVector const& x,
+ GlobalVector const& xdot,
+ const double dxdot_dx, const double dx_dx,
+ GlobalMatrix& M, GlobalMatrix& K,
+ GlobalVector& b, GlobalMatrix& Jac)
{
- assembleJacobianConcreteProcess(t, x, xdot, dxdot_dx, M, dx_dx, K, Jac);
-
- // TODO In this method one could check if the user wants to use an
- // analytical or a numerical Jacobian. Then the right
- // assembleJacobianConcreteProcess() method will be chosen.
- // Additionally in the default implementation of said method one
- // could provide a fallback to a numerical Jacobian. However, that
- // would be in a sense implicit behaviour and it might be better to
- // just abort, as is currently the case.
- // In order to implement the Jacobian assembly entirely, in addition
- // to the operator() in VectorMatrixAssembler there has to be a method
- // that dispatches the Jacobian assembly.
- // Similarly, then the NeumannBC specialization of VectorMatrixAssembler
- // probably can be merged into the main class s.t. one has only one
- // type of VectorMatrixAssembler (for each equation type) with the
- // three methods assemble(), assembleJacobian() and assembleNeumannBC().
- // That list can be extended, e.g. by methods for the assembly of
- // source terms.
- // UPDATE: Probably it is better to keep a separate NeumannBC version of the
- // VectoMatrixAssembler since that will work for all kinds of processes.
-}
+ assembleWithJacobianConcreteProcess(t, x, xdot, dxdot_dx, dx_dx, M, K, b, Jac);
-void Process::assembleJacobianConcreteProcess(
- const double /*t*/, GlobalVector const& /*x*/, GlobalVector const& /*xdot*/,
- const double /*dxdot_dx*/, GlobalMatrix const& /*M*/,
- const double /*dx_dx*/, GlobalMatrix const& /*K*/, GlobalMatrix& /*Jac*/)
-{
- OGS_FATAL(
- "The concrete implementation of this Process did not override the"
- " assembleJacobianConcreteProcess() method."
- " Hence, no analytical Jacobian is provided for this process"
- " and the Newton-Raphson method cannot be used to solve it.");
+ // TODO apply BCs to Jacobian.
+ _boundary_conditions.applyNaturalBC(t, x, K, b);
}
void Process::constructDofTable()
diff --git a/ProcessLib/Process.h b/ProcessLib/Process.h
index ef9c60dc754..9ba6f93afa6 100644
--- a/ProcessLib/Process.h
+++ b/ProcessLib/Process.h
@@ -40,6 +40,8 @@ public:
Process(
MeshLib::Mesh& mesh,
+ std::unique_ptr<ProcessLib::AbstractJacobianAssembler>&&
+ jacobian_assembler,
std::vector<std::unique_ptr<ParameterBase>> const& parameters,
std::vector<std::reference_wrapper<ProcessVariable>>&&
process_variables,
@@ -68,14 +70,14 @@ public:
void assemble(const double t, GlobalVector const& x, GlobalMatrix& M,
GlobalMatrix& K, GlobalVector& b) override final;
- void assembleJacobian(const double t, GlobalVector const& x,
- GlobalVector const& xdot, const double dxdot_dx,
- GlobalMatrix const& M, const double dx_dx,
- GlobalMatrix const& K,
- GlobalMatrix& Jac) override final;
+ void assembleWithJacobian(const double t, GlobalVector const& x,
+ GlobalVector const& xdot, const double dxdot_dx,
+ const double dx_dx, GlobalMatrix& M,
+ GlobalMatrix& K, GlobalVector& b,
+ GlobalMatrix& Jac) override final;
+
std::vector<NumLib::IndexValueVector<GlobalIndexType>> const*
getKnownSolutions(
-
double const t) const override final
{
return _boundary_conditions.getKnownSolutions(t);
@@ -121,6 +123,12 @@ private:
GlobalMatrix& M, GlobalMatrix& K,
GlobalVector& b) = 0;
+ virtual void assembleWithJacobianConcreteProcess(
+ const double t, GlobalVector const& x,
+ GlobalVector const& xdot, const double dxdot_dx,
+ const double dx_dx, GlobalMatrix& M, GlobalMatrix& K,
+ GlobalVector& b, GlobalMatrix& Jac) = 0;
+
virtual void preTimestepConcreteProcess(GlobalVector const& /*x*/,
const double /*t*/,
const double /*delta_t*/)
@@ -169,6 +177,8 @@ protected:
_cached_secondary_variables;
SecondaryVariableContext _secondary_variable_context;
+ VectorMatrixAssembler _global_assembler;
+
private:
unsigned const _integration_order = 2;
GlobalSparsityPattern _sparsity_pattern;
--
GitLab