Skip to content

Commit

Permalink
Add coupling validation routines (idaholab#5595)
Browse files Browse the repository at this point in the history
  • Loading branch information
dschwen committed Aug 28, 2015
1 parent a21b6bf commit dc04a01
Show file tree
Hide file tree
Showing 17 changed files with 216 additions and 8 deletions.
138 changes: 137 additions & 1 deletion framework/include/materials/DerivativeMaterialInterface.h
Expand Up @@ -16,6 +16,7 @@

#include "Material.h"
#include "MaterialProperty.h"
#include "KernelBase.h"
#include "FEProblem.h"
#include "BlockRestrictable.h"
#include "BoundaryRestrictable.h"
Expand All @@ -37,7 +38,7 @@ void mooseSetToZero(T & v)
v = 0;
}
template<typename T>
void mooseSetToZero(T* & )
void mooseSetToZero(T* &)
{
mooseError("Cannot use pointer types for MaterialProperty derivatives.");
}
Expand Down Expand Up @@ -104,6 +105,19 @@ class DerivativeMaterialInterface :
const MaterialProperty<U> & getMaterialPropertyDerivativeByName(const MaterialPropertyName & base, const VariableName & c1, const VariableName & c2 = "", const VariableName & c3 = "");
///@}

///@{
/**
* check if derivatives of the passed in material property exist w.r.t a variable
* that is _not_ coupled in to the current object
*/
template<typename U>
void validateCoupling(const MaterialPropertyName & base, const std::vector<VariableName> & c, bool validate_aux = true);
template<typename U>
void validateCoupling(const MaterialPropertyName & base, const VariableName & c1 = "", const VariableName & c2 = "", const VariableName & c3 = "");
template<typename U>
void validateNonlinearCoupling(const MaterialPropertyName & base, const VariableName & c1 = "", const VariableName & c2 = "", const VariableName & c3 = "");
///@}

private:
/// Return a constant zero property
template<typename U>
Expand All @@ -113,6 +127,16 @@ class DerivativeMaterialInterface :
template<typename U>
bool haveMaterialProperty(const std::string & prop_name);

/// helper method to combine multiple VariableNames into a vector (if they are != "")
std::vector<VariableName> buildVariableVector(const VariableName & c1, const VariableName & c2, const VariableName & c3);

/// helper method to compile list of missing coupled variables for a given system
template<typename U>
void validateCouplingHelper(const MaterialPropertyName & base, const std::vector<VariableName> & c, const System & system, std::vector<VariableName> & missing);

// check if the speciified variable name is not the variable this kernel is acting on (always true for any other type of object)
bool isNotKernelVariable(const VariableName & name);

/// Reference to FEProblem
FEProblem & _dmi_fe_problem;

Expand Down Expand Up @@ -305,4 +329,116 @@ DerivativeMaterialInterface<T>::getMaterialPropertyDerivativeByName(const Materi
return getDefaultMaterialPropertyByName<U>(propertyNameFirst(base, c1));
}

template<class T>
template<typename U>
void
DerivativeMaterialInterface<T>::validateCouplingHelper(const MaterialPropertyName & base, const std::vector<VariableName> & c, const System & system, std::vector<VariableName> & missing)
{
unsigned int ncoupled = this->_coupled_moose_vars.size();

// iterate over all variables in the current system (in groups)
for (unsigned int i = 0; i < system.n_variable_groups(); ++i)
{
const VariableGroup & vg = system.variable_group(i);
for (unsigned int j = 0; j < vg.n_variables(); ++j)
{
std::vector<VariableName> cj(c);
VariableName jname = vg.name(j);
cj.push_back(jname);

// if the derivative exists make sure the variable is coupled
if (haveMaterialProperty<U>(propertyName(base, cj)))
{
// kernels to not have the variable they are acting on in coupled_moose_vars
bool is_missing = isNotKernelVariable(jname);

for (unsigned int k = 0; k < ncoupled; ++k)
if (this->_coupled_moose_vars[k]->name() == jname)
{
is_missing = false;
break;
}

if (is_missing)
missing.push_back(jname);
}
}
}
}

