Calculate directional derivative of the right hand side and compare it with tangent matrix derivative. More...

#include "src/interfaces/OperatorsTester.hpp"

Inheritance diagram for MoFEM::OperatorsTester:

Collaboration diagram for MoFEM::OperatorsTester:

Public Types
using	RandomFieldData = std::pair< std::string, std::array< double, 2 > >

Public Member Functions
MoFEMErrorCode	query_interface (boost::typeindex::type_index type_index, UnknownInterface **iface) const

	OperatorsTester (const MoFEM::Core &core)

virtual	~OperatorsTester ()=default

SmartPetscObj< Vec >	setRandomFields (SmartPetscObj< DM > dm, std::vector< RandomFieldData > random_fields, boost::shared_ptr< Range > ents=nullptr, RowColData r=RowColData::COL)
	Generate random fields.

SmartPetscObj< Vec >	assembleVec (SmartPetscObj< DM > dm, std::string fe_name, boost::shared_ptr< FEMethod > pipeline, SmartPetscObj< Vec > x, SmartPetscObj< Vec > delta_x, SmartPetscObj< Vec > delta2_x, double time, double delta_t, CacheTupleWeakPtr cache_ptr)
	Assemble the right hand side vector.

SmartPetscObj< Mat >	assembleMat (SmartPetscObj< DM > dm, std::string fe_name, boost::shared_ptr< FEMethod > pipeline, SmartPetscObj< Vec > x, SmartPetscObj< Vec > delta_x, SmartPetscObj< Vec > delta2_x, double time, double delta_t, CacheTupleWeakPtr cache_ptr)
	Assemble the left hand side vector.

SmartPetscObj< Vec >	directionalCentralFiniteDifference (SmartPetscObj< DM > dm, std::string fe_name, boost::shared_ptr< FEMethod > pipeline, SmartPetscObj< Vec > x, SmartPetscObj< Vec > delta_x, SmartPetscObj< Vec > delta2_x, SmartPetscObj< Vec > diff_x, double time, double delta_t, double eps, CacheTupleWeakPtr cache_ptr=CacheTupleSharedPtr())
	Calculate directional directive using finite difference.

SmartPetscObj< Vec >	checkCentralFiniteDifference (SmartPetscObj< DM > dm, std::string fe_name, boost::shared_ptr< FEMethod > pipeline_rhs, boost::shared_ptr< FEMethod > pipeline_lhs, SmartPetscObj< Vec > x, SmartPetscObj< Vec > delta_x, SmartPetscObj< Vec > delta2_x, SmartPetscObj< Vec > diff_x, double time, double delta_t, double eps, CacheTupleWeakPtr cache_ptr=CacheTupleSharedPtr())
	Check consistency between directional derivative with matrix.

Public Member Functions inherited from MoFEM::UnknownInterface
template<class IFACE >
MoFEMErrorCode	registerInterface (bool error_if_registration_failed=true)
	Register interface.

template<class IFACE >
MoFEMErrorCode	getInterface (IFACE *&iface) const
	Get interface reference to pointer of interface.

template<class IFACE >
MoFEMErrorCode	getInterface (IFACE **const iface) const
	Get interface pointer to pointer of interface.

template<class IFACE , typename boost::enable_if< boost::is_pointer< IFACE >, int >::type = 0>
IFACE	getInterface () const
	Get interface pointer to pointer of interface.

template<class IFACE , typename boost::enable_if< boost::is_reference< IFACE >, int >::type = 0>
IFACE	getInterface () const
	Get reference to interface.

template<class IFACE >
IFACE *	getInterface () const
	Function returning pointer to interface.

virtual	~UnknownInterface ()=default

Private Member Functions
std::pair< SmartPetscObj< Vec >, SmartPetscObj< Vec > >	setPipelineX (boost::shared_ptr< FEMethod > pipeline, SmartPetscObj< Vec > x, SmartPetscObj< Vec > delta_x, SmartPetscObj< Vec > delta2_x, double delta_t)
	Set vectors, x, x_t, and x_tt to finite element instance.

Private Attributes
MoFEM::Core &	cOre

