Detailed Description

Assembles the Jacobian matrix using central differences.

Definition at line 24 of file CentralDifferencesJacobianAssembler.h.

#include <CentralDifferencesJacobianAssembler.h>

Inheritance diagram for ProcessLib::CentralDifferencesJacobianAssembler:

Collaboration diagram for ProcessLib::CentralDifferencesJacobianAssembler:

Public Member Functions
	CentralDifferencesJacobianAssembler (std::vector< double > &&absolute_epsilons)

void	assembleWithJacobian (LocalAssemblerInterface &local_assembler, double const t, double const dt, 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) override

Public Member Functions inherited from ProcessLib::AbstractJacobianAssembler
virtual void	assembleWithJacobianForStaggeredScheme (LocalAssemblerInterface &, double const, double const, Eigen::VectorXd const &, Eigen::VectorXd const &, const double, const double, int const, std::vector< double > &, std::vector< double > &, std::vector< double > &, std::vector< double > &)

virtual	~AbstractJacobianAssembler ()=default

Private Attributes
std::vector< double > const	_absolute_epsilons

std::vector< double >	_local_M_data

std::vector< double >	_local_K_data

std::vector< double >	_local_b_data

std::vector< double >	_local_x_perturbed_data

std::vector< double >	_local_xdot_perturbed_data

Constructor & Destructor Documentation

◆ CentralDifferencesJacobianAssembler()

ProcessLib::CentralDifferencesJacobianAssembler::CentralDifferencesJacobianAssembler ( std::vector< double > && absolute_epsilons )

explicit

Constructs a new instance.

Parameters

absolute_epsilons perturbations of the components of the local solution vector used for evaluating the finite differences.

Note: The size of absolute_epsilons defines the "number of components" of the local solution vector (This is not the number of elements of the vector!). Therefore the size of the local solution vector must be divisible by the size of absolute_epsilons. This is the only consistency check performed. It is not checked whether said "number of components" is sensible. E.g., one could pass one epsilon per node, which would be valid but would not make sense at all.

Definition at line 20 of file CentralDifferencesJacobianAssembler.cpp.

     : _absolute_epsilons(std::move(absolute_epsilons))
 {
     if (_absolute_epsilons.empty())
     {
         OGS_FATAL("No values for the absolute epsilons have been given.");
     }
 }

References _absolute_epsilons, and OGS_FATAL.

Member Function Documentation

◆ assembleWithJacobian()

void ProcessLib::CentralDifferencesJacobianAssembler::assembleWithJacobian	(	LocalAssemblerInterface &	local_assembler,
		double const	t,
		double const	dt,
		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
	)

overridevirtual

Assembles the Jacobian, the matrices \(M\) and \(K\), and the vector \(b\). For the assembly the assemble() method of the given local_assembler is called several times and the Jacobian is built from finite differences. The number of calls of the assemble() method is \(2N+1\) if \(N\) is the size of local_x.

Attention: It is assumed that the local vectors and matrices are ordered by component.

Implements ProcessLib::AbstractJacobianAssembler.