template<class T>
template<typename U>
void
DerivativeMaterialInterface<T>::validateCoupling(const MaterialPropertyName & base, const std::vector<VariableName> & c, bool validate_aux)
{
// get the base property name
std::string prop_name = this->deducePropertyName(base);
// list of potentially missing coupled variables
std::vector<VariableName> missing;

// iterate over all variables in the both the non-linear and auxiliary system (optional)
validateCouplingHelper<U>(prop_name, c, _dmi_fe_problem.getNonlinearSystem().system(), missing);
if (validate_aux)
validateCouplingHelper<U>(prop_name, c, _dmi_fe_problem.getAuxiliarySystem().system(), missing);

if (missing.size() > 0)
{
// join list of missing variable names
std::string list = missing[0];
for (unsigned int i = 1; i < missing.size(); ++i)
list += ", " + missing[i];

mooseError("Missing coupled variables {" << list << "} (add them to args parameter of " << this->name() << ")");
}
}

template<class T>
std::vector<VariableName>
DerivativeMaterialInterface<T>::buildVariableVector(const VariableName & c1, const VariableName & c2, const VariableName & c3)
{
std::vector<VariableName> c;
if (c1 != "")
{
c.push_back(c1);
if (c2 != "")
{
c.push_back(c2);
if (c3 != "")
c.push_back(c3);
}
}
return c;
}

template<class T>
template<typename U>
void
DerivativeMaterialInterface<T>::validateCoupling(const MaterialPropertyName & base, const VariableName & c1, const VariableName & c2, const VariableName & c3)
{
validateCoupling<U>(base, buildVariableVector(c1, c2, c3), true);
}

template<class T>
template<typename U>
void
DerivativeMaterialInterface<T>::validateNonlinearCoupling(const MaterialPropertyName & base, const VariableName & c1, const VariableName & c2, const VariableName & c3)
{
validateCoupling<U>(base, buildVariableVector(c1, c2, c3), false);
}

template<class T>
inline bool
DerivativeMaterialInterface<T>::isNotKernelVariable(const VariableName & name)
{
// try to cast this to a Kernel pointer
KernelBase * k = dynamic_cast<KernelBase *>(this);

// This interface is not templated on a class derived from Kernel
if (k == NULL)
return true;

// We are templated on a kernel class, so we check if the kernel variable
return k->variable().name() != name;
}

#endif //DERIVATIVEMATERIALINTERFACE_H
2 changes: 2 additions & 0 deletions modules/phase_field/include/kernels/ACBulk.h
Expand Up @@ -21,6 +21,8 @@ class ACBulk : public DerivativeMaterialInterface<JvarMapInterface<KernelValue>
public:
ACBulk(const InputParameters & parameters);

virtual void initialSetup();

protected:

/// Enum used with computeDFDOP function
Expand Down
2 changes: 2 additions & 0 deletions modules/phase_field/include/kernels/ACInterface.h
Expand Up @@ -21,6 +21,8 @@ class ACInterface : public DerivativeMaterialInterface<JvarMapInterface<KernelGr
public:
ACInterface(const InputParameters & parameters);

virtual void initialSetup();

protected:

/// Enum of computeDFDOP inputs
Expand Down
2 changes: 2 additions & 0 deletions modules/phase_field/include/kernels/ACParsed.h
Expand Up @@ -25,6 +25,8 @@ class ACParsed : public ACBulk
public:
ACParsed(const InputParameters & parameters);

virtual void initialSetup();

protected:
virtual Real computeDFDOP(PFFunctionType type);
virtual Real computeQpOffDiagJacobian(unsigned int jvar);
Expand Down
13 changes: 13 additions & 0 deletions modules/phase_field/include/kernels/CahnHilliardBase.h
Expand Up @@ -21,6 +21,7 @@ class CahnHilliardBase : public CHBulk<T>
CahnHilliardBase(const InputParameters & parameters);

static InputParameters validParams();
virtual void initialSetup();