Additional Inherited Members
Static Public Member Functions inherited from MoFEM::UnknownInterface
static MoFEMErrorCode	getLibVersion (Version &version)
	Get library version.

static MoFEMErrorCode	getFileVersion (moab::Interface &moab, Version &version)
	Get database major version.

static MoFEMErrorCode	setFileVersion (moab::Interface &moab, Version version=Version(MoFEM_VERSION_MAJOR, MoFEM_VERSION_MINOR, MoFEM_VERSION_BUILD))
	Get database major version.

static MoFEMErrorCode	getInterfaceVersion (Version &version)
	Get database major version.

Detailed Description

Calculate directional derivative of the right hand side and compare it with tangent matrix derivative.

Examples: mofem/atom_tests/operators_tests.cpp, mofem/atom_tests/tensor_divergence_operator.cpp, mofem/tutorials/adv-0_plasticity/plastic.cpp, mofem/tutorials/adv-3_level_set/level_set.cpp, and plastic.cpp.

Definition at line 21 of file OperatorsTester.hpp.

Member Typedef Documentation

◆ RandomFieldData

using MoFEM::OperatorsTester::RandomFieldData = std::pair<std::string, std::array<double, 2> >

Definition at line 29 of file OperatorsTester.hpp.

Constructor & Destructor Documentation

◆ OperatorsTester()

MoFEM::OperatorsTester::OperatorsTester ( const MoFEM::Core & core )

Definition at line 13 of file OperatorsTester.cpp.

14 : cOre(const_cast<Core &>(core)) {}

MoFEM::Core

CoreTmp< 0 > Core

Definition Core.hpp:1148

MoFEM::OperatorsTester::cOre

MoFEM::Core & cOre

Definition OperatorsTester.hpp:146

◆ ~OperatorsTester()

virtual MoFEM::OperatorsTester::~OperatorsTester ( )

virtualdefault

Member Function Documentation

◆ assembleMat()

SmartPetscObj< Mat > MoFEM::OperatorsTester::assembleMat	(	SmartPetscObj< DM >	dm,
		std::string	fe_name,
		boost::shared_ptr< FEMethod >	pipeline,
		SmartPetscObj< Vec >	x,
		SmartPetscObj< Vec >	delta_x,
		SmartPetscObj< Vec >	delta2_x,
		double	time,
		double	delta_t,
		CacheTupleWeakPtr	cache_ptr
	)

Assemble the left hand side vector.

Parameters

dm
fe_name	// fe name
pipeline	// pipeline, i.e. fe instance
x	// problem (dm) vector
delta_x	// vector for x rate, can be null, i.e, SmartPetscObj<Vec>()
delta2_x	// vector for x second rate, i.e. acceleration
time	// time
delta_t	// time increment
cache_ptr	// finite element data cache, can be null

Returns: SmartPetscObj<Vec>

Definition at line 148 of file OperatorsTester.cpp.

                                                                          {
  auto m = createDMMatrix(dm);
 
  pipeline->data_ctx =
      PetscData::CtxSetA | PetscData::CtxSetB | PetscData::CtxSetTime;
  pipeline->A = m;
  pipeline->B = m;
 
  pipeline->ts_t = time;
 
  pipeline->ts_dt = delta_t;
  pipeline->ts_a = 1. / delta_t;
  pipeline->ts_aa = 1. / pow(delta_t, 2);
 
  pipeline->ksp_ctx = KspMethod::KSPContext::CTX_OPERATORS;
  pipeline->snes_ctx = SnesMethod::SNESContext::CTX_SNESSETJACOBIAN;
  pipeline->ts_ctx = TSMethod::TSContext::CTX_TSSETIJACOBIAN;
 
  auto [v, a] = setPipelineX(pipeline, x, delta_x, delta2_x, delta_t);
  std::ignore = a;
  std::ignore = v;
 
  CHK_THROW_MESSAGE(
      DMoFEMMeshToLocalVector(dm, x, INSERT_VALUES, SCATTER_REVERSE),
      "loc to vec");
  CHK_THROW_MESSAGE(DMoFEMLoopFiniteElements(dm, fe_name, pipeline, cache_ptr),
                    "Run assemble pipeline");
  CHK_THROW_MESSAGE(MatAssemblyBegin(m, MAT_FINAL_ASSEMBLY), "Asmb mat");
  CHK_THROW_MESSAGE(MatAssemblyEnd(m, MAT_FINAL_ASSEMBLY), "Asmb mat");
 
  auto get_norm = [](auto m) {
    double fnorm;
    MatNorm(m, NORM_1, &fnorm);
    return fnorm;
  };
  MOFEM_LOG_CHANNEL("WORLD");
  MOFEM_TAG_AND_LOG_C("WORLD", Sev::noisy, "OpTester", "Matrix norm %3.4e",
                      get_norm(m));
 
  return m;
}