Definition at line 30 of file CentralDifferencesJacobianAssembler.cpp.

 {
     // TODO do not check in every call.
     if (local_x_data.size() % _absolute_epsilons.size() != 0)
     {
         OGS_FATAL(
             "The number of specified epsilons ({:d}) and the number of local "
             "d.o.f.s ({:d}) do not match, i.e., the latter is not divisible by "
             "the former.",
             _absolute_epsilons.size(), local_x_data.size());
     }
  
     auto const num_r_c =
         static_cast<Eigen::MatrixXd::Index>(local_x_data.size());
  
     auto const local_x =
         MathLib::toVector<Eigen::VectorXd>(local_x_data, num_r_c);
     auto const local_xdot =
         MathLib::toVector<Eigen::VectorXd>(local_xdot_data, num_r_c);
  
     auto local_Jac =
         MathLib::createZeroedMatrix(local_Jac_data, num_r_c, num_r_c);
     _local_x_perturbed_data = local_x_data;
     _local_xdot_perturbed_data = local_xdot_data;
  
     auto const num_dofs_per_component =
         local_x_data.size() / _absolute_epsilons.size();
  
     // Residual  res := M xdot + K x - b
     // Computing Jac := dres/dx
     //                = M dxdot/dx + dM/dx xdot + K dx/dx + dK/dx x - db/dx
     //                  (Note: dM/dx and dK/dx actually have the second and
     //                  third index transposed.)
     // The loop computes the dM/dx, dK/dx and db/dx terms, the rest is computed
     // afterwards.
     for (Eigen::MatrixXd::Index i = 0; i < num_r_c; ++i)
     {
         // assume that local_x_data is ordered by component.
         auto const component = i / num_dofs_per_component;
         auto const eps = _absolute_epsilons[component];
  
         _local_x_perturbed_data[i] += eps;
         _local_xdot_perturbed_data[i] = local_xdot_data[i] + eps / dt;
         local_assembler.assemble(t, dt, _local_x_perturbed_data,
                                  _local_xdot_perturbed_data, local_M_data,
                                  local_K_data, local_b_data);
  
         _local_x_perturbed_data[i] = local_x_data[i] - eps;
         _local_xdot_perturbed_data[i] = local_xdot_data[i] - eps / dt;
         local_assembler.assemble(t, dt, _local_x_perturbed_data,
                                  _local_xdot_perturbed_data, _local_M_data,
                                  _local_K_data, _local_b_data);
  
         _local_x_perturbed_data[i] = local_x_data[i];
         _local_xdot_perturbed_data[i] = local_xdot_data[i];
  
         if (!local_M_data.empty())
         {
             auto const local_M_p =
                 MathLib::toMatrix(local_M_data, num_r_c, num_r_c);
             auto const local_M_m =
                 MathLib::toMatrix(_local_M_data, num_r_c, num_r_c);
             local_Jac.col(i).noalias() +=
                 // dM/dxi * x_dot
                 (local_M_p - local_M_m) * local_xdot / (2.0 * eps);
             local_M_data.clear();
             _local_M_data.clear();
         }
         if (!local_K_data.empty())
         {
             auto const local_K_p =
                 MathLib::toMatrix(local_K_data, num_r_c, num_r_c);
             auto const local_K_m =
                 MathLib::toMatrix(_local_K_data, num_r_c, num_r_c);
             local_Jac.col(i).noalias() +=
                 // dK/dxi * x
                 (local_K_p - local_K_m) * local_x / (2.0 * eps);
             local_K_data.clear();
             _local_K_data.clear();
         }
         if (!local_b_data.empty())
         {
             auto const local_b_p =
                 MathLib::toVector<Eigen::VectorXd>(local_b_data, num_r_c);
             auto const local_b_m =
                 MathLib::toVector<Eigen::VectorXd>(_local_b_data, num_r_c);
             local_Jac.col(i).noalias() -=
                 // db/dxi
                 (local_b_p - local_b_m) / (2.0 * eps);
             local_b_data.clear();
             _local_b_data.clear();
         }
     }
  
     // Assemble with unperturbed local x.
     local_assembler.assemble(t, dt, local_x_data, local_xdot_data, local_M_data,
                              local_K_data, local_b_data);
  
     // Compute remaining terms of the Jacobian.
     if (dxdot_dx != 0.0 && !local_M_data.empty())
     {
         auto local_M = MathLib::toMatrix(local_M_data, num_r_c, num_r_c);
         local_Jac.noalias() += local_M * dxdot_dx;
     }
     if (dx_dx != 0.0 && !local_K_data.empty())
     {
         auto local_K = MathLib::toMatrix(local_K_data, num_r_c, num_r_c);
         local_Jac.noalias() += local_K * dx_dx;
     }
 }