protected:
typedef typename CHBulk<T>::PFFunctionType PFFunctionType;
Expand Down Expand Up @@ -77,6 +78,18 @@ CahnHilliardBase<T>::CahnHilliardBase(const InputParameters & parameters) :
}
}

template<typename T>
void
CahnHilliardBase<T>:: initialSetup()
{
/**
* Check if both the non-linear as well as the auxiliary variables variables
* are coupled. Derivatives with respect to both types of variables contribute
* the residual.
*/
this->template validateCoupling<Real>("f_name", this->_var.name());
}

template<typename T>
RealGradient
CahnHilliardBase<T>::computeGradDFDCons(PFFunctionType type)
Expand Down
2 changes: 2 additions & 0 deletions modules/phase_field/include/kernels/SplitCHParsed.h
Expand Up @@ -28,6 +28,8 @@ class SplitCHParsed : public DerivativeMaterialInterface<JvarMapInterface<SplitC
public:
SplitCHParsed(const InputParameters & parameters);

virtual void initialSetup();

protected:
virtual Real computeDFDC(PFFunctionType type);
virtual Real computeQpOffDiagJacobian(unsigned int jvar);
Expand Down
Expand Up @@ -26,6 +26,8 @@ class DerivativeMultiPhaseBase : public DerivativeFunctionMaterialBase
public:
DerivativeMultiPhaseBase(const InputParameters & parameters);

virtual void initialSetup();

protected:
virtual Real computeF();

Expand Down
2 changes: 2 additions & 0 deletions modules/phase_field/include/materials/DerivativeSumMaterial.h
Expand Up @@ -19,6 +19,8 @@ class DerivativeSumMaterial : public DerivativeFunctionMaterialBase
public:
DerivativeSumMaterial(const InputParameters & parameters);

virtual void initialSetup();

protected:
virtual void computeProperties();

Expand Down
Expand Up @@ -26,6 +26,8 @@ class DerivativeTwoPhaseMaterial : public DerivativeFunctionMaterialBase
public:
DerivativeTwoPhaseMaterial(const InputParameters & parameters);

virtual void initialSetup();

protected:
virtual Real computeF();
virtual Real computeDF(unsigned int);
Expand Down
7 changes: 6 additions & 1 deletion modules/phase_field/src/kernels/ACBulk.C
Expand Up @@ -32,6 +32,12 @@ ACBulk::ACBulk(const InputParameters & parameters) :
_dLdarg[i] = &getMaterialPropertyDerivative<Real>("mob_name", _coupled_moose_vars[i]->name());
}

void
ACBulk::initialSetup()
{
validateNonlinearCoupling<Real>("mob_name");
}

Real
ACBulk::precomputeQpResidual()
{
Expand Down Expand Up @@ -65,4 +71,3 @@ ACBulk::computeQpOffDiagJacobian(unsigned int jvar)
// Set off-diagonal Jacobian term from mobility derivatives
return (*_dLdarg[cvar])[_qp] * _phi[_j][_qp] * computeDFDOP(Residual) * _test[_i][_qp];
}

7 changes: 6 additions & 1 deletion modules/phase_field/src/kernels/ACInterface.C
Expand Up @@ -34,6 +34,12 @@ ACInterface::ACInterface(const InputParameters & parameters) :
_dLdarg[i] = &getMaterialPropertyDerivative<Real>("mob_name", _coupled_moose_vars[i]->name());
}

void
ACInterface::initialSetup()
{
validateNonlinearCoupling<Real>("mob_name");
}

RealGradient
ACInterface::precomputeQpResidual()
{
Expand All @@ -59,4 +65,3 @@ ACInterface::computeQpOffDiagJacobian(unsigned int jvar)
// Set off-diagonal jaocbian terms from mobility dependence
return _kappa[_qp] * (*_dLdarg[cvar])[_qp] * _phi[_j][_qp] * _grad_u[_qp] * _grad_test[_i][_qp];
}