◆ assembleVec()

SmartPetscObj< Vec > MoFEM::OperatorsTester::assembleVec	(	SmartPetscObj< DM >	dm,
		std::string	fe_name,
		boost::shared_ptr< FEMethod >	pipeline,
		SmartPetscObj< Vec >	x,
		SmartPetscObj< Vec >	delta_x,
		SmartPetscObj< Vec >	delta2_x,
		double	time,
		double	delta_t,
		CacheTupleWeakPtr	cache_ptr
	)

Assemble the right hand side vector.

Parameters

dm
fe_name	// fe name
pipeline	// pipeline, i.e. fe instance
x	// problem (dm) vector
delta_x	// vector for x rate, can be null, i.e, SmartPetscObj<Vec>()
delta2_x	// vector for x second rate, i.e. acceleration
time	// time
delta_t	// time increment
cache_ptr	// finite element data cache, can be null

Returns: SmartPetscObj<Vec>

Definition at line 98 of file OperatorsTester.cpp.

                                                                          {
 
  pipeline->data_ctx = PetscData::CtxSetF | PetscData::CtxSetTime;
  auto f = createDMVector(dm);
  pipeline->f = f;
  pipeline->ts_t = time;
  pipeline->ts_dt = delta_t;
 
  pipeline->ksp_ctx = KspMethod::KSPContext::CTX_SETFUNCTION;
  pipeline->snes_ctx = SnesMethod::SNESContext::CTX_SNESSETFUNCTION;
  pipeline->ts_ctx = TSMethod::TSContext::CTX_TSSETIFUNCTION;
 
  auto [v, a] = setPipelineX(pipeline, x, delta_x, delta2_x, delta_t);
  std::ignore = a;
  std::ignore = v;
 
  CHK_THROW_MESSAGE(
      DMoFEMMeshToLocalVector(dm, x, INSERT_VALUES, SCATTER_REVERSE),
      "loc to vec");
  CHK_THROW_MESSAGE(DMoFEMLoopFiniteElements(dm, fe_name, pipeline, cache_ptr),
                    "Run assemble pipeline");
  CHK_THROW_MESSAGE(VecAssemblyBegin(f), "Asmb vector");
  CHK_THROW_MESSAGE(VecAssemblyEnd(f), "Asmb vector");
 
  CHK_THROW_MESSAGE(VecGhostUpdateBegin(f, ADD_VALUES, SCATTER_REVERSE),
                    "Update vec");
  CHK_THROW_MESSAGE(VecGhostUpdateEnd(f, ADD_VALUES, SCATTER_REVERSE),
                    "Update vec");
  CHK_THROW_MESSAGE(VecGhostUpdateBegin(f, INSERT_VALUES, SCATTER_FORWARD),
                    "Update vec");
  CHK_THROW_MESSAGE(VecGhostUpdateEnd(f, INSERT_VALUES, SCATTER_FORWARD),
                    "Update vec");
 
  auto get_norm = [](auto f) {
    double fnorm;
    VecNorm(f, NORM_2, &fnorm);
    return fnorm;
  };
  MOFEM_LOG_CHANNEL("WORLD");
  MOFEM_TAG_AND_LOG_C("WORLD", Sev::noisy, "OpTester", "Vec norm %3.4e",
                      get_norm(f));
 
  return f;
}

◆ checkCentralFiniteDifference()

SmartPetscObj< Vec > MoFEM::OperatorsTester::checkCentralFiniteDifference	(	SmartPetscObj< DM >	dm,
		std::string	fe_name,
		boost::shared_ptr< FEMethod >	pipeline_rhs,
		boost::shared_ptr< FEMethod >	pipeline_lhs,
		SmartPetscObj< Vec >	x,
		SmartPetscObj< Vec >	delta_x,
		SmartPetscObj< Vec >	delta2_x,
		SmartPetscObj< Vec >	diff_x,
		double	time,
		double	delta_t,
		double	eps,
		CacheTupleWeakPtr	cache_ptr = `CacheTupleSharedPtr()`
	)

Check consistency between directional derivative with matrix.

Parameters

dm
fe_name
pipeline_rhs
pipeline_lhs
x
delta_x
delta2_x
diff_x
time
delta_t
eps
cache_ptr

Returns: SmartPetscObj<Vec>

Definition at line 251 of file OperatorsTester.cpp.

                                 {
 
  auto fd_diff = directionalCentralFiniteDifference(
      dm, fe_name, pipeline_rhs, x, delta_x, delta2_x, diff_x, time, delta_t,
      eps, cache_ptr);
 
  auto m = assembleMat(dm, fe_name, pipeline_lhs, x, delta_x, delta2_x, time,
                       delta_t, cache_ptr);
 
  auto fm = vectorDuplicate(fd_diff);
  CHK_THROW_MESSAGE(MatMult(m, diff_x, fm), "Mat mult");
  CHK_THROW_MESSAGE(VecAXPY(fd_diff, -1., fm), "Add");
 
  return fd_diff;
}

◆ directionalCentralFiniteDifference()

SmartPetscObj< Vec > MoFEM::OperatorsTester::directionalCentralFiniteDifference	(	SmartPetscObj< DM >	dm,
		std::string	fe_name,
		boost::shared_ptr< FEMethod >	pipeline,
		SmartPetscObj< Vec >	x,
		SmartPetscObj< Vec >	delta_x,
		SmartPetscObj< Vec >	delta2_x,
		SmartPetscObj< Vec >	diff_x,
		double	time,
		double	delta_t,
		double	eps,
		CacheTupleWeakPtr	cache_ptr = `CacheTupleSharedPtr()`
	)

Calculate directional directive using finite difference.

Parameters

dm
fe_name
pipeline
x
delta_x
delta2_x
diff_x	// direction of derivative
time
delta_t
eps
cache_ptr

Returns: SmartPetscObj<Vec>

Definition at line 194 of file OperatorsTester.cpp.

                                 {
 
  auto axpy = [&](auto v0, auto diff_v, auto eps) {
    SmartPetscObj<Vec> v;
    if (v0.use_count()) {
      v = vectorDuplicate(v0);
      CHK_THROW_MESSAGE(VecCopy(v0, v), "Cpy");
      CHK_THROW_MESSAGE(VecAXPY(v, eps, diff_v), "Add");
      CHK_THROW_MESSAGE(VecGhostUpdateBegin(v, INSERT_VALUES, SCATTER_FORWARD),
                        "ghost");
      CHK_THROW_MESSAGE(VecGhostUpdateEnd(v, INSERT_VALUES, SCATTER_FORWARD),
                        "ghost");
    }
   return v;
  };
 
  auto fm1 = assembleVec(dm, fe_name, pipeline,
                         //
                         axpy(x, diff_x, -eps),
                         //
                         axpy(delta_x, diff_x, -eps),
                         //
                         axpy(delta2_x, diff_x, -eps),
                         //
                         time, delta_t, cache_ptr);
 
  auto fp1 = assembleVec(dm, fe_name, pipeline,
                         //
                         axpy(x, diff_x, eps),
                         //
                         axpy(delta_x, diff_x, eps),
                         //
                         axpy(delta2_x, diff_x, eps),
                         //
                         time, delta_t, cache_ptr);
 
  auto calculate_derivative = [eps](auto fm1, auto fp1) {
    int size;
    VecGetLocalSize(fp1, &size);
    VecAXPY(fp1, -1, fm1);
    double *array;
    VecGetArray(fp1, &array);
    for (int j = 0; j != size; ++j) {
      array[j] /= 2 * eps;
    }
    VecRestoreArray(fp1, &array);
    return fp1;
  };
 
  return calculate_derivative(fm1, fp1);
}