8 changes: 7 additions & 1 deletion modules/phase_field/src/kernels/ACParsed.C
Expand Up @@ -29,6 +29,13 @@ ACParsed::ACParsed(const InputParameters & parameters) :
_d2FdEtadarg[i] = &getMaterialPropertyDerivative<Real>("f_name", _var.name(), _coupled_moose_vars[i]->name());
}

void
ACParsed::initialSetup()
{
ACBulk::initialSetup();
validateNonlinearCoupling<Real>("f_name");
}

Real
ACParsed::computeDFDOP(PFFunctionType type)
{
Expand All @@ -55,4 +62,3 @@ ACParsed::computeQpOffDiagJacobian(unsigned int jvar)
return ACBulk::computeQpOffDiagJacobian(jvar) +
_L[_qp] * (*_d2FdEtadarg[cvar])[_qp] * _phi[_j][_qp] * _test[_i][_qp];
}

12 changes: 11 additions & 1 deletion modules/phase_field/src/kernels/SplitCHParsed.C
Expand Up @@ -30,6 +30,17 @@ SplitCHParsed::SplitCHParsed(const InputParameters & parameters) :
_d2Fdcdarg[i] = &getMaterialPropertyDerivative<Real>("f_name", _var.name(), _coupled_moose_vars[i]->name());
}

void
SplitCHParsed::initialSetup()
{
/**
* We are only interested if the necessary non-linear variables are coupled,
* as those are thge ones used in constructing the Jacobian. The AuxVariables
* do not have Jacobian entries.
*/
validateNonlinearCoupling<Real>("f_name", _var.name());
}

Real
SplitCHParsed::computeDFDC(PFFunctionType type)
{
Expand Down Expand Up @@ -58,4 +69,3 @@ SplitCHParsed::computeQpOffDiagJacobian(unsigned int jvar)

return (*_d2Fdcdarg[cvar])[_qp] * _phi[_j][_qp] * _test[_i][_qp];
}

8 changes: 7 additions & 1 deletion modules/phase_field/src/materials/DerivativeMultiPhaseBase.C
Expand Up @@ -115,6 +115,13 @@ DerivativeMultiPhaseBase::DerivativeMultiPhaseBase(const InputParameters & param
}
}

void
DerivativeMultiPhaseBase::initialSetup()
{
for (unsigned int n = 0; n < _num_fi; ++n)
validateCoupling<Real>(_fi_names[n]);
}

Real
DerivativeMultiPhaseBase::computeF()
{
Expand All @@ -123,4 +130,3 @@ DerivativeMultiPhaseBase::computeF()
F += (*_hi[n])[_qp] * (*_prop_Fi[n])[_qp];
return F + _W * _g[_qp];
}

8 changes: 7 additions & 1 deletion modules/phase_field/src/materials/DerivativeSumMaterial.C
Expand Up @@ -81,6 +81,13 @@ DerivativeSumMaterial::DerivativeSumMaterial(const InputParameters & parameters)
}
}

void
DerivativeSumMaterial::initialSetup()
{
for (unsigned int n = 0; n < _num_materials; ++n)
validateCoupling<Real>(_sum_materials[n]);
}

void
DerivativeSumMaterial::computeProperties()
{
Expand Down Expand Up @@ -127,4 +134,3 @@ DerivativeSumMaterial::computeProperties()
}
}
}

Expand Up @@ -87,6 +87,13 @@ DerivativeTwoPhaseMaterial::DerivativeTwoPhaseMaterial(const InputParameters & p
}
}

void
DerivativeTwoPhaseMaterial::initialSetup()
{
validateCoupling<Real>("fa_name");
validateCoupling<Real>("fb_name");
}

Real
DerivativeTwoPhaseMaterial::computeF()
{
Expand Down Expand Up @@ -148,4 +155,3 @@ DerivativeTwoPhaseMaterial::computeD3F(unsigned int i_var, unsigned int j_var, u

return _h[_qp] * (*_prop_d3Fb[i][j][k])[_qp] + (1.0 - _h[_qp]) * (*_prop_d3Fa[i][j][k])[_qp];
}

0 comments on commit dc04a01

Please sign in to comment.