◆ query_interface()

MoFEMErrorCode MoFEM::OperatorsTester::query_interface	(	boost::typeindex::type_index	type_index,
		UnknownInterface **	iface
	)		const

virtual

Implements MoFEM::UnknownInterface.

Definition at line 17 of file OperatorsTester.cpp.

                                                                 {
  *iface = const_cast<OperatorsTester *>(this);
  return 0;
}

◆ setPipelineX()

std::pair< SmartPetscObj< Vec >, SmartPetscObj< Vec > > MoFEM::OperatorsTester::setPipelineX	(	boost::shared_ptr< FEMethod >	pipeline,
		SmartPetscObj< Vec >	x,
		SmartPetscObj< Vec >	delta_x,
		SmartPetscObj< Vec >	delta2_x,
		double	delta_t
	)

private

Set vectors, x, x_t, and x_tt to finite element instance.

Finite element instance is a pipeline. x_t and x_tt are evaluated for given delta_x, delta2_x and delta_t.

Parameters

pipeline
x
delta_x
delta2_x
delta_t

Returns: std::pair<SmartPetscObj<Vec>, SmartPetscObj<Vec>>

Definition at line 274 of file OperatorsTester.cpp.

                                                                           {
 
  // Set dofs vector to finite element instance
  pipeline->x = x;
  pipeline->data_ctx |= PetscData::CTX_SET_X;
 
  SmartPetscObj<Vec> x_t, x_tt;
 
  // Set velocity dofs vector to finite element instance
 
  if (delta_x.use_count()) {
    x_t = vectorDuplicate(x);
    VecCopy(delta_x, x_t);
    VecScale(x_t, 1. / delta_t);
    CHK_THROW_MESSAGE(VecGhostUpdateBegin(x_t, INSERT_VALUES, SCATTER_FORWARD),
                      "ghost");
    CHK_THROW_MESSAGE(VecGhostUpdateEnd(x_t, INSERT_VALUES, SCATTER_FORWARD),
                      "ghost");
    pipeline->x_t = x_t;
    pipeline->data_ctx |= PetscData::CTX_SET_X_T;
  } else {
    pipeline->x_t = PETSC_NULLPTR;
  }
 
  // Set acceleration dofs vector to finite element instance
 
  if (delta2_x.use_count()) {
    x_tt = vectorDuplicate(x);
    VecCopy(delta2_x, x_tt);
    VecScale(x_tt, 1. / pow(delta_t, 2));
    CHK_THROW_MESSAGE(VecGhostUpdateBegin(x_tt, INSERT_VALUES, SCATTER_FORWARD),
                      "ghost");
    CHK_THROW_MESSAGE(VecGhostUpdateEnd(x_tt, INSERT_VALUES, SCATTER_FORWARD),
                      "ghost");
    pipeline->x_tt = x_tt;
    pipeline->data_ctx |= PetscData::CTX_SET_X_TT;
  } else {
    pipeline->x_tt = PETSC_NULLPTR;
  }
 
  return std::make_pair(x_t, x_tt);
}

◆ setRandomFields()

SmartPetscObj< Vec > MoFEM::OperatorsTester::setRandomFields	(	SmartPetscObj< DM >	dm,
		std::vector< RandomFieldData >	random_fields,
		boost::shared_ptr< Range >	ents = `nullptr`,
		RowColData	r = `RowColData::COL`
	)

Generate random fields.

Example: generate random vector for DM (problem) from simple interface, where FIELD random values of DOFs are in range from -1 to 1, and FIELD2 random values are in range from 0 to 1.

auto x = opt->setRandomFields(simple->getDM(),

{{"FIELD1", {-1, 1}}, {"FIELD2", {0,1}}});

simple

void simple(double P1[], double P2[], double P3[], double c[], const int N)

Definition acoustic.cpp:69

TODO: Set random field to specific entities, and potentially order for testing proposes to dissect error in tangent matrix.

Parameters

dm
random_fields	look at definition RandomFieldData

Returns: SmartPetscObj<Vec> smart vector

Examples: mofem/atom_tests/tensor_divergence_operator.cpp.

Definition at line 23 of file OperatorsTester.cpp.

                                                   {
 
  MoFEM::Interface &m_field = cOre;
 
  auto prb_ptr = getProblemPtr(dm);
 
  auto get_random_number = [](auto &range) {
    const auto r = static_cast<double>(std::rand()) / RAND_MAX;
    return range[0] + r * (range[1] - range[0]);
  };
 
  auto v = createDMVector(dm, r);
  double *array;
  CHKERR VecGetArray(v, &array);
 
  int size;
  VecGetLocalSize(v, &size);
 
  MOFEM_LOG_CHANNEL("WORLD");
 
  for (auto &rf : random_fields) {
    const auto field_id = m_field.get_field_bit_number(rf.first);
    MOFEM_TAG_AND_LOG("WORLD", Sev::noisy, "OperatorsTester")
        << "Set random field " << rf.first << " " << static_cast<int>(field_id);
 
    auto get_numbered_dofs = [&]() {
      switch (r) {
      case RowColData::ROW:
        return prb_ptr->getNumeredRowDofsPtr();
      case RowColData::COL:
        return prb_ptr->getNumeredColDofsPtr();
      default:
        throw std::runtime_error("Invalid RowColData");
      }
    };
 
    if (!ents_ptr) {
      auto lo = get_numbered_dofs()->lower_bound(
          FieldEntity::getLoBitNumberUId(field_id));
      const auto hi = get_numbered_dofs()->upper_bound(
          FieldEntity::getHiBitNumberUId(field_id));
      for (; lo != hi; ++lo) {
        const auto idx = (*lo)->getPetscLocalDofIdx();
        if (idx >= 0 && idx < size) {
          array[idx] = get_random_number(rf.second);
        }
      }
    } else {
      for (auto pit = ents_ptr->const_pair_begin();
           pit != ents_ptr->const_pair_end(); ++pit) {
        auto lo = get_numbered_dofs()->get<Unique_mi_tag>().lower_bound(
            DofEntity::getLoFieldEntityUId(field_id, pit->first));
        auto hi = get_numbered_dofs()->get<Unique_mi_tag>().upper_bound(
            DofEntity::getHiFieldEntityUId(field_id, pit->second));
        for (; lo != hi; ++lo) {
          const auto idx = (*lo)->getPetscLocalDofIdx();
          if (idx >= 0 && idx < size) {
            array[idx] = get_random_number(rf.second);
          }
        }
      }
    }
  }
 
  CHKERR VecRestoreArray(v, &array);
 
  CHKERR VecGhostUpdateBegin(v, INSERT_VALUES, SCATTER_FORWARD);
  CHKERR VecGhostUpdateEnd(v, INSERT_VALUES, SCATTER_FORWARD);
 
  return v;
};

Member Data Documentation

◆ cOre

MoFEM::Core& MoFEM::OperatorsTester::cOre

private

Definition at line 146 of file OperatorsTester.hpp.

The documentation for this struct was generated from the following files:

src/interfaces/OperatorsTester.hpp
src/interfaces/impl/OperatorsTester.cpp

Public Types

Public Member Functions

Private Member Functions

Private Attributes

Additional Inherited Members

Detailed Description

Member Typedef Documentation

◆ RandomFieldData

Constructor & Destructor Documentation

◆ OperatorsTester()

◆ ~OperatorsTester()

Member Function Documentation

◆ assembleMat()

◆ assembleVec()

◆ checkCentralFiniteDifference()

◆ directionalCentralFiniteDifference()

◆ query_interface()

◆ setPipelineX()

◆ setRandomFields()

Member Data Documentation

◆ cOre