www.mooseframework.org
Public Member Functions | Static Public Member Functions | Protected Member Functions | Protected Attributes | Private Member Functions | List of all members
MultiPlasticityDebugger Class Reference

MultiPlasticityDebugger computes various finite-difference things to help developers remove bugs in their derivatives, etc. More...

#include <MultiPlasticityDebugger.h>

Inheritance diagram for MultiPlasticityDebugger:
[legend]

Public Member Functions

 MultiPlasticityDebugger (const MooseObject *moose_object)
 
void outputAndCheckDebugParameters ()
 Outputs the debug parameters: _fspb_debug_stress, _fspd_debug_pm, etc and checks that they are sized correctly. More...
 
void checkDerivatives ()
 Checks the derivatives, eg dyieldFunction_dstress by using finite difference approximations. More...
 
void checkJacobian (const RankFourTensor &E_inv, const std::vector< Real > &intnl_old)
 Checks the full Jacobian, which is just certain linear combinations of the dyieldFunction_dstress, etc, by using finite difference approximations. More...
 
void checkSolution (const RankFourTensor &E_inv)
 Checks that Ax does equal b in the NR procedure. More...
 
UserObjectName getUserObjectName (const std::string &param_name) const
 
const T & getUserObject (const std::string &param_name, bool is_dependency=true) const
 
const T & getUserObjectByName (const UserObjectName &object_name, bool is_dependency=true) const
 
const UserObjectgetUserObjectBase (const std::string &param_name, bool is_dependency=true) const
 
const UserObjectgetUserObjectBaseByName (const UserObjectName &object_name, bool is_dependency=true) const
 
bool hasUserObject (const std::string &param_name) const
 
bool hasUserObject (const std::string &param_name) const
 
bool hasUserObject (const std::string &param_name) const
 
bool hasUserObject (const std::string &param_name) const
 
bool hasUserObjectByName (const UserObjectName &object_name) const
 
bool hasUserObjectByName (const UserObjectName &object_name) const
 
bool hasUserObjectByName (const UserObjectName &object_name) const
 
bool hasUserObjectByName (const UserObjectName &object_name) const
 

Static Public Member Functions

static InputParameters validParams ()
 

Protected Member Functions

virtual void calculateConstraints (const RankTwoTensor &stress, const std::vector< Real > &intnl_old, const std::vector< Real > &intnl, const std::vector< Real > &pm, const RankTwoTensor &delta_dp, std::vector< Real > &f, std::vector< RankTwoTensor > &r, RankTwoTensor &epp, std::vector< Real > &ic, const std::vector< bool > &active)
 The constraints. More...
 
virtual void calculateRHS (const RankTwoTensor &stress, const std::vector< Real > &intnl_old, const std::vector< Real > &intnl, const std::vector< Real > &pm, const RankTwoTensor &delta_dp, std::vector< Real > &rhs, const std::vector< bool > &active, bool eliminate_ld, std::vector< bool > &deactivated_due_to_ld)
 Calculate the RHS which is rhs = -(epp(0,0), epp(1,0), epp(1,1), epp(2,0), epp(2,1), epp(2,2), f[0], f[1], ..., f[num_f], ic[0], ic[1], ..., ic[num_ic]) More...
 
virtual void calculateJacobian (const RankTwoTensor &stress, const std::vector< Real > &intnl, const std::vector< Real > &pm, const RankFourTensor &E_inv, const std::vector< bool > &active, const std::vector< bool > &deactivated_due_to_ld, std::vector< std::vector< Real >> &jac)
 d(rhs)/d(dof) More...
 
virtual void nrStep (const RankTwoTensor &stress, const std::vector< Real > &intnl_old, const std::vector< Real > &intnl, const std::vector< Real > &pm, const RankFourTensor &E_inv, const RankTwoTensor &delta_dp, RankTwoTensor &dstress, std::vector< Real > &dpm, std::vector< Real > &dintnl, const std::vector< bool > &active, std::vector< bool > &deactivated_due_to_ld)
 Performs one Newton-Raphson step. More...
 
virtual void yieldFunction (const RankTwoTensor &stress, const std::vector< Real > &intnl, const std::vector< bool > &active, std::vector< Real > &f)
 The active yield function(s) More...
 
virtual void dyieldFunction_dstress (const RankTwoTensor &stress, const std::vector< Real > &intnl, const std::vector< bool > &active, std::vector< RankTwoTensor > &df_dstress)
 The derivative of the active yield function(s) with respect to stress. More...
 
virtual void dyieldFunction_dintnl (const RankTwoTensor &stress, const std::vector< Real > &intnl, const std::vector< bool > &active, std::vector< Real > &df_dintnl)
 The derivative of active yield function(s) with respect to their internal parameters (the user objects assume there is exactly one internal param per yield function) More...
 
virtual void flowPotential (const RankTwoTensor &stress, const std::vector< Real > &intnl, const std::vector< bool > &active, std::vector< RankTwoTensor > &r)
 The active flow potential(s) - one for each yield function. More...
 
virtual void dflowPotential_dstress (const RankTwoTensor &stress, const std::vector< Real > &intnl, const std::vector< bool > &active, std::vector< RankFourTensor > &dr_dstress)
 The derivative of the active flow potential(s) with respect to stress. More...
 
virtual void dflowPotential_dintnl (const RankTwoTensor &stress, const std::vector< Real > &intnl, const std::vector< bool > &active, std::vector< RankTwoTensor > &dr_dintnl)
 The derivative of the active flow potentials with respect to the active internal parameters The UserObjects explicitly assume that r[alpha] is only dependent on intnl[alpha]. More...
 
virtual void hardPotential (const RankTwoTensor &stress, const std::vector< Real > &intnl, const std::vector< bool > &active, std::vector< Real > &h)
 The active hardening potentials (one for each internal parameter and for each yield function) by assumption in the Userobjects, the h[a][alpha] is nonzero only if the surface alpha is part of model a, so we only calculate those here. More...
 
virtual void dhardPotential_dstress (const RankTwoTensor &stress, const std::vector< Real > &intnl, const std::vector< bool > &active, std::vector< RankTwoTensor > &dh_dstress)
 The derivative of the active hardening potentials with respect to stress By assumption in the Userobjects, the h[a][alpha] is nonzero only for a = alpha, so we only calculate those here. More...
 
virtual void dhardPotential_dintnl (const RankTwoTensor &stress, const std::vector< Real > &intnl, const std::vector< bool > &active, std::vector< Real > &dh_dintnl)
 The derivative of the active hardening potentials with respect to the active internal parameters. More...
 
virtual void buildActiveConstraints (const std::vector< Real > &f, const RankTwoTensor &stress, const std::vector< Real > &intnl, const RankFourTensor &Eijkl, std::vector< bool > &act)
 Constructs a set of active constraints, given the yield functions, f. More...
 
unsigned int modelNumber (unsigned int surface)
 returns the model number, given the surface number More...
 
bool anyActiveSurfaces (int model, const std::vector< bool > &active)
 returns true if any internal surfaces of the given model are active according to 'active' More...
 
void activeModelSurfaces (int model, const std::vector< bool > &active, std::vector< unsigned int > &active_surfaces_of_model)
 Returns the internal surface number(s) of the active surfaces of the given model This may be of size=0 if there are no active surfaces of the given model. More...
 
void activeSurfaces (int model, const std::vector< bool > &active, std::vector< unsigned int > &active_surfaces)
 Returns the external surface number(s) of the active surfaces of the given model This may be of size=0 if there are no active surfaces of the given model. More...
 
bool returnMapAll (const RankTwoTensor &trial_stress, const std::vector< Real > &intnl_old, const RankFourTensor &E_ijkl, Real ep_plastic_tolerance, RankTwoTensor &stress, std::vector< Real > &intnl, std::vector< Real > &pm, std::vector< Real > &cumulative_pm, RankTwoTensor &delta_dp, std::vector< Real > &yf, unsigned &num_successful_plastic_returns, unsigned &custom_model)
 Performs a returnMap for each plastic model using their inbuilt returnMap functions. More...
 
virtual void addUserObjectDependencyHelper (const UserObject &) const
 

Protected Attributes

MooseEnum _fspb_debug
 none - don't do any debugging crash - currently inactive jacobian - check the jacobian entries jacobian_and_linear_system - check entire jacobian and check that Ax=b More...
 
RankTwoTensor _fspb_debug_stress
 Debug the Jacobian entries at this stress. More...
 
std::vector< Real_fspb_debug_pm
 Debug the Jacobian entires at these plastic multipliers. More...
 
std::vector< Real_fspb_debug_intnl
 Debug the Jacobian entires at these internal parameters. More...
 
Real _fspb_debug_stress_change
 Debug finite-differencing parameter for the stress. More...
 
std::vector< Real_fspb_debug_pm_change
 Debug finite-differencing parameters for the plastic multipliers. More...
 
std::vector< Real_fspb_debug_intnl_change
 Debug finite-differencing parameters for the internal parameters. More...
 
Real _svd_tol
 Tolerance on the minimum ratio of singular values before flow-directions are deemed linearly dependent. More...
 
Real _min_f_tol
 Minimum value of the _f_tol parameters for the Yield Function User Objects. More...
 
const InputParameters_params
 
unsigned int _num_models
 Number of plastic models for this material. More...
 
unsigned int _num_surfaces
 Number of surfaces within the plastic models. More...
 
std::vector< std::vector< unsigned int > > _surfaces_given_model
 _surfaces_given_model[model_number] = vector of surface numbers for this model More...
 
MooseEnum _specialIC
 Allows initial set of active constraints to be chosen optimally. More...
 
std::vector< const SolidMechanicsPlasticModel * > _f
 User objects that define the yield functions, flow potentials, etc. More...
 

Private Member Functions

void fddyieldFunction_dstress (const RankTwoTensor &stress, const std::vector< Real > &intnl, std::vector< RankTwoTensor > &df_dstress)
 The finite-difference derivative of yield function(s) with respect to stress. More...
 
void fddyieldFunction_dintnl (const RankTwoTensor &stress, const std::vector< Real > &intnl, std::vector< Real > &df_dintnl)
 The finite-difference derivative of yield function(s) with respect to internal parameter(s) More...
 
virtual void fddflowPotential_dstress (const RankTwoTensor &stress, const std::vector< Real > &intnl, std::vector< RankFourTensor > &dr_dstress)
 The finite-difference derivative of the flow potential(s) with respect to stress. More...
 
virtual void fddflowPotential_dintnl (const RankTwoTensor &stress, const std::vector< Real > &intnl, std::vector< RankTwoTensor > &dr_dintnl)
 The finite-difference derivative of the flow potentials with respect to internal parameters. More...
 
virtual void fdJacobian (const RankTwoTensor &stress, const std::vector< Real > &intnl_old, const std::vector< Real > &intnl, const std::vector< Real > &pm, const RankTwoTensor &delta_dp, const RankFourTensor &E_inv, bool eliminate_ld, std::vector< std::vector< Real >> &jac)
 The Jacobian calculated using finite differences. More...
 
bool dof_included (unsigned int dof, const std::vector< bool > &deactivated_due_to_ld)
 

Detailed Description

MultiPlasticityDebugger computes various finite-difference things to help developers remove bugs in their derivatives, etc.

Definition at line 19 of file MultiPlasticityDebugger.h.

Constructor & Destructor Documentation

◆ MultiPlasticityDebugger()

MultiPlasticityDebugger::MultiPlasticityDebugger ( const MooseObject moose_object)

Definition at line 40 of file MultiPlasticityDebugger.C.

41  : MultiPlasticityLinearSystem(moose_object),
42  _fspb_debug(_params.get<MooseEnum>("debug_fspb")),
43  _fspb_debug_stress(_params.get<RealTensorValue>("debug_jac_at_stress")),
44  _fspb_debug_pm(_params.get<std::vector<Real>>("debug_jac_at_pm")),
45  _fspb_debug_intnl(_params.get<std::vector<Real>>("debug_jac_at_intnl")),
46  _fspb_debug_stress_change(_params.get<Real>("debug_stress_change")),
47  _fspb_debug_pm_change(_params.get<std::vector<Real>>("debug_pm_change")),
48  _fspb_debug_intnl_change(_params.get<std::vector<Real>>("debug_intnl_change"))
49 {
50 }
MooseEnum _fspb_debug
none - don&#39;t do any debugging crash - currently inactive jacobian - check the jacobian entries jacobi...
std::vector< std::pair< R1, R2 > > get(const std::string &param1, const std::string &param2) const
Real _fspb_debug_stress_change
Debug finite-differencing parameter for the stress.
std::vector< Real > _fspb_debug_intnl_change
Debug finite-differencing parameters for the internal parameters.
std::vector< Real > _fspb_debug_intnl
Debug the Jacobian entires at these internal parameters.
TensorValue< Real > RealTensorValue
MultiPlasticityLinearSystem(const MooseObject *moose_object)
RankTwoTensor _fspb_debug_stress
Debug the Jacobian entries at this stress.
std::vector< Real > _fspb_debug_pm_change
Debug finite-differencing parameters for the plastic multipliers.
std::vector< Real > _fspb_debug_pm
Debug the Jacobian entires at these plastic multipliers.

Member Function Documentation

◆ activeModelSurfaces()

void MultiPlasticityRawComponentAssembler::activeModelSurfaces ( int  model,
const std::vector< bool > &  active,
std::vector< unsigned int > &  active_surfaces_of_model 
)
protectedinherited

Returns the internal surface number(s) of the active surfaces of the given model This may be of size=0 if there are no active surfaces of the given model.

Parameters
modelthe model number
activearray with entries being 'true' if the surface is active
[out]active_surfaces_of_modelthe output

Definition at line 809 of file MultiPlasticityRawComponentAssembler.C.

Referenced by MultiPlasticityRawComponentAssembler::dflowPotential_dintnl(), MultiPlasticityRawComponentAssembler::dflowPotential_dstress(), MultiPlasticityRawComponentAssembler::dhardPotential_dintnl(), MultiPlasticityRawComponentAssembler::dhardPotential_dstress(), MultiPlasticityRawComponentAssembler::dyieldFunction_dintnl(), MultiPlasticityRawComponentAssembler::dyieldFunction_dstress(), MultiPlasticityRawComponentAssembler::flowPotential(), MultiPlasticityRawComponentAssembler::hardPotential(), and MultiPlasticityRawComponentAssembler::yieldFunction().

813 {
814  active_surfaces_of_model.resize(0);
815  for (unsigned model_surface = 0; model_surface < _f[model]->numberSurfaces(); ++model_surface)
816  if (active[_surfaces_given_model[model][model_surface]])
817  active_surfaces_of_model.push_back(model_surface);
818 }
std::vector< std::vector< unsigned int > > _surfaces_given_model
_surfaces_given_model[model_number] = vector of surface numbers for this model
const PertinentGeochemicalSystem model(database, {"H2O", "H+", "HCO3-", "O2(aq)", "Ca++", ">(s)FeOH", "radius_neg1", "radius_neg1.5"}, {"Calcite"}, {}, {"Calcite_asdf"}, {"CH4(aq)"}, {">(s)FeOCa+"}, "O2(aq)", "e-")
std::vector< const SolidMechanicsPlasticModel * > _f
User objects that define the yield functions, flow potentials, etc.

◆ activeSurfaces()

void MultiPlasticityRawComponentAssembler::activeSurfaces ( int  model,
const std::vector< bool > &  active,
std::vector< unsigned int > &  active_surfaces 
)
protectedinherited

Returns the external surface number(s) of the active surfaces of the given model This may be of size=0 if there are no active surfaces of the given model.

Parameters
modelthe model number
activearray with entries being 'true' if the surface is active
[out]active_surfacesthe output

Definition at line 798 of file MultiPlasticityRawComponentAssembler.C.

Referenced by MultiPlasticityLinearSystem::calculateConstraints().

801 {
802  active_surfaces.resize(0);
803  for (unsigned model_surface = 0; model_surface < _f[model]->numberSurfaces(); ++model_surface)
804  if (active[_surfaces_given_model[model][model_surface]])
805  active_surfaces.push_back(_surfaces_given_model[model][model_surface]);
806 }
std::vector< std::vector< unsigned int > > _surfaces_given_model
_surfaces_given_model[model_number] = vector of surface numbers for this model
const PertinentGeochemicalSystem model(database, {"H2O", "H+", "HCO3-", "O2(aq)", "Ca++", ">(s)FeOH", "radius_neg1", "radius_neg1.5"}, {"Calcite"}, {}, {"Calcite_asdf"}, {"CH4(aq)"}, {">(s)FeOCa+"}, "O2(aq)", "e-")
std::vector< const SolidMechanicsPlasticModel * > _f
User objects that define the yield functions, flow potentials, etc.

◆ anyActiveSurfaces()

bool MultiPlasticityRawComponentAssembler::anyActiveSurfaces ( int  model,
const std::vector< bool > &  active 
)
protectedinherited

returns true if any internal surfaces of the given model are active according to 'active'

Definition at line 789 of file MultiPlasticityRawComponentAssembler.C.

Referenced by MultiPlasticityLinearSystem::calculateJacobian(), MultiPlasticityLinearSystem::calculateRHS(), checkSolution(), dof_included(), MultiPlasticityLinearSystem::nrStep(), and ComputeMultiPlasticityStress::residual2().

790 {
791  for (unsigned model_surface = 0; model_surface < _f[model]->numberSurfaces(); ++model_surface)
792  if (active[_surfaces_given_model[model][model_surface]])
793  return true;
794  return false;
795 }
std::vector< std::vector< unsigned int > > _surfaces_given_model
_surfaces_given_model[model_number] = vector of surface numbers for this model
const PertinentGeochemicalSystem model(database, {"H2O", "H+", "HCO3-", "O2(aq)", "Ca++", ">(s)FeOH", "radius_neg1", "radius_neg1.5"}, {"Calcite"}, {}, {"Calcite_asdf"}, {"CH4(aq)"}, {">(s)FeOCa+"}, "O2(aq)", "e-")
std::vector< const SolidMechanicsPlasticModel * > _f
User objects that define the yield functions, flow potentials, etc.

◆ buildActiveConstraints()

void MultiPlasticityRawComponentAssembler::buildActiveConstraints ( const std::vector< Real > &  f,
const RankTwoTensor stress,
const std::vector< Real > &  intnl,
const RankFourTensor Eijkl,
std::vector< bool > &  act 
)
protectedvirtualinherited

Constructs a set of active constraints, given the yield functions, f.

This uses SolidMechanicsPlasticModel::activeConstraints to identify the active constraints for each model.

Parameters
fyield functions (should be _num_surfaces of these)
stressstress tensor
intnlinternal parameters
Eijklelasticity tensor (stress = Eijkl*strain)
[out]actthe set of active constraints (will be resized to _num_surfaces)

Definition at line 342 of file MultiPlasticityRawComponentAssembler.C.

Referenced by ComputeMultiPlasticityStress::returnMap().

347 {
348  mooseAssert(f.size() == _num_surfaces,
349  "buildActiveConstraints called with f.size = " << f.size() << " while there are "
350  << _num_surfaces << " surfaces");
351  mooseAssert(intnl.size() == _num_models,
352  "buildActiveConstraints called with intnl.size = "
353  << intnl.size() << " while there are " << _num_models << " models");
354 
355  if (_specialIC == "rock")
356  buildActiveConstraintsRock(f, stress, intnl, Eijkl, act);
357  else if (_specialIC == "joint")
358  buildActiveConstraintsJoint(f, stress, intnl, Eijkl, act);
359  else // no specialIC
360  {
361  act.resize(0);
362  unsigned ind = 0;
363  for (unsigned model = 0; model < _num_models; ++model)
364  {
365  std::vector<Real> model_f(0);
366  for (unsigned model_surface = 0; model_surface < _f[model]->numberSurfaces(); ++model_surface)
367  model_f.push_back(f[ind++]);
368  std::vector<bool> model_act;
369  RankTwoTensor returned_stress;
370  _f[model]->activeConstraints(
371  model_f, stress, intnl[model], Eijkl, model_act, returned_stress);
372  for (unsigned model_surface = 0; model_surface < _f[model]->numberSurfaces(); ++model_surface)
373  act.push_back(model_act[model_surface]);
374  }
375  }
376 }
MooseEnum _specialIC
Allows initial set of active constraints to be chosen optimally.
void buildActiveConstraintsRock(const std::vector< Real > &f, const RankTwoTensor &stress, const std::vector< Real > &intnl, const RankFourTensor &Eijkl, std::vector< bool > &act)
"Rock" version Constructs a set of active constraints, given the yield functions, f...
void buildActiveConstraintsJoint(const std::vector< Real > &f, const RankTwoTensor &stress, const std::vector< Real > &intnl, const RankFourTensor &Eijkl, std::vector< bool > &act)
"Joint" version Constructs a set of active constraints, given the yield functions, f.
unsigned int _num_surfaces
Number of surfaces within the plastic models.
Real f(Real x)
Test function for Brents method.
const PertinentGeochemicalSystem model(database, {"H2O", "H+", "HCO3-", "O2(aq)", "Ca++", ">(s)FeOH", "radius_neg1", "radius_neg1.5"}, {"Calcite"}, {}, {"Calcite_asdf"}, {"CH4(aq)"}, {">(s)FeOCa+"}, "O2(aq)", "e-")
unsigned int _num_models
Number of plastic models for this material.
std::vector< const SolidMechanicsPlasticModel * > _f
User objects that define the yield functions, flow potentials, etc.

◆ calculateConstraints()

void MultiPlasticityLinearSystem::calculateConstraints ( const RankTwoTensor stress,
const std::vector< Real > &  intnl_old,
const std::vector< Real > &  intnl,
const std::vector< Real > &  pm,
const RankTwoTensor delta_dp,
std::vector< Real > &  f,
std::vector< RankTwoTensor > &  r,
RankTwoTensor epp,
std::vector< Real > &  ic,
const std::vector< bool > &  active 
)
protectedvirtualinherited

The constraints.

These are set to zero (or <=0 in the case of the yield functions) by the Newton-Raphson process, except in the case of linear-dependence which complicates things.

Parameters
stressThe stress
intnl_oldold values of the internal parameters
intnlinternal parameters
pmCurrent value(s) of the plasticity multiplier(s) (consistency parameters)
delta_dpChange in plastic strain incurred so far during the return
[out]fActive yield function(s)
[out]rActive flow directions
[out]eppPlastic-strain increment constraint
[out]icActive internal-parameter constraint
activeThe active constraints.

Definition at line 227 of file MultiPlasticityLinearSystem.C.

Referenced by MultiPlasticityLinearSystem::calculateRHS(), and ComputeMultiPlasticityStress::lineSearch().

237 {
238  // see comments at the start of .h file
239 
240  mooseAssert(intnl_old.size() == _num_models,
241  "Size of intnl_old is " << intnl_old.size()
242  << " which is incorrect in calculateConstraints");
243  mooseAssert(intnl.size() == _num_models,
244  "Size of intnl is " << intnl.size() << " which is incorrect in calculateConstraints");
245  mooseAssert(pm.size() == _num_surfaces,
246  "Size of pm is " << pm.size() << " which is incorrect in calculateConstraints");
247  mooseAssert(active.size() == _num_surfaces,
248  "Size of active is " << active.size()
249  << " which is incorrect in calculateConstraints");
250 
251  // yield functions
252  yieldFunction(stress, intnl, active, f);
253 
254  // flow directions and "epp"
255  flowPotential(stress, intnl, active, r);
256  epp = RankTwoTensor();
257  unsigned ind = 0;
258  for (unsigned surface = 0; surface < _num_surfaces; ++surface)
259  if (active[surface])
260  epp += pm[surface] * r[ind++]; // note, even the deactivated_due_to_ld must get added in
261  epp -= delta_dp;
262 
263  // internal constraints
264  std::vector<Real> h;
265  hardPotential(stress, intnl, active, h);
266  ic.resize(0);
267  ind = 0;
268  std::vector<unsigned int> active_surfaces;
269  std::vector<unsigned int>::iterator active_surface;
270  for (unsigned model = 0; model < _num_models; ++model)
271  {
272  activeSurfaces(model, active, active_surfaces);
273  if (active_surfaces.size() > 0)
274  {
275  // some surfaces are active in this model, so must form an internal constraint
276  ic.push_back(intnl[model] - intnl_old[model]);
277  for (active_surface = active_surfaces.begin(); active_surface != active_surfaces.end();
278  ++active_surface)
279  ic[ic.size() - 1] += pm[*active_surface] * h[ind++]; // we know the correct one is h[ind]
280  // since it was constructed in the same
281  // manner
282  }
283  }
284 }
virtual void yieldFunction(const RankTwoTensor &stress, const std::vector< Real > &intnl, const std::vector< bool > &active, std::vector< Real > &f)
The active yield function(s)
unsigned int _num_surfaces
Number of surfaces within the plastic models.
Real f(Real x)
Test function for Brents method.
virtual void flowPotential(const RankTwoTensor &stress, const std::vector< Real > &intnl, const std::vector< bool > &active, std::vector< RankTwoTensor > &r)
The active flow potential(s) - one for each yield function.
const PertinentGeochemicalSystem model(database, {"H2O", "H+", "HCO3-", "O2(aq)", "Ca++", ">(s)FeOH", "radius_neg1", "radius_neg1.5"}, {"Calcite"}, {}, {"Calcite_asdf"}, {"CH4(aq)"}, {">(s)FeOCa+"}, "O2(aq)", "e-")
unsigned int _num_models
Number of plastic models for this material.
virtual void hardPotential(const RankTwoTensor &stress, const std::vector< Real > &intnl, const std::vector< bool > &active, std::vector< Real > &h)
The active hardening potentials (one for each internal parameter and for each yield function) by assu...
void activeSurfaces(int model, const std::vector< bool > &active, std::vector< unsigned int > &active_surfaces)
Returns the external surface number(s) of the active surfaces of the given model This may be of size=...

◆ calculateJacobian()

void MultiPlasticityLinearSystem::calculateJacobian ( const RankTwoTensor stress,
const std::vector< Real > &  intnl,
const std::vector< Real > &  pm,
const RankFourTensor E_inv,
const std::vector< bool > &  active,
const std::vector< bool > &  deactivated_due_to_ld,
std::vector< std::vector< Real >> &  jac 
)
protectedvirtualinherited

d(rhs)/d(dof)

Definition at line 366 of file MultiPlasticityLinearSystem.C.

Referenced by checkJacobian(), checkSolution(), and MultiPlasticityLinearSystem::nrStep().

373 {
374  // see comments at the start of .h file
375 
376  mooseAssert(intnl.size() == _num_models,
377  "Size of intnl is " << intnl.size() << " which is incorrect in calculateJacobian");
378  mooseAssert(pm.size() == _num_surfaces,
379  "Size of pm is " << pm.size() << " which is incorrect in calculateJacobian");
380  mooseAssert(active.size() == _num_surfaces,
381  "Size of active is " << active.size() << " which is incorrect in calculateJacobian");
382  mooseAssert(deactivated_due_to_ld.size() == _num_surfaces,
383  "Size of deactivated_due_to_ld is " << deactivated_due_to_ld.size()
384  << " which is incorrect in calculateJacobian");
385 
386  unsigned ind = 0;
387  unsigned active_surface_ind = 0;
388 
389  std::vector<bool> active_surface(_num_surfaces); // active and not deactivated_due_to_ld
390  for (unsigned surface = 0; surface < _num_surfaces; ++surface)
391  active_surface[surface] = (active[surface] && !deactivated_due_to_ld[surface]);
392  unsigned num_active_surface = 0;
393  for (unsigned surface = 0; surface < _num_surfaces; ++surface)
394  if (active_surface[surface])
395  num_active_surface++;
396 
397  std::vector<bool> active_model(
398  _num_models); // whether a model has surfaces that are active and not deactivated_due_to_ld
399  for (unsigned model = 0; model < _num_models; ++model)
400  active_model[model] = anyActiveSurfaces(model, active_surface);
401 
402  unsigned num_active_model = 0;
403  for (unsigned model = 0; model < _num_models; ++model)
404  if (active_model[model])
405  num_active_model++;
406 
407  ind = 0;
408  std::vector<unsigned int> active_model_index(_num_models);
409  for (unsigned model = 0; model < _num_models; ++model)
410  if (active_model[model])
411  active_model_index[model] = ind++;
412  else
413  active_model_index[model] =
414  _num_models + 1; // just a dummy, that will probably cause a crash if something goes wrong
415 
416  std::vector<RankTwoTensor> df_dstress;
417  dyieldFunction_dstress(stress, intnl, active_surface, df_dstress);
418 
419  std::vector<Real> df_dintnl;
420  dyieldFunction_dintnl(stress, intnl, active_surface, df_dintnl);
421 
422  std::vector<RankTwoTensor> r;
423  flowPotential(stress, intnl, active, r);
424 
425  std::vector<RankFourTensor> dr_dstress;
426  dflowPotential_dstress(stress, intnl, active, dr_dstress);
427 
428  std::vector<RankTwoTensor> dr_dintnl;
429  dflowPotential_dintnl(stress, intnl, active, dr_dintnl);
430 
431  std::vector<Real> h;
432  hardPotential(stress, intnl, active, h);
433 
434  std::vector<RankTwoTensor> dh_dstress;
435  dhardPotential_dstress(stress, intnl, active, dh_dstress);
436 
437  std::vector<Real> dh_dintnl;
438  dhardPotential_dintnl(stress, intnl, active, dh_dintnl);
439 
440  // d(epp)/dstress = sum_{active alpha} pm[alpha]*dr_dstress
441  RankFourTensor depp_dstress;
442  ind = 0;
443  for (unsigned surface = 0; surface < _num_surfaces; ++surface)
444  if (active[surface]) // includes deactivated_due_to_ld
445  depp_dstress += pm[surface] * dr_dstress[ind++];
446  depp_dstress += E_inv;
447 
448  // d(epp)/dpm_{active_surface_index} = r_{active_surface_index}
449  std::vector<RankTwoTensor> depp_dpm;
450  depp_dpm.resize(num_active_surface);
451  ind = 0;
452  active_surface_ind = 0;
453  for (unsigned surface = 0; surface < _num_surfaces; ++surface)
454  {
455  if (active[surface])
456  {
457  if (active_surface[surface]) // do not include the deactived_due_to_ld, since their pm are not
458  // dofs in the NR
459  depp_dpm[active_surface_ind++] = r[ind];
460  ind++;
461  }
462  }
463 
464  // d(epp)/dintnl_{active_model_index} = sum(pm[asdf]*dr_dintnl[fdsa])
465  std::vector<RankTwoTensor> depp_dintnl;
466  depp_dintnl.assign(num_active_model, RankTwoTensor());
467  ind = 0;
468  for (unsigned surface = 0; surface < _num_surfaces; ++surface)
469  {
470  if (active[surface])
471  {
472  unsigned int model_num = modelNumber(surface);
473  if (active_model[model_num]) // only include models with surfaces which are still active after
474  // deactivated_due_to_ld
475  depp_dintnl[active_model_index[model_num]] += pm[surface] * dr_dintnl[ind];
476  ind++;
477  }
478  }
479 
480  // df_dstress has been calculated above
481  // df_dpm is always zero
482  // df_dintnl has been calculated above, but only the active_surface+active_model stuff needs to be
483  // included in Jacobian: see below
484 
485  std::vector<RankTwoTensor> dic_dstress;
486  dic_dstress.assign(num_active_model, RankTwoTensor());
487  ind = 0;
488  for (unsigned surface = 0; surface < _num_surfaces; ++surface)
489  {
490  if (active[surface])
491  {
492  unsigned int model_num = modelNumber(surface);
493  if (active_model[model_num]) // only include ic for models with active_surface (ie, if model
494  // only contains deactivated_due_to_ld don't include it)
495  dic_dstress[active_model_index[model_num]] += pm[surface] * dh_dstress[ind];
496  ind++;
497  }
498  }
499 
500  std::vector<std::vector<Real>> dic_dpm;
501  dic_dpm.resize(num_active_model);
502  ind = 0;
503  active_surface_ind = 0;
504  for (unsigned model = 0; model < num_active_model; ++model)
505  dic_dpm[model].assign(num_active_surface, 0);
506  for (unsigned surface = 0; surface < _num_surfaces; ++surface)
507  {
508  if (active[surface])
509  {
510  if (active_surface[surface]) // only take derivs wrt active-but-not-deactivated_due_to_ld pm
511  {
512  unsigned int model_num = modelNumber(surface);
513  // if (active_model[model_num]) // do not need this check as if the surface has
514  // active_surface, the model must be deemed active!
515  dic_dpm[active_model_index[model_num]][active_surface_ind] = h[ind];
516  active_surface_ind++;
517  }
518  ind++;
519  }
520  }
521 
522  std::vector<std::vector<Real>> dic_dintnl;
523  dic_dintnl.resize(num_active_model);
524  for (unsigned model = 0; model < num_active_model; ++model)
525  {
526  dic_dintnl[model].assign(num_active_model, 0);
527  dic_dintnl[model][model] = 1; // deriv wrt internal parameter
528  }
529  ind = 0;
530  for (unsigned surface = 0; surface < _num_surfaces; ++surface)
531  {
532  if (active[surface])
533  {
534  unsigned int model_num = modelNumber(surface);
535  if (active_model[model_num]) // only the models that contain surfaces that are still active
536  // after deactivation_due_to_ld
537  dic_dintnl[active_model_index[model_num]][active_model_index[model_num]] +=
538  pm[surface] * dh_dintnl[ind];
539  ind++;
540  }
541  }
542 
543  unsigned int dim = 3;
544  unsigned int system_size =
545  6 + num_active_surface + num_active_model; // "6" comes from symmeterizing epp
546  jac.resize(system_size);
547  for (unsigned i = 0; i < system_size; ++i)
548  jac[i].assign(system_size, 0);
549 
550  unsigned int row_num = 0;
551  unsigned int col_num = 0;
552  for (unsigned i = 0; i < dim; ++i)
553  for (unsigned j = 0; j <= i; ++j)
554  {
555  for (unsigned k = 0; k < dim; ++k)
556  for (unsigned l = 0; l <= k; ++l)
557  jac[col_num][row_num++] =
558  depp_dstress(i, j, k, l) +
559  (k != l ? depp_dstress(i, j, l, k)
560  : 0); // extra part is needed because i assume dstress(i, j) = dstress(j, i)
561  for (unsigned surface = 0; surface < num_active_surface; ++surface)
562  jac[col_num][row_num++] = depp_dpm[surface](i, j);
563  for (unsigned a = 0; a < num_active_model; ++a)
564  jac[col_num][row_num++] = depp_dintnl[a](i, j);
565  row_num = 0;
566  col_num++;
567  }
568 
569  ind = 0;
570  for (unsigned surface = 0; surface < _num_surfaces; ++surface)
571  if (active_surface[surface])
572  {
573  for (unsigned k = 0; k < dim; ++k)
574  for (unsigned l = 0; l <= k; ++l)
575  jac[col_num][row_num++] =
576  df_dstress[ind](k, l) +
577  (k != l ? df_dstress[ind](l, k)
578  : 0); // extra part is needed because i assume dstress(i, j) = dstress(j, i)
579  for (unsigned beta = 0; beta < num_active_surface; ++beta)
580  jac[col_num][row_num++] = 0; // df_dpm
581  for (unsigned model = 0; model < _num_models; ++model)
582  if (active_model[model]) // only use df_dintnl for models in active_model
583  {
584  if (modelNumber(surface) == model)
585  jac[col_num][row_num++] = df_dintnl[ind];
586  else
587  jac[col_num][row_num++] = 0;
588  }
589  ind++;
590  row_num = 0;
591  col_num++;
592  }
593 
594  for (unsigned a = 0; a < num_active_model; ++a)
595  {
596  for (unsigned k = 0; k < dim; ++k)
597  for (unsigned l = 0; l <= k; ++l)
598  jac[col_num][row_num++] =
599  dic_dstress[a](k, l) +
600  (k != l ? dic_dstress[a](l, k)
601  : 0); // extra part is needed because i assume dstress(i, j) = dstress(j, i)
602  for (unsigned alpha = 0; alpha < num_active_surface; ++alpha)
603  jac[col_num][row_num++] = dic_dpm[a][alpha];
604  for (unsigned b = 0; b < num_active_model; ++b)
605  jac[col_num][row_num++] = dic_dintnl[a][b];
606  row_num = 0;
607  col_num++;
608  }
609 
610  mooseAssert(col_num == system_size, "Incorrect filling of cols in Jacobian");
611 }
virtual void dyieldFunction_dstress(const RankTwoTensor &stress, const std::vector< Real > &intnl, const std::vector< bool > &active, std::vector< RankTwoTensor > &df_dstress)
The derivative of the active yield function(s) with respect to stress.
unsigned int modelNumber(unsigned int surface)
returns the model number, given the surface number
unsigned int dim
virtual void dhardPotential_dintnl(const RankTwoTensor &stress, const std::vector< Real > &intnl, const std::vector< bool > &active, std::vector< Real > &dh_dintnl)
The derivative of the active hardening potentials with respect to the active internal parameters...
virtual void dhardPotential_dstress(const RankTwoTensor &stress, const std::vector< Real > &intnl, const std::vector< bool > &active, std::vector< RankTwoTensor > &dh_dstress)
The derivative of the active hardening potentials with respect to stress By assumption in the Userobj...
unsigned int _num_surfaces
Number of surfaces within the plastic models.
virtual void dflowPotential_dintnl(const RankTwoTensor &stress, const std::vector< Real > &intnl, const std::vector< bool > &active, std::vector< RankTwoTensor > &dr_dintnl)
The derivative of the active flow potentials with respect to the active internal parameters The UserO...
bool anyActiveSurfaces(int model, const std::vector< bool > &active)
returns true if any internal surfaces of the given model are active according to &#39;active&#39; ...
virtual void flowPotential(const RankTwoTensor &stress, const std::vector< Real > &intnl, const std::vector< bool > &active, std::vector< RankTwoTensor > &r)
The active flow potential(s) - one for each yield function.
const PertinentGeochemicalSystem model(database, {"H2O", "H+", "HCO3-", "O2(aq)", "Ca++", ">(s)FeOH", "radius_neg1", "radius_neg1.5"}, {"Calcite"}, {}, {"Calcite_asdf"}, {"CH4(aq)"}, {">(s)FeOCa+"}, "O2(aq)", "e-")
unsigned int _num_models
Number of plastic models for this material.
static const std::string alpha
Definition: NS.h:126
virtual void dflowPotential_dstress(const RankTwoTensor &stress, const std::vector< Real > &intnl, const std::vector< bool > &active, std::vector< RankFourTensor > &dr_dstress)
The derivative of the active flow potential(s) with respect to stress.
virtual void hardPotential(const RankTwoTensor &stress, const std::vector< Real > &intnl, const std::vector< bool > &active, std::vector< Real > &h)
The active hardening potentials (one for each internal parameter and for each yield function) by assu...
static const std::complex< double > j(0, 1)
Complex number "j" (also known as "i")
virtual void dyieldFunction_dintnl(const RankTwoTensor &stress, const std::vector< Real > &intnl, const std::vector< bool > &active, std::vector< Real > &df_dintnl)
The derivative of active yield function(s) with respect to their internal parameters (the user object...
static const std::string k
Definition: NS.h:124

◆ calculateRHS()

void MultiPlasticityLinearSystem::calculateRHS ( const RankTwoTensor stress,
const std::vector< Real > &  intnl_old,
const std::vector< Real > &  intnl,
const std::vector< Real > &  pm,
const RankTwoTensor delta_dp,
std::vector< Real > &  rhs,
const std::vector< bool > &  active,
bool  eliminate_ld,
std::vector< bool > &  deactivated_due_to_ld 
)
protectedvirtualinherited

Calculate the RHS which is rhs = -(epp(0,0), epp(1,0), epp(1,1), epp(2,0), epp(2,1), epp(2,2), f[0], f[1], ..., f[num_f], ic[0], ic[1], ..., ic[num_ic])

Note that the 'epp' components only contain the upper diagonal. These contain flow directions and plasticity-multipliers for all active surfaces, even the deactivated_due_to_ld surfaces. Note that the 'f' components only contain the active and not deactivated_due_to_ld surfaces Note that the 'ic' components only contain the internal constraints for models which contain active and not deactivated_due_to_ld surfaces. They contain hardening-potentials and plasticity-multipliers for the active surfaces, even the deactivated_due_to_ld surfaces

Parameters
stressThe stress
intnl_oldold values of the internal parameters
intnlinternal parameters
pmCurrent value(s) of the plasticity multiplier(s) (consistency parameters)
delta_dpChange in plastic strain incurred so far during the return
[out]rhsthe rhs
activeThe active constraints.
eliminate_ldCheck for linear dependence of constraints and put the results into deactivated_due_to_ld. Usually this should be true, but for certain debug operations it should be false
[out]deactivated_due_to_ldconstraints deactivated due to linear-dependence of flow directions

Definition at line 287 of file MultiPlasticityLinearSystem.C.

Referenced by checkSolution(), fdJacobian(), and MultiPlasticityLinearSystem::nrStep().

296 {
297  // see comments at the start of .h file
298 
299  mooseAssert(intnl_old.size() == _num_models,
300  "Size of intnl_old is " << intnl_old.size() << " which is incorrect in calculateRHS");
301  mooseAssert(intnl.size() == _num_models,
302  "Size of intnl is " << intnl.size() << " which is incorrect in calculateRHS");
303  mooseAssert(pm.size() == _num_surfaces,
304  "Size of pm is " << pm.size() << " which is incorrect in calculateRHS");
305  mooseAssert(active.size() == _num_surfaces,
306  "Size of active is " << active.size() << " which is incorrect in calculateRHS");
307 
308  std::vector<Real> f; // the yield functions
309  RankTwoTensor epp; // the plastic-strain constraint ("direction constraint")
310  std::vector<Real> ic; // the "internal constraints"
311 
312  std::vector<RankTwoTensor> r;
313  calculateConstraints(stress, intnl_old, intnl, pm, delta_dp, f, r, epp, ic, active);
314 
315  if (eliminate_ld)
316  eliminateLinearDependence(stress, intnl, f, r, active, deactivated_due_to_ld);
317  else
318  deactivated_due_to_ld.assign(_num_surfaces, false);
319 
320  std::vector<bool> active_not_deact(_num_surfaces);
321  for (unsigned surface = 0; surface < _num_surfaces; ++surface)
322  active_not_deact[surface] = (active[surface] && !deactivated_due_to_ld[surface]);
323 
324  unsigned num_active_f = 0;
325  for (unsigned surface = 0; surface < _num_surfaces; ++surface)
326  if (active_not_deact[surface])
327  num_active_f++;
328 
329  unsigned num_active_ic = 0;
330  for (unsigned model = 0; model < _num_models; ++model)
331  if (anyActiveSurfaces(model, active_not_deact))
332  num_active_ic++;
333 
334  unsigned int dim = 3;
335  unsigned int system_size = 6 + num_active_f + num_active_ic; // "6" comes from symmeterizing epp,
336  // num_active_f comes from "f",
337  // num_active_f comes from "ic"
338 
339  rhs.resize(system_size);
340 
341  unsigned ind = 0;
342  for (unsigned i = 0; i < dim; ++i)
343  for (unsigned j = 0; j <= i; ++j)
344  rhs[ind++] = -epp(i, j);
345  unsigned active_surface = 0;
346  for (unsigned surface = 0; surface < _num_surfaces; ++surface)
347  if (active[surface])
348  {
349  if (!deactivated_due_to_ld[surface])
350  rhs[ind++] = -f[active_surface];
351  active_surface++;
352  }
353  unsigned active_model = 0;
354  for (unsigned model = 0; model < _num_models; ++model)
355  if (anyActiveSurfaces(model, active))
356  {
357  if (anyActiveSurfaces(model, active_not_deact))
358  rhs[ind++] = -ic[active_model];
359  active_model++;
360  }
361 
362  mooseAssert(ind == system_size, "Incorrect filling of the rhs in calculateRHS");
363 }
unsigned int dim
virtual void eliminateLinearDependence(const RankTwoTensor &stress, const std::vector< Real > &intnl, const std::vector< Real > &f, const std::vector< RankTwoTensor > &r, const std::vector< bool > &active, std::vector< bool > &deactivated_due_to_ld)
Performs a number of singular-value decompositions to check for linear-dependence of the active direc...
unsigned int _num_surfaces
Number of surfaces within the plastic models.
Real f(Real x)
Test function for Brents method.
bool anyActiveSurfaces(int model, const std::vector< bool > &active)
returns true if any internal surfaces of the given model are active according to &#39;active&#39; ...
const PertinentGeochemicalSystem model(database, {"H2O", "H+", "HCO3-", "O2(aq)", "Ca++", ">(s)FeOH", "radius_neg1", "radius_neg1.5"}, {"Calcite"}, {}, {"Calcite_asdf"}, {"CH4(aq)"}, {">(s)FeOCa+"}, "O2(aq)", "e-")
unsigned int _num_models
Number of plastic models for this material.
static const std::complex< double > j(0, 1)
Complex number "j" (also known as "i")
virtual void calculateConstraints(const RankTwoTensor &stress, const std::vector< Real > &intnl_old, const std::vector< Real > &intnl, const std::vector< Real > &pm, const RankTwoTensor &delta_dp, std::vector< Real > &f, std::vector< RankTwoTensor > &r, RankTwoTensor &epp, std::vector< Real > &ic, const std::vector< bool > &active)
The constraints.

◆ checkDerivatives()

void MultiPlasticityDebugger::checkDerivatives ( )

Checks the derivatives, eg dyieldFunction_dstress by using finite difference approximations.

Definition at line 85 of file MultiPlasticityDebugger.C.

Referenced by ComputeMultiPlasticityStress::initQpStatefulProperties(), and ComputeMultiPlasticityStress::plasticStep().

86 {
87  Moose::err
88  << "\n\n++++++++++++++++++++++++\nChecking the derivatives\n++++++++++++++++++++++++\n";
90 
91  std::vector<bool> act;
92  act.assign(_num_surfaces, true);
93 
94  Moose::err << "\ndyieldFunction_dstress. Relative L2 norms.\n";
95  std::vector<RankTwoTensor> df_dstress;
96  std::vector<RankTwoTensor> fddf_dstress;
99  for (unsigned surface = 0; surface < _num_surfaces; ++surface)
100  {
101  Moose::err << "surface = " << surface << " Relative L2norm = "
102  << 2 * (df_dstress[surface] - fddf_dstress[surface]).L2norm() /
103  (df_dstress[surface] + fddf_dstress[surface]).L2norm()
104  << "\n";
105  Moose::err << "Coded:\n";
106  df_dstress[surface].print();
107  Moose::err << "Finite difference:\n";
108  fddf_dstress[surface].print();
109  }
110 
111  Moose::err << "\ndyieldFunction_dintnl.\n";
112  std::vector<Real> df_dintnl;
114  Moose::err << "Coded:\n";
115  for (unsigned surface = 0; surface < _num_surfaces; ++surface)
116  Moose::err << df_dintnl[surface] << " ";
117  Moose::err << "\n";
118  std::vector<Real> fddf_dintnl;
120  Moose::err << "Finite difference:\n";
121  for (unsigned surface = 0; surface < _num_surfaces; ++surface)
122  Moose::err << fddf_dintnl[surface] << " ";
123  Moose::err << "\n";
124 
125  Moose::err << "\ndflowPotential_dstress. Relative L2 norms.\n";
126  std::vector<RankFourTensor> dr_dstress;
127  std::vector<RankFourTensor> fddr_dstress;
130  for (unsigned surface = 0; surface < _num_surfaces; ++surface)
131  {
132  Moose::err << "surface = " << surface << " Relative L2norm = "
133  << 2 * (dr_dstress[surface] - fddr_dstress[surface]).L2norm() /
134  (dr_dstress[surface] + fddr_dstress[surface]).L2norm()
135  << "\n";
136  Moose::err << "Coded:\n";
137  dr_dstress[surface].print();
138  Moose::err << "Finite difference:\n";
139  fddr_dstress[surface].print();
140  }
141 
142  Moose::err << "\ndflowPotential_dintnl. Relative L2 norms.\n";
143  std::vector<RankTwoTensor> dr_dintnl;
144  std::vector<RankTwoTensor> fddr_dintnl;
147  for (unsigned surface = 0; surface < _num_surfaces; ++surface)
148  {
149  Moose::err << "surface = " << surface << " Relative L2norm = "
150  << 2 * (dr_dintnl[surface] - fddr_dintnl[surface]).L2norm() /
151  (dr_dintnl[surface] + fddr_dintnl[surface]).L2norm()
152  << "\n";
153  Moose::err << "Coded:\n";
154  dr_dintnl[surface].print();
155  Moose::err << "Finite difference:\n";
156  fddr_dintnl[surface].print();
157  }
158 
159  Moose::err << std::flush;
160 }
void fddyieldFunction_dintnl(const RankTwoTensor &stress, const std::vector< Real > &intnl, std::vector< Real > &df_dintnl)
The finite-difference derivative of yield function(s) with respect to internal parameter(s) ...
void fddyieldFunction_dstress(const RankTwoTensor &stress, const std::vector< Real > &intnl, std::vector< RankTwoTensor > &df_dstress)
The finite-difference derivative of yield function(s) with respect to stress.
virtual void dyieldFunction_dstress(const RankTwoTensor &stress, const std::vector< Real > &intnl, const std::vector< bool > &active, std::vector< RankTwoTensor > &df_dstress)
The derivative of the active yield function(s) with respect to stress.
T L2norm(const RankTwoTensorTempl< T > &r2tensor)
virtual void fddflowPotential_dintnl(const RankTwoTensor &stress, const std::vector< Real > &intnl, std::vector< RankTwoTensor > &dr_dintnl)
The finite-difference derivative of the flow potentials with respect to internal parameters.
virtual void fddflowPotential_dstress(const RankTwoTensor &stress, const std::vector< Real > &intnl, std::vector< RankFourTensor > &dr_dstress)
The finite-difference derivative of the flow potential(s) with respect to stress. ...
std::vector< Real > _fspb_debug_intnl
Debug the Jacobian entires at these internal parameters.
unsigned int _num_surfaces
Number of surfaces within the plastic models.
void outputAndCheckDebugParameters()
Outputs the debug parameters: _fspb_debug_stress, _fspd_debug_pm, etc and checks that they are sized ...
virtual void dflowPotential_dintnl(const RankTwoTensor &stress, const std::vector< Real > &intnl, const std::vector< bool > &active, std::vector< RankTwoTensor > &dr_dintnl)
The derivative of the active flow potentials with respect to the active internal parameters The UserO...
RankTwoTensor _fspb_debug_stress
Debug the Jacobian entries at this stress.
virtual void dflowPotential_dstress(const RankTwoTensor &stress, const std::vector< Real > &intnl, const std::vector< bool > &active, std::vector< RankFourTensor > &dr_dstress)
The derivative of the active flow potential(s) with respect to stress.
virtual void dyieldFunction_dintnl(const RankTwoTensor &stress, const std::vector< Real > &intnl, const std::vector< bool > &active, std::vector< Real > &df_dintnl)
The derivative of active yield function(s) with respect to their internal parameters (the user object...

◆ checkJacobian()

void MultiPlasticityDebugger::checkJacobian ( const RankFourTensor E_inv,
const std::vector< Real > &  intnl_old 
)

Checks the full Jacobian, which is just certain linear combinations of the dyieldFunction_dstress, etc, by using finite difference approximations.

Definition at line 163 of file MultiPlasticityDebugger.C.

Referenced by ComputeMultiPlasticityStress::computeQpStress(), and ComputeMultiPlasticityStress::plasticStep().

165 {
166  Moose::err << "\n\n+++++++++++++++++++++\nChecking the Jacobian\n+++++++++++++++++++++\n";
168 
169  std::vector<bool> act;
170  act.assign(_num_surfaces, true);
171  std::vector<bool> deactivated_due_to_ld;
172  deactivated_due_to_ld.assign(_num_surfaces, false);
173 
174  RankTwoTensor delta_dp = -E_inv * _fspb_debug_stress;
175 
176  std::vector<std::vector<Real>> jac;
180  E_inv,
181  act,
182  deactivated_due_to_ld,
183  jac);
184 
185  std::vector<std::vector<Real>> fdjac;
187  intnl_old,
190  delta_dp,
191  E_inv,
192  false,
193  fdjac);
194 
195  Real L2_numer = 0;
196  Real L2_denom = 0;
197  for (unsigned row = 0; row < jac.size(); ++row)
198  for (unsigned col = 0; col < jac.size(); ++col)
199  {
200  L2_numer += Utility::pow<2>(jac[row][col] - fdjac[row][col]);
201  L2_denom += Utility::pow<2>(jac[row][col] + fdjac[row][col]);
202  }
203  Moose::err << "\nRelative L2norm = " << std::sqrt(L2_numer / L2_denom) / 0.5 << "\n";
204 
205  Moose::err << "\nHand-coded Jacobian:\n";
206  for (unsigned row = 0; row < jac.size(); ++row)
207  {
208  for (unsigned col = 0; col < jac.size(); ++col)
209  Moose::err << jac[row][col] << " ";
210  Moose::err << "\n";
211  }
212 
213  Moose::err << "Finite difference Jacobian:\n";
214  for (unsigned row = 0; row < fdjac.size(); ++row)
215  {
216  for (unsigned col = 0; col < fdjac.size(); ++col)
217  Moose::err << fdjac[row][col] << " ";
218  Moose::err << "\n";
219  }
220 
221  Moose::err << std::flush;
222 }
virtual void fdJacobian(const RankTwoTensor &stress, const std::vector< Real > &intnl_old, const std::vector< Real > &intnl, const std::vector< Real > &pm, const RankTwoTensor &delta_dp, const RankFourTensor &E_inv, bool eliminate_ld, std::vector< std::vector< Real >> &jac)
The Jacobian calculated using finite differences.
virtual void calculateJacobian(const RankTwoTensor &stress, const std::vector< Real > &intnl, const std::vector< Real > &pm, const RankFourTensor &E_inv, const std::vector< bool > &active, const std::vector< bool > &deactivated_due_to_ld, std::vector< std::vector< Real >> &jac)
d(rhs)/d(dof)
ADRealEigenVector< T, D, asd > sqrt(const ADRealEigenVector< T, D, asd > &)
std::vector< Real > _fspb_debug_intnl
Debug the Jacobian entires at these internal parameters.
unsigned int _num_surfaces
Number of surfaces within the plastic models.
void outputAndCheckDebugParameters()
Outputs the debug parameters: _fspb_debug_stress, _fspd_debug_pm, etc and checks that they are sized ...
RankTwoTensor _fspb_debug_stress
Debug the Jacobian entries at this stress.
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
std::vector< Real > _fspb_debug_pm
Debug the Jacobian entires at these plastic multipliers.

◆ checkSolution()

void MultiPlasticityDebugger::checkSolution ( const RankFourTensor E_inv)

Checks that Ax does equal b in the NR procedure.

Definition at line 371 of file MultiPlasticityDebugger.C.

Referenced by ComputeMultiPlasticityStress::computeQpStress(), and ComputeMultiPlasticityStress::plasticStep().

372 {
373  Moose::err << "\n\n+++++++++++++++++++++\nChecking the Solution\n";
374  Moose::err << "(Ie, checking Ax = b)\n+++++++++++++++++++++\n";
376 
377  std::vector<bool> act;
378  act.assign(_num_surfaces, true);
379  std::vector<bool> deactivated_due_to_ld;
380  deactivated_due_to_ld.assign(_num_surfaces, false);
381 
382  RankTwoTensor delta_dp = -E_inv * _fspb_debug_stress;
383 
384  std::vector<Real> orig_rhs;
389  delta_dp,
390  orig_rhs,
391  act,
392  true,
393  deactivated_due_to_ld);
394 
395  Moose::err << "\nb = ";
396  for (unsigned i = 0; i < orig_rhs.size(); ++i)
397  Moose::err << orig_rhs[i] << " ";
398  Moose::err << "\n\n";
399 
400  std::vector<std::vector<Real>> jac_coded;
404  E_inv,
405  act,
406  deactivated_due_to_ld,
407  jac_coded);
408 
409  Moose::err
410  << "Before checking Ax=b is correct, check that the Jacobians given below are equal.\n";
411  Moose::err
412  << "The hand-coded Jacobian is used in calculating the solution 'x', given 'b' above.\n";
413  Moose::err << "Note that this only includes degrees of freedom that aren't deactivated due to "
414  "linear dependence.\n";
415  Moose::err << "Hand-coded Jacobian:\n";
416  for (unsigned row = 0; row < jac_coded.size(); ++row)
417  {
418  for (unsigned col = 0; col < jac_coded.size(); ++col)
419  Moose::err << jac_coded[row][col] << " ";
420  Moose::err << "\n";
421  }
422 
423  deactivated_due_to_ld.assign(_num_surfaces,
424  false); // this potentially gets changed by nrStep, below
425  RankTwoTensor dstress;
426  std::vector<Real> dpm;
427  std::vector<Real> dintnl;
432  E_inv,
433  delta_dp,
434  dstress,
435  dpm,
436  dintnl,
437  act,
438  deactivated_due_to_ld);
439 
440  std::vector<bool> active_not_deact(_num_surfaces);
441  for (unsigned surface = 0; surface < _num_surfaces; ++surface)
442  active_not_deact[surface] = !deactivated_due_to_ld[surface];
443 
444  std::vector<Real> x;
445  x.assign(orig_rhs.size(), 0);
446  unsigned ind = 0;
447  for (unsigned i = 0; i < 3; ++i)
448  for (unsigned j = 0; j <= i; ++j)
449  x[ind++] = dstress(i, j);
450  for (unsigned surface = 0; surface < _num_surfaces; ++surface)
451  if (active_not_deact[surface])
452  x[ind++] = dpm[surface];
453  for (unsigned model = 0; model < _num_models; ++model)
454  if (anyActiveSurfaces(model, active_not_deact))
455  x[ind++] = dintnl[model];
456 
457  mooseAssert(ind == orig_rhs.size(),
458  "Incorrect extracting of changes from NR solution in the "
459  "finite-difference checking of nrStep");
460 
461  Moose::err << "\nThis yields x =";
462  for (unsigned i = 0; i < orig_rhs.size(); ++i)
463  Moose::err << x[i] << " ";
464  Moose::err << "\n";
465 
466  std::vector<std::vector<Real>> jac_fd;
471  delta_dp,
472  E_inv,
473  true,
474  jac_fd);
475 
476  Moose::err << "\nThe finite-difference Jacobian is used to multiply by this 'x',\n";
477  Moose::err << "in order to check that the solution is correct\n";
478  Moose::err << "Finite-difference Jacobian:\n";
479  for (unsigned row = 0; row < jac_fd.size(); ++row)
480  {
481  for (unsigned col = 0; col < jac_fd.size(); ++col)
482  Moose::err << jac_fd[row][col] << " ";
483  Moose::err << "\n";
484  }
485 
486  Real L2_numer = 0;
487  Real L2_denom = 0;
488  for (unsigned row = 0; row < jac_coded.size(); ++row)
489  for (unsigned col = 0; col < jac_coded.size(); ++col)
490  {
491  L2_numer += Utility::pow<2>(jac_coded[row][col] - jac_fd[row][col]);
492  L2_denom += Utility::pow<2>(jac_coded[row][col] + jac_fd[row][col]);
493  }
494  Moose::err << "Relative L2norm of the hand-coded and finite-difference Jacobian is "
495  << std::sqrt(L2_numer / L2_denom) / 0.5 << "\n";
496 
497  std::vector<Real> fd_times_x;
498  fd_times_x.assign(orig_rhs.size(), 0);
499  for (unsigned row = 0; row < orig_rhs.size(); ++row)
500  for (unsigned col = 0; col < orig_rhs.size(); ++col)
501  fd_times_x[row] += jac_fd[row][col] * x[col];
502 
503  Moose::err << "\n(Finite-difference Jacobian)*x =\n";
504  for (unsigned i = 0; i < orig_rhs.size(); ++i)
505  Moose::err << fd_times_x[i] << " ";
506  Moose::err << "\n";
507  Moose::err << "Recall that b = \n";
508  for (unsigned i = 0; i < orig_rhs.size(); ++i)
509  Moose::err << orig_rhs[i] << " ";
510  Moose::err << "\n";
511 
512  L2_numer = 0;
513  L2_denom = 0;
514  for (unsigned i = 0; i < orig_rhs.size(); ++i)
515  {
516  L2_numer += Utility::pow<2>(orig_rhs[i] - fd_times_x[i]);
517  L2_denom += Utility::pow<2>(orig_rhs[i] + fd_times_x[i]);
518  }
519  Moose::err << "\nRelative L2norm of these is " << std::sqrt(L2_numer / L2_denom) / 0.5
520  << std::endl;
521 }
virtual void fdJacobian(const RankTwoTensor &stress, const std::vector< Real > &intnl_old, const std::vector< Real > &intnl, const std::vector< Real > &pm, const RankTwoTensor &delta_dp, const RankFourTensor &E_inv, bool eliminate_ld, std::vector< std::vector< Real >> &jac)
The Jacobian calculated using finite differences.
virtual void calculateRHS(const RankTwoTensor &stress, const std::vector< Real > &intnl_old, const std::vector< Real > &intnl, const std::vector< Real > &pm, const RankTwoTensor &delta_dp, std::vector< Real > &rhs, const std::vector< bool > &active, bool eliminate_ld, std::vector< bool > &deactivated_due_to_ld)
Calculate the RHS which is rhs = -(epp(0,0), epp(1,0), epp(1,1), epp(2,0), epp(2,1), epp(2,2), f[0], f[1], ..., f[num_f], ic[0], ic[1], ..., ic[num_ic])
virtual void calculateJacobian(const RankTwoTensor &stress, const std::vector< Real > &intnl, const std::vector< Real > &pm, const RankFourTensor &E_inv, const std::vector< bool > &active, const std::vector< bool > &deactivated_due_to_ld, std::vector< std::vector< Real >> &jac)
d(rhs)/d(dof)
ADRealEigenVector< T, D, asd > sqrt(const ADRealEigenVector< T, D, asd > &)
std::vector< Real > _fspb_debug_intnl
Debug the Jacobian entires at these internal parameters.
unsigned int _num_surfaces
Number of surfaces within the plastic models.
const std::vector< double > x
void outputAndCheckDebugParameters()
Outputs the debug parameters: _fspb_debug_stress, _fspd_debug_pm, etc and checks that they are sized ...
RankTwoTensor _fspb_debug_stress
Debug the Jacobian entries at this stress.
bool anyActiveSurfaces(int model, const std::vector< bool > &active)
returns true if any internal surfaces of the given model are active according to &#39;active&#39; ...
const PertinentGeochemicalSystem model(database, {"H2O", "H+", "HCO3-", "O2(aq)", "Ca++", ">(s)FeOH", "radius_neg1", "radius_neg1.5"}, {"Calcite"}, {}, {"Calcite_asdf"}, {"CH4(aq)"}, {">(s)FeOCa+"}, "O2(aq)", "e-")
unsigned int _num_models
Number of plastic models for this material.
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
std::vector< Real > _fspb_debug_pm
Debug the Jacobian entires at these plastic multipliers.
static const std::complex< double > j(0, 1)
Complex number "j" (also known as "i")
virtual void nrStep(const RankTwoTensor &stress, const std::vector< Real > &intnl_old, const std::vector< Real > &intnl, const std::vector< Real > &pm, const RankFourTensor &E_inv, const RankTwoTensor &delta_dp, RankTwoTensor &dstress, std::vector< Real > &dpm, std::vector< Real > &dintnl, const std::vector< bool > &active, std::vector< bool > &deactivated_due_to_ld)
Performs one Newton-Raphson step.

◆ dflowPotential_dintnl()

void MultiPlasticityRawComponentAssembler::dflowPotential_dintnl ( const RankTwoTensor stress,
const std::vector< Real > &  intnl,
const std::vector< bool > &  active,
std::vector< RankTwoTensor > &  dr_dintnl 
)
protectedvirtualinherited

The derivative of the active flow potentials with respect to the active internal parameters The UserObjects explicitly assume that r[alpha] is only dependent on intnl[alpha].

Parameters
stressthe stress at which to calculate the flow potential
intnlvector of internal parameters
activeset of active constraints - only the active derivatives are put into "dr_dintnl"
[out]dr_dintnlthe derivatives. dr_dintnl[alpha](i, j) = dr[alpha](i, j)/dintnl[alpha]

Definition at line 233 of file MultiPlasticityRawComponentAssembler.C.

Referenced by MultiPlasticityLinearSystem::calculateJacobian(), checkDerivatives(), and ComputeMultiPlasticityStress::consistentTangentOperator().

237 {
238  mooseAssert(intnl.size() == _num_models, "Incorrect size of internal parameters");
239  mooseAssert(active.size() == _num_surfaces, "Incorrect size of active");
240 
241  dr_dintnl.resize(0);
242  std::vector<unsigned int> active_surfaces_of_model;
243  std::vector<unsigned int>::iterator active_surface;
244  std::vector<RankTwoTensor> model_dr_dintnl;
245  for (unsigned model = 0; model < _num_models; ++model)
246  {
247  activeModelSurfaces(model, active, active_surfaces_of_model);
248  if (active_surfaces_of_model.size() > 0)
249  {
250  _f[model]->dflowPotential_dintnlV(stress, intnl[model], model_dr_dintnl);
251  for (active_surface = active_surfaces_of_model.begin();
252  active_surface != active_surfaces_of_model.end();
253  ++active_surface)
254  dr_dintnl.push_back(model_dr_dintnl[*active_surface]);
255  }
256  }
257 }
void activeModelSurfaces(int model, const std::vector< bool > &active, std::vector< unsigned int > &active_surfaces_of_model)
Returns the internal surface number(s) of the active surfaces of the given model This may be of size=...
unsigned int _num_surfaces
Number of surfaces within the plastic models.
const PertinentGeochemicalSystem model(database, {"H2O", "H+", "HCO3-", "O2(aq)", "Ca++", ">(s)FeOH", "radius_neg1", "radius_neg1.5"}, {"Calcite"}, {}, {"Calcite_asdf"}, {"CH4(aq)"}, {">(s)FeOCa+"}, "O2(aq)", "e-")
unsigned int _num_models
Number of plastic models for this material.
std::vector< const SolidMechanicsPlasticModel * > _f
User objects that define the yield functions, flow potentials, etc.

◆ dflowPotential_dstress()

void MultiPlasticityRawComponentAssembler::dflowPotential_dstress ( const RankTwoTensor stress,
const std::vector< Real > &  intnl,
const std::vector< bool > &  active,
std::vector< RankFourTensor > &  dr_dstress 
)
protectedvirtualinherited

The derivative of the active flow potential(s) with respect to stress.

Parameters
stressthe stress at which to calculate the flow potential
intnlvector of internal parameters
activeset of active constraints - only the active derivatives are put into "dr_dstress"
[out]dr_dstressthe derivative. dr_dstress[alpha](i, j, k, l) = dr[alpha](i, j)/dstress(k, l)

Definition at line 205 of file MultiPlasticityRawComponentAssembler.C.

Referenced by MultiPlasticityLinearSystem::calculateJacobian(), checkDerivatives(), and ComputeMultiPlasticityStress::consistentTangentOperator().

210 {
211  mooseAssert(intnl.size() == _num_models, "Incorrect size of internal parameters");
212  mooseAssert(active.size() == _num_surfaces, "Incorrect size of active");
213 
214  dr_dstress.resize(0);
215  std::vector<unsigned int> active_surfaces_of_model;
216  std::vector<unsigned int>::iterator active_surface;
217  std::vector<RankFourTensor> model_dr_dstress;
218  for (unsigned model = 0; model < _num_models; ++model)
219  {
220  activeModelSurfaces(model, active, active_surfaces_of_model);
221  if (active_surfaces_of_model.size() > 0)
222  {
223  _f[model]->dflowPotential_dstressV(stress, intnl[model], model_dr_dstress);
224  for (active_surface = active_surfaces_of_model.begin();
225  active_surface != active_surfaces_of_model.end();
226  ++active_surface)
227  dr_dstress.push_back(model_dr_dstress[*active_surface]);
228  }
229  }
230 }
void activeModelSurfaces(int model, const std::vector< bool > &active, std::vector< unsigned int > &active_surfaces_of_model)
Returns the internal surface number(s) of the active surfaces of the given model This may be of size=...
unsigned int _num_surfaces
Number of surfaces within the plastic models.
const PertinentGeochemicalSystem model(database, {"H2O", "H+", "HCO3-", "O2(aq)", "Ca++", ">(s)FeOH", "radius_neg1", "radius_neg1.5"}, {"Calcite"}, {}, {"Calcite_asdf"}, {"CH4(aq)"}, {">(s)FeOCa+"}, "O2(aq)", "e-")
unsigned int _num_models
Number of plastic models for this material.
std::vector< const SolidMechanicsPlasticModel * > _f
User objects that define the yield functions, flow potentials, etc.

◆ dhardPotential_dintnl()

void MultiPlasticityRawComponentAssembler::dhardPotential_dintnl ( const RankTwoTensor stress,
const std::vector< Real > &  intnl,
const std::vector< bool > &  active,
std::vector< Real > &  dh_dintnl 
)
protectedvirtualinherited

The derivative of the active hardening potentials with respect to the active internal parameters.

Parameters
stressthe stress at which to calculate the hardening potentials
intnlvector of internal parameters
activeset of active constraints - only the active derivatives are put into "dh_dintnl"
[out]dh_dintnlthe derivatives. dh_dintnl[a][alpha][b] = dh[a][alpha]/dintnl[b]. Note that the userobjects assume that there is exactly one internal parameter per yield function, so the derivative is only nonzero for a=alpha=b, so that is all we calculate

Definition at line 315 of file MultiPlasticityRawComponentAssembler.C.

Referenced by MultiPlasticityLinearSystem::calculateJacobian().

319 {
320  mooseAssert(intnl.size() == _num_models, "Incorrect size of internal parameters");
321  mooseAssert(active.size() == _num_surfaces, "Incorrect size of active");
322 
323  dh_dintnl.resize(0);
324  std::vector<unsigned int> active_surfaces_of_model;
325  std::vector<unsigned int>::iterator active_surface;
326  std::vector<Real> model_dh_dintnl;
327  for (unsigned model = 0; model < _num_models; ++model)
328  {
329  activeModelSurfaces(model, active, active_surfaces_of_model);
330  if (active_surfaces_of_model.size() > 0)
331  {
332  _f[model]->dhardPotential_dintnlV(stress, intnl[model], model_dh_dintnl);
333  for (active_surface = active_surfaces_of_model.begin();
334  active_surface != active_surfaces_of_model.end();
335  ++active_surface)
336  dh_dintnl.push_back(model_dh_dintnl[*active_surface]);
337  }
338  }
339 }
void activeModelSurfaces(int model, const std::vector< bool > &active, std::vector< unsigned int > &active_surfaces_of_model)
Returns the internal surface number(s) of the active surfaces of the given model This may be of size=...
unsigned int _num_surfaces
Number of surfaces within the plastic models.
const PertinentGeochemicalSystem model(database, {"H2O", "H+", "HCO3-", "O2(aq)", "Ca++", ">(s)FeOH", "radius_neg1", "radius_neg1.5"}, {"Calcite"}, {}, {"Calcite_asdf"}, {"CH4(aq)"}, {">(s)FeOCa+"}, "O2(aq)", "e-")
unsigned int _num_models
Number of plastic models for this material.
std::vector< const SolidMechanicsPlasticModel * > _f
User objects that define the yield functions, flow potentials, etc.

◆ dhardPotential_dstress()

void MultiPlasticityRawComponentAssembler::dhardPotential_dstress ( const RankTwoTensor stress,
const std::vector< Real > &  intnl,
const std::vector< bool > &  active,
std::vector< RankTwoTensor > &  dh_dstress 
)
protectedvirtualinherited

The derivative of the active hardening potentials with respect to stress By assumption in the Userobjects, the h[a][alpha] is nonzero only for a = alpha, so we only calculate those here.

Parameters
stressthe stress at which to calculate the hardening potentials
intnlvector of internal parameters
activeset of active constraints - only the active derivatives are put into "dh_dstress"
[out]dh_dstressthe derivative. dh_dstress[a](i, j) = dh[a]/dstress(k, l)

Definition at line 287 of file MultiPlasticityRawComponentAssembler.C.

Referenced by MultiPlasticityLinearSystem::calculateJacobian().

292 {
293  mooseAssert(intnl.size() == _num_models, "Incorrect size of internal parameters");
294  mooseAssert(active.size() == _num_surfaces, "Incorrect size of active");
295 
296  dh_dstress.resize(0);
297  std::vector<unsigned int> active_surfaces_of_model;
298  std::vector<unsigned int>::iterator active_surface;
299  std::vector<RankTwoTensor> model_dh_dstress;
300  for (unsigned model = 0; model < _num_models; ++model)
301  {
302  activeModelSurfaces(model, active, active_surfaces_of_model);
303  if (active_surfaces_of_model.size() > 0)
304  {
305  _f[model]->dhardPotential_dstressV(stress, intnl[model], model_dh_dstress);
306  for (active_surface = active_surfaces_of_model.begin();
307  active_surface != active_surfaces_of_model.end();
308  ++active_surface)
309  dh_dstress.push_back(model_dh_dstress[*active_surface]);
310  }
311  }
312 }
void activeModelSurfaces(int model, const std::vector< bool > &active, std::vector< unsigned int > &active_surfaces_of_model)
Returns the internal surface number(s) of the active surfaces of the given model This may be of size=...
unsigned int _num_surfaces
Number of surfaces within the plastic models.
const PertinentGeochemicalSystem model(database, {"H2O", "H+", "HCO3-", "O2(aq)", "Ca++", ">(s)FeOH", "radius_neg1", "radius_neg1.5"}, {"Calcite"}, {}, {"Calcite_asdf"}, {"CH4(aq)"}, {">(s)FeOCa+"}, "O2(aq)", "e-")
unsigned int _num_models
Number of plastic models for this material.
std::vector< const SolidMechanicsPlasticModel * > _f
User objects that define the yield functions, flow potentials, etc.

◆ dof_included()

bool MultiPlasticityDebugger::dof_included ( unsigned int  dof,
const std::vector< bool > &  deactivated_due_to_ld 
)
private

Definition at line 353 of file MultiPlasticityDebugger.C.

Referenced by fdJacobian().

355 {
356  if (dof < unsigned(6))
357  // these are the stress components
358  return true;
359  unsigned eff_dof = dof - 6;
360  if (eff_dof < _num_surfaces)
361  // these are the plastic multipliers, pm
362  return !deactivated_due_to_ld[eff_dof];
363  eff_dof -= _num_surfaces; // now we know the dof is an intnl parameter
364  std::vector<bool> active_surface(_num_surfaces);
365  for (unsigned surface = 0; surface < _num_surfaces; ++surface)
366  active_surface[surface] = !deactivated_due_to_ld[surface];
367  return anyActiveSurfaces(eff_dof, active_surface);
368 }
unsigned int _num_surfaces
Number of surfaces within the plastic models.
bool anyActiveSurfaces(int model, const std::vector< bool > &active)
returns true if any internal surfaces of the given model are active according to &#39;active&#39; ...

◆ dyieldFunction_dintnl()

void MultiPlasticityRawComponentAssembler::dyieldFunction_dintnl ( const RankTwoTensor stress,
const std::vector< Real > &  intnl,
const std::vector< bool > &  active,
std::vector< Real > &  df_dintnl 
)
protectedvirtualinherited

The derivative of active yield function(s) with respect to their internal parameters (the user objects assume there is exactly one internal param per yield function)

Parameters
stressthe stress at which to calculate the yield function
intnlvector of internal parameters
activeset of active constraints - only the active derivatives are put into "df_dintnl"
[out]df_dintnlthe derivatives. df_dstress[alpha] = dyieldFunction[alpha]/dintnl[alpha]

Definition at line 151 of file MultiPlasticityRawComponentAssembler.C.

Referenced by MultiPlasticityLinearSystem::calculateJacobian(), checkDerivatives(), and ComputeMultiPlasticityStress::consistentTangentOperator().

155 {
156  mooseAssert(intnl.size() == _num_models, "Incorrect size of internal parameters");
157  mooseAssert(active.size() == _num_surfaces, "Incorrect size of active");
158 
159  df_dintnl.resize(0);
160  std::vector<unsigned int> active_surfaces_of_model;
161  std::vector<unsigned int>::iterator active_surface;
162  std::vector<Real> model_df_dintnl;
163  for (unsigned model = 0; model < _num_models; ++model)
164  {
165  activeModelSurfaces(model, active, active_surfaces_of_model);
166  if (active_surfaces_of_model.size() > 0)
167  {
168  _f[model]->dyieldFunction_dintnlV(stress, intnl[model], model_df_dintnl);
169  for (active_surface = active_surfaces_of_model.begin();
170  active_surface != active_surfaces_of_model.end();
171  ++active_surface)
172  df_dintnl.push_back(model_df_dintnl[*active_surface]);
173  }
174  }
175 }
void activeModelSurfaces(int model, const std::vector< bool > &active, std::vector< unsigned int > &active_surfaces_of_model)
Returns the internal surface number(s) of the active surfaces of the given model This may be of size=...
unsigned int _num_surfaces
Number of surfaces within the plastic models.
const PertinentGeochemicalSystem model(database, {"H2O", "H+", "HCO3-", "O2(aq)", "Ca++", ">(s)FeOH", "radius_neg1", "radius_neg1.5"}, {"Calcite"}, {}, {"Calcite_asdf"}, {"CH4(aq)"}, {">(s)FeOCa+"}, "O2(aq)", "e-")
unsigned int _num_models
Number of plastic models for this material.
std::vector< const SolidMechanicsPlasticModel * > _f
User objects that define the yield functions, flow potentials, etc.

◆ dyieldFunction_dstress()

void MultiPlasticityRawComponentAssembler::dyieldFunction_dstress ( const RankTwoTensor stress,
const std::vector< Real > &  intnl,
const std::vector< bool > &  active,
std::vector< RankTwoTensor > &  df_dstress 
)
protectedvirtualinherited

The derivative of the active yield function(s) with respect to stress.

Parameters
stressthe stress at which to calculate the yield function
intnlvector of internal parameters
activeset of active constraints - only the active derivatives are put into "df_dstress"
[out]df_dstressthe derivative (or derivatives in the case of multisurface plasticity). df_dstress[alpha](i, j) = dyieldFunction[alpha]/dstress(i, j)

Definition at line 123 of file MultiPlasticityRawComponentAssembler.C.

Referenced by ComputeMultiPlasticityStress::buildDumbOrder(), MultiPlasticityLinearSystem::calculateJacobian(), checkDerivatives(), ComputeMultiPlasticityStress::consistentTangentOperator(), and MultiPlasticityLinearSystem::eliminateLinearDependence().

128 {
129  mooseAssert(intnl.size() == _num_models, "Incorrect size of internal parameters");
130  mooseAssert(active.size() == _num_surfaces, "Incorrect size of active");
131 
132  df_dstress.resize(0);
133  std::vector<unsigned int> active_surfaces_of_model;
134  std::vector<unsigned int>::iterator active_surface;
135  std::vector<RankTwoTensor> model_df_dstress;
136  for (unsigned model = 0; model < _num_models; ++model)
137  {
138  activeModelSurfaces(model, active, active_surfaces_of_model);
139  if (active_surfaces_of_model.size() > 0)
140  {
141  _f[model]->dyieldFunction_dstressV(stress, intnl[model], model_df_dstress);
142  for (active_surface = active_surfaces_of_model.begin();
143  active_surface != active_surfaces_of_model.end();
144  ++active_surface)
145  df_dstress.push_back(model_df_dstress[*active_surface]);
146  }
147  }
148 }
void activeModelSurfaces(int model, const std::vector< bool > &active, std::vector< unsigned int > &active_surfaces_of_model)
Returns the internal surface number(s) of the active surfaces of the given model This may be of size=...
unsigned int _num_surfaces
Number of surfaces within the plastic models.
const PertinentGeochemicalSystem model(database, {"H2O", "H+", "HCO3-", "O2(aq)", "Ca++", ">(s)FeOH", "radius_neg1", "radius_neg1.5"}, {"Calcite"}, {}, {"Calcite_asdf"}, {"CH4(aq)"}, {">(s)FeOCa+"}, "O2(aq)", "e-")
unsigned int _num_models
Number of plastic models for this material.
std::vector< const SolidMechanicsPlasticModel * > _f
User objects that define the yield functions, flow potentials, etc.

◆ fddflowPotential_dintnl()

void MultiPlasticityDebugger::fddflowPotential_dintnl ( const RankTwoTensor stress,
const std::vector< Real > &  intnl,
std::vector< RankTwoTensor > &  dr_dintnl 
)
privatevirtual

The finite-difference derivative of the flow potentials with respect to internal parameters.

Parameters
stressthe stress at which to calculate the flow potential
intnlvector of internal parameters
[out]dr_dintnlthe derivatives. dr_dintnl[alpha](i, j) = dr[alpha](i, j)/dintnl[alpha]

Definition at line 612 of file MultiPlasticityDebugger.C.

Referenced by checkDerivatives().

615 {
616  dr_dintnl.resize(_num_surfaces);
617 
618  std::vector<bool> act;
619  act.assign(_num_surfaces, true);
620 
621  std::vector<RankTwoTensor> origr;
622  flowPotential(stress, intnl, act, origr);
623 
624  std::vector<Real> intnlep;
625  intnlep.resize(_num_models);
626  for (unsigned model = 0; model < _num_models; ++model)
627  intnlep[model] = intnl[model];
628  Real ep;
629  std::vector<RankTwoTensor> rep;
630  unsigned int model;
631  for (unsigned surface = 0; surface < _num_surfaces; ++surface)
632  {
633  model = modelNumber(surface);
635  intnlep[model] += ep;
636  flowPotential(stress, intnlep, act, rep);
637  dr_dintnl[surface] = (rep[surface] - origr[surface]) / ep;
638  intnlep[model] -= ep;
639  }
640 }
unsigned int modelNumber(unsigned int surface)
returns the model number, given the surface number
std::vector< Real > _fspb_debug_intnl_change
Debug finite-differencing parameters for the internal parameters.
unsigned int _num_surfaces
Number of surfaces within the plastic models.
virtual void flowPotential(const RankTwoTensor &stress, const std::vector< Real > &intnl, const std::vector< bool > &active, std::vector< RankTwoTensor > &r)
The active flow potential(s) - one for each yield function.
const PertinentGeochemicalSystem model(database, {"H2O", "H+", "HCO3-", "O2(aq)", "Ca++", ">(s)FeOH", "radius_neg1", "radius_neg1.5"}, {"Calcite"}, {}, {"Calcite_asdf"}, {"CH4(aq)"}, {">(s)FeOCa+"}, "O2(aq)", "e-")
unsigned int _num_models
Number of plastic models for this material.
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real

◆ fddflowPotential_dstress()

void MultiPlasticityDebugger::fddflowPotential_dstress ( const RankTwoTensor stress,
const std::vector< Real > &  intnl,
std::vector< RankFourTensor > &  dr_dstress 
)
privatevirtual

The finite-difference derivative of the flow potential(s) with respect to stress.

Parameters
stressthe stress at which to calculate the flow potential
intnlvector of internal parameters
[out]dr_dstressthe derivative. dr_dstress[alpha](i, j, k, l) = dr[alpha](i, j)/dstress(k, l)

Definition at line 583 of file MultiPlasticityDebugger.C.

Referenced by checkDerivatives().

586 {
587  dr_dstress.assign(_num_surfaces, RankFourTensor());
588 
589  std::vector<bool> act;
590  act.assign(_num_surfaces, true);
591 
593  RankTwoTensor stressep;
594  std::vector<RankTwoTensor> rep, rep_minus;
595  for (unsigned i = 0; i < 3; ++i)
596  for (unsigned j = 0; j < 3; ++j)
597  {
598  stressep = stress;
599  // do a central difference
600  stressep(i, j) += ep / 2.0;
601  flowPotential(stressep, intnl, act, rep);
602  stressep(i, j) -= ep;
603  flowPotential(stressep, intnl, act, rep_minus);
604  for (unsigned surface = 0; surface < _num_surfaces; ++surface)
605  for (unsigned k = 0; k < 3; ++k)
606  for (unsigned l = 0; l < 3; ++l)
607  dr_dstress[surface](k, l, i, j) = (rep[surface](k, l) - rep_minus[surface](k, l)) / ep;
608  }
609 }
Real _fspb_debug_stress_change
Debug finite-differencing parameter for the stress.
unsigned int _num_surfaces
Number of surfaces within the plastic models.
virtual void flowPotential(const RankTwoTensor &stress, const std::vector< Real > &intnl, const std::vector< bool > &active, std::vector< RankTwoTensor > &r)
The active flow potential(s) - one for each yield function.
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
static const std::complex< double > j(0, 1)
Complex number "j" (also known as "i")
static const std::string k
Definition: NS.h:124

◆ fddyieldFunction_dintnl()

void MultiPlasticityDebugger::fddyieldFunction_dintnl ( const RankTwoTensor stress,
const std::vector< Real > &  intnl,
std::vector< Real > &  df_dintnl 
)
private

The finite-difference derivative of yield function(s) with respect to internal parameter(s)

Parameters
stressthe stress at which to calculate the yield function
intnlvector of internal parameters
[out]df_dintnlthe derivative (or derivatives in the case of multisurface plasticity). df_dintnl[alpha] = dyieldFunction[alpha]/dintnl[alpha]

Definition at line 552 of file MultiPlasticityDebugger.C.

Referenced by checkDerivatives().

555 {
556  df_dintnl.resize(_num_surfaces);
557 
558  std::vector<bool> act;
559  act.assign(_num_surfaces, true);
560 
561  std::vector<Real> origf;
562  yieldFunction(stress, intnl, act, origf);
563 
564  std::vector<Real> intnlep;
565  intnlep.resize(_num_models);
566  for (unsigned model = 0; model < _num_models; ++model)
567  intnlep[model] = intnl[model];
568  Real ep;
569  std::vector<Real> fep;
570  unsigned int model;
571  for (unsigned surface = 0; surface < _num_surfaces; ++surface)
572  {
573  model = modelNumber(surface);
575  intnlep[model] += ep;
576  yieldFunction(stress, intnlep, act, fep);
577  df_dintnl[surface] = (fep[surface] - origf[surface]) / ep;
578  intnlep[model] -= ep;
579  }
580 }
unsigned int modelNumber(unsigned int surface)
returns the model number, given the surface number
std::vector< Real > _fspb_debug_intnl_change
Debug finite-differencing parameters for the internal parameters.
virtual void yieldFunction(const RankTwoTensor &stress, const std::vector< Real > &intnl, const std::vector< bool > &active, std::vector< Real > &f)
The active yield function(s)
unsigned int _num_surfaces
Number of surfaces within the plastic models.
const PertinentGeochemicalSystem model(database, {"H2O", "H+", "HCO3-", "O2(aq)", "Ca++", ">(s)FeOH", "radius_neg1", "radius_neg1.5"}, {"Calcite"}, {}, {"Calcite_asdf"}, {"CH4(aq)"}, {">(s)FeOCa+"}, "O2(aq)", "e-")
unsigned int _num_models
Number of plastic models for this material.
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real

◆ fddyieldFunction_dstress()

void MultiPlasticityDebugger::fddyieldFunction_dstress ( const RankTwoTensor stress,
const std::vector< Real > &  intnl,
std::vector< RankTwoTensor > &  df_dstress 
)
private

The finite-difference derivative of yield function(s) with respect to stress.

Parameters
stressthe stress at which to calculate the yield function
intnlvector of internal parameters
[out]df_dstressthe derivative (or derivatives in the case of multisurface plasticity). df_dstress[alpha](i, j) = dyieldFunction[alpha]/dstress(i, j)

Definition at line 524 of file MultiPlasticityDebugger.C.

Referenced by checkDerivatives().

527 {
528  df_dstress.assign(_num_surfaces, RankTwoTensor());
529 
530  std::vector<bool> act;
531  act.assign(_num_surfaces, true);
532 
534  RankTwoTensor stressep;
535  std::vector<Real> fep, fep_minus;
536  for (unsigned i = 0; i < 3; ++i)
537  for (unsigned j = 0; j < 3; ++j)
538  {
539  stressep = stress;
540  // do a central difference to attempt to capture discontinuities
541  // such as those encountered in tensile and Mohr-Coulomb
542  stressep(i, j) += ep / 2.0;
543  yieldFunction(stressep, intnl, act, fep);
544  stressep(i, j) -= ep;
545  yieldFunction(stressep, intnl, act, fep_minus);
546  for (unsigned surface = 0; surface < _num_surfaces; ++surface)
547  df_dstress[surface](i, j) = (fep[surface] - fep_minus[surface]) / ep;
548  }
549 }
Real _fspb_debug_stress_change
Debug finite-differencing parameter for the stress.
virtual void yieldFunction(const RankTwoTensor &stress, const std::vector< Real > &intnl, const std::vector< bool > &active, std::vector< Real > &f)
The active yield function(s)
unsigned int _num_surfaces
Number of surfaces within the plastic models.
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
static const std::complex< double > j(0, 1)
Complex number "j" (also known as "i")

◆ fdJacobian()

void MultiPlasticityDebugger::fdJacobian ( const RankTwoTensor stress,
const std::vector< Real > &  intnl_old,
const std::vector< Real > &  intnl,
const std::vector< Real > &  pm,
const RankTwoTensor delta_dp,
const RankFourTensor E_inv,
bool  eliminate_ld,
std::vector< std::vector< Real >> &  jac 
)
privatevirtual

The Jacobian calculated using finite differences.

The output should be equal to calculateJacobian(...) if everything is coded correctly.

Parameters
stressthe stress at which to calculate the Jacobian
intnl_oldthe old values of internal variables (jacobian is inependent of these, but they are needed to do the finite-differencing cleanly)
intnlthe vector of internal parameters at which to calculate the Jacobian
pmthe plasticity multipliers at which to calculate the Jacobian
delta_dpplastic_strain - plastic_strain_old (Jacobian is independent of this, but it is needed to do the finite-differencing cleanly)
E_invinverse of the elasticity tensor
eliminate_ldonly calculate the Jacobian for the linearly independent constraints
[out]jacthe finite-difference Jacobian

Definition at line 225 of file MultiPlasticityDebugger.C.

Referenced by checkJacobian(), and checkSolution().

233 {
234  std::vector<bool> active;
235  active.assign(_num_surfaces, true);
236 
237  std::vector<bool> deactivated_due_to_ld;
238  std::vector<bool> deactivated_due_to_ld_ep;
239 
240  std::vector<Real> orig_rhs;
241  calculateRHS(stress,
242  intnl_old,
243  intnl,
244  pm,
245  delta_dp,
246  orig_rhs,
247  active,
248  eliminate_ld,
249  deactivated_due_to_ld); // this calculates RHS, and also set deactivated_due_to_ld.
250  // The latter stays fixed for the rest of this routine
251 
252  unsigned int whole_system_size = 6 + _num_surfaces + _num_models;
253  unsigned int system_size =
254  orig_rhs.size(); // will be = whole_system_size if eliminate_ld = false, since all active=true
255  jac.resize(system_size);
256  for (unsigned row = 0; row < system_size; ++row)
257  jac[row].assign(system_size, 0);
258 
259  std::vector<Real> rhs_ep;
260  unsigned col = 0;
261 
262  RankTwoTensor stressep;
263  RankTwoTensor delta_dpep;
265  for (unsigned i = 0; i < 3; ++i)
266  for (unsigned j = 0; j <= i; ++j)
267  {
268  stressep = stress;
269  stressep(i, j) += ep;
270  if (i != j)
271  stressep(j, i) += ep;
272  delta_dpep = delta_dp;
273  for (unsigned k = 0; k < 3; ++k)
274  for (unsigned l = 0; l < 3; ++l)
275  {
276  delta_dpep(k, l) -= E_inv(k, l, i, j) * ep;
277  if (i != j)
278  delta_dpep(k, l) -= E_inv(k, l, j, i) * ep;
279  }
280  active.assign(_num_surfaces, true);
281  calculateRHS(stressep,
282  intnl_old,
283  intnl,
284  pm,
285  delta_dpep,
286  rhs_ep,
287  active,
288  false,
289  deactivated_due_to_ld_ep);
290  unsigned row = 0;
291  for (unsigned dof = 0; dof < whole_system_size; ++dof)
292  if (dof_included(dof, deactivated_due_to_ld))
293  {
294  jac[row][col] =
295  -(rhs_ep[dof] - orig_rhs[row]) / ep; // remember jacobian = -d(rhs)/d(something)
296  row++;
297  }
298  col++; // all of the first 6 columns are dof_included since they're stresses
299  }
300 
301  std::vector<Real> pmep;
302  pmep.resize(_num_surfaces);
303  for (unsigned surface = 0; surface < _num_surfaces; ++surface)
304  pmep[surface] = pm[surface];
305  for (unsigned surface = 0; surface < _num_surfaces; ++surface)
306  {
307  if (!dof_included(6 + surface, deactivated_due_to_ld))
308  continue;
309  ep = _fspb_debug_pm_change[surface];
310  pmep[surface] += ep;
311  active.assign(_num_surfaces, true);
312  calculateRHS(
313  stress, intnl_old, intnl, pmep, delta_dp, rhs_ep, active, false, deactivated_due_to_ld_ep);
314  unsigned row = 0;
315  for (unsigned dof = 0; dof < whole_system_size; ++dof)
316  if (dof_included(dof, deactivated_due_to_ld))
317  {
318  jac[row][col] =
319  -(rhs_ep[dof] - orig_rhs[row]) / ep; // remember jacobian = -d(rhs)/d(something)
320  row++;
321  }
322  pmep[surface] -= ep;
323  col++;
324  }
325 
326  std::vector<Real> intnlep;
327  intnlep.resize(_num_models);
328  for (unsigned model = 0; model < _num_models; ++model)
329  intnlep[model] = intnl[model];
330  for (unsigned model = 0; model < _num_models; ++model)
331  {
332  if (!dof_included(6 + _num_surfaces + model, deactivated_due_to_ld))
333  continue;
335  intnlep[model] += ep;
336  active.assign(_num_surfaces, true);
337  calculateRHS(
338  stress, intnl_old, intnlep, pm, delta_dp, rhs_ep, active, false, deactivated_due_to_ld_ep);
339  unsigned row = 0;
340  for (unsigned dof = 0; dof < whole_system_size; ++dof)
341  if (dof_included(dof, deactivated_due_to_ld))
342  {
343  jac[row][col] =
344  -(rhs_ep[dof] - orig_rhs[row]) / ep; // remember jacobian = -d(rhs)/d(something)
345  row++;
346  }
347  intnlep[model] -= ep;
348  col++;
349  }
350 }
virtual void calculateRHS(const RankTwoTensor &stress, const std::vector< Real > &intnl_old, const std::vector< Real > &intnl, const std::vector< Real > &pm, const RankTwoTensor &delta_dp, std::vector< Real > &rhs, const std::vector< bool > &active, bool eliminate_ld, std::vector< bool > &deactivated_due_to_ld)
Calculate the RHS which is rhs = -(epp(0,0), epp(1,0), epp(1,1), epp(2,0), epp(2,1), epp(2,2), f[0], f[1], ..., f[num_f], ic[0], ic[1], ..., ic[num_ic])
Real _fspb_debug_stress_change
Debug finite-differencing parameter for the stress.
std::vector< Real > _fspb_debug_intnl_change
Debug finite-differencing parameters for the internal parameters.
unsigned int _num_surfaces
Number of surfaces within the plastic models.
const PertinentGeochemicalSystem model(database, {"H2O", "H+", "HCO3-", "O2(aq)", "Ca++", ">(s)FeOH", "radius_neg1", "radius_neg1.5"}, {"Calcite"}, {}, {"Calcite_asdf"}, {"CH4(aq)"}, {">(s)FeOCa+"}, "O2(aq)", "e-")
bool dof_included(unsigned int dof, const std::vector< bool > &deactivated_due_to_ld)
std::vector< Real > _fspb_debug_pm_change
Debug finite-differencing parameters for the plastic multipliers.
unsigned int _num_models
Number of plastic models for this material.
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
static const std::complex< double > j(0, 1)
Complex number "j" (also known as "i")
static const std::string k
Definition: NS.h:124

◆ flowPotential()

void MultiPlasticityRawComponentAssembler::flowPotential ( const RankTwoTensor stress,
const std::vector< Real > &  intnl,
const std::vector< bool > &  active,
std::vector< RankTwoTensor > &  r 
)
protectedvirtualinherited

The active flow potential(s) - one for each yield function.

Parameters
stressthe stress at which to calculate the flow potential
intnlvector of internal parameters
activeset of active constraints - only the active flow potentials are put into "r"
[out]rthe flow potential (flow potentials in the multi-surface case)

Definition at line 178 of file MultiPlasticityRawComponentAssembler.C.

Referenced by MultiPlasticityLinearSystem::calculateConstraints(), MultiPlasticityLinearSystem::calculateJacobian(), ComputeMultiPlasticityStress::consistentTangentOperator(), fddflowPotential_dintnl(), and fddflowPotential_dstress().

182 {
183  mooseAssert(intnl.size() == _num_models, "Incorrect size of internal parameters");
184  mooseAssert(active.size() == _num_surfaces, "Incorrect size of active");
185 
186  r.resize(0);
187  std::vector<unsigned int> active_surfaces_of_model;
188  std::vector<unsigned int>::iterator active_surface;
189  std::vector<RankTwoTensor> model_r;
190  for (unsigned model = 0; model < _num_models; ++model)
191  {
192  activeModelSurfaces(model, active, active_surfaces_of_model);
193  if (active_surfaces_of_model.size() > 0)
194  {
195  _f[model]->flowPotentialV(stress, intnl[model], model_r);
196  for (active_surface = active_surfaces_of_model.begin();
197  active_surface != active_surfaces_of_model.end();
198  ++active_surface)
199  r.push_back(model_r[*active_surface]);
200  }
201  }
202 }
void activeModelSurfaces(int model, const std::vector< bool > &active, std::vector< unsigned int > &active_surfaces_of_model)
Returns the internal surface number(s) of the active surfaces of the given model This may be of size=...
unsigned int _num_surfaces
Number of surfaces within the plastic models.
const PertinentGeochemicalSystem model(database, {"H2O", "H+", "HCO3-", "O2(aq)", "Ca++", ">(s)FeOH", "radius_neg1", "radius_neg1.5"}, {"Calcite"}, {}, {"Calcite_asdf"}, {"CH4(aq)"}, {">(s)FeOCa+"}, "O2(aq)", "e-")
unsigned int _num_models
Number of plastic models for this material.
std::vector< const SolidMechanicsPlasticModel * > _f
User objects that define the yield functions, flow potentials, etc.

◆ hardPotential()

void MultiPlasticityRawComponentAssembler::hardPotential ( const RankTwoTensor stress,
const std::vector< Real > &  intnl,
const std::vector< bool > &  active,
std::vector< Real > &  h 
)
protectedvirtualinherited

The active hardening potentials (one for each internal parameter and for each yield function) by assumption in the Userobjects, the h[a][alpha] is nonzero only if the surface alpha is part of model a, so we only calculate those here.

Parameters
stressthe stress at which to calculate the hardening potential
intnlvector of internal parameters
activeset of active constraints - only the active hardening potentials are put into "h"
[out]hthe hardening potentials. h[alpha] = hardening potential for yield fcn alpha (and, by the above assumption we know which hardening parameter, a, this belongs to)

Definition at line 260 of file MultiPlasticityRawComponentAssembler.C.

Referenced by MultiPlasticityLinearSystem::calculateConstraints(), MultiPlasticityLinearSystem::calculateJacobian(), and ComputeMultiPlasticityStress::consistentTangentOperator().

264 {
265  mooseAssert(intnl.size() == _num_models, "Incorrect size of internal parameters");
266  mooseAssert(active.size() == _num_surfaces, "Incorrect size of active");
267 
268  h.resize(0);
269  std::vector<unsigned int> active_surfaces_of_model;
270  std::vector<unsigned int>::iterator active_surface;
271  std::vector<Real> model_h;
272  for (unsigned model = 0; model < _num_models; ++model)
273  {
274  activeModelSurfaces(model, active, active_surfaces_of_model);
275  if (active_surfaces_of_model.size() > 0)
276  {
277  _f[model]->hardPotentialV(stress, intnl[model], model_h);
278  for (active_surface = active_surfaces_of_model.begin();
279  active_surface != active_surfaces_of_model.end();
280  ++active_surface)
281  h.push_back(model_h[*active_surface]);
282  }
283  }
284 }
void activeModelSurfaces(int model, const std::vector< bool > &active, std::vector< unsigned int > &active_surfaces_of_model)
Returns the internal surface number(s) of the active surfaces of the given model This may be of size=...
unsigned int _num_surfaces
Number of surfaces within the plastic models.
const PertinentGeochemicalSystem model(database, {"H2O", "H+", "HCO3-", "O2(aq)", "Ca++", ">(s)FeOH", "radius_neg1", "radius_neg1.5"}, {"Calcite"}, {}, {"Calcite_asdf"}, {"CH4(aq)"}, {">(s)FeOCa+"}, "O2(aq)", "e-")
unsigned int _num_models
Number of plastic models for this material.
std::vector< const SolidMechanicsPlasticModel * > _f
User objects that define the yield functions, flow potentials, etc.

◆ modelNumber()

unsigned int MultiPlasticityRawComponentAssembler::modelNumber ( unsigned int  surface)
protectedinherited

◆ nrStep()

void MultiPlasticityLinearSystem::nrStep ( const RankTwoTensor stress,
const std::vector< Real > &  intnl_old,
const std::vector< Real > &  intnl,
const std::vector< Real > &  pm,
const RankFourTensor E_inv,
const RankTwoTensor delta_dp,
RankTwoTensor dstress,
std::vector< Real > &  dpm,
std::vector< Real > &  dintnl,
const std::vector< bool > &  active,
std::vector< bool > &  deactivated_due_to_ld 
)
protectedvirtualinherited

Performs one Newton-Raphson step.

The purpose here is to find the changes, dstress, dpm and dintnl according to the Newton-Raphson procedure

Parameters
stressCurrent value of stress
intnl_oldThe internal variables at the previous "time" step
intnlCurrent value of the internal variables
pmCurrent value of the plasticity multipliers (consistency parameters)
E_invinverse of the elasticity tensor
delta_dpCurrent value of the plastic-strain increment (ie plastic_strain - plastic_strain_old)
[out]dstressThe change in stress for a full Newton step
[out]dpmThe change in all plasticity multipliers for a full Newton step
[out]dintnlThe change in all internal variables for a full Newton step
activeThe active constraints
[out]deactivated_due_to_ldThe constraints deactivated due to linear-dependence of the flow directions

Definition at line 614 of file MultiPlasticityLinearSystem.C.

Referenced by checkSolution(), and ComputeMultiPlasticityStress::singleStep().

625 {
626  // Calculate RHS and Jacobian
627  std::vector<Real> rhs;
628  calculateRHS(stress, intnl_old, intnl, pm, delta_dp, rhs, active, true, deactivated_due_to_ld);
629 
630  std::vector<std::vector<Real>> jac;
631  calculateJacobian(stress, intnl, pm, E_inv, active, deactivated_due_to_ld, jac);
632 
633  // prepare for LAPACKgesv_ routine provided by PETSc
634  PetscBLASInt system_size = rhs.size();
635 
636  std::vector<double> a(system_size * system_size);
637  // Fill in the a "matrix" by going down columns
638  unsigned ind = 0;
639  for (int col = 0; col < system_size; ++col)
640  for (int row = 0; row < system_size; ++row)
641  a[ind++] = jac[row][col];
642 
643  PetscBLASInt nrhs = 1;
644  std::vector<PetscBLASInt> ipiv(system_size);
645  PetscBLASInt info;
646  LAPACKgesv_(&system_size, &nrhs, &a[0], &system_size, &ipiv[0], &rhs[0], &system_size, &info);
647 
648  if (info != 0)
649  mooseError("In solving the linear system in a Newton-Raphson process, the PETSC LAPACK gsev "
650  "routine returned with error code ",
651  info);
652 
653  // Extract the results back to dstress, dpm and dintnl
654  std::vector<bool> active_not_deact(_num_surfaces);
655  for (unsigned surface = 0; surface < _num_surfaces; ++surface)
656  active_not_deact[surface] = (active[surface] && !deactivated_due_to_ld[surface]);
657 
658  unsigned int dim = 3;
659  ind = 0;
660 
661  for (unsigned i = 0; i < dim; ++i)
662  for (unsigned j = 0; j <= i; ++j)
663  dstress(i, j) = dstress(j, i) = rhs[ind++];
664  dpm.assign(_num_surfaces, 0);
665  for (unsigned surface = 0; surface < _num_surfaces; ++surface)
666  if (active_not_deact[surface])
667  dpm[surface] = rhs[ind++];
668  dintnl.assign(_num_models, 0);
669  for (unsigned model = 0; model < _num_models; ++model)
670  if (anyActiveSurfaces(model, active_not_deact))
671  dintnl[model] = rhs[ind++];
672 
673  mooseAssert(static_cast<int>(ind) == system_size,
674  "Incorrect extracting of changes from NR solution in nrStep");
675 }
virtual void calculateRHS(const RankTwoTensor &stress, const std::vector< Real > &intnl_old, const std::vector< Real > &intnl, const std::vector< Real > &pm, const RankTwoTensor &delta_dp, std::vector< Real > &rhs, const std::vector< bool > &active, bool eliminate_ld, std::vector< bool > &deactivated_due_to_ld)
Calculate the RHS which is rhs = -(epp(0,0), epp(1,0), epp(1,1), epp(2,0), epp(2,1), epp(2,2), f[0], f[1], ..., f[num_f], ic[0], ic[1], ..., ic[num_ic])
MPI_Info info
void mooseError(Args &&... args)
unsigned int dim
virtual void calculateJacobian(const RankTwoTensor &stress, const std::vector< Real > &intnl, const std::vector< Real > &pm, const RankFourTensor &E_inv, const std::vector< bool > &active, const std::vector< bool > &deactivated_due_to_ld, std::vector< std::vector< Real >> &jac)
d(rhs)/d(dof)
unsigned int _num_surfaces
Number of surfaces within the plastic models.
void assign(const TypeTensor< T2 > &)
bool anyActiveSurfaces(int model, const std::vector< bool > &active)
returns true if any internal surfaces of the given model are active according to &#39;active&#39; ...
const PertinentGeochemicalSystem model(database, {"H2O", "H+", "HCO3-", "O2(aq)", "Ca++", ">(s)FeOH", "radius_neg1", "radius_neg1.5"}, {"Calcite"}, {}, {"Calcite_asdf"}, {"CH4(aq)"}, {">(s)FeOCa+"}, "O2(aq)", "e-")
unsigned int _num_models
Number of plastic models for this material.
static const std::complex< double > j(0, 1)
Complex number "j" (also known as "i")

◆ outputAndCheckDebugParameters()

void MultiPlasticityDebugger::outputAndCheckDebugParameters ( )

Outputs the debug parameters: _fspb_debug_stress, _fspd_debug_pm, etc and checks that they are sized correctly.

Definition at line 53 of file MultiPlasticityDebugger.C.

Referenced by checkDerivatives(), checkJacobian(), and checkSolution().

54 {
55  Moose::err << "Debug Parameters are as follows\n";
56  Moose::err << "stress = \n";
58 
59  if (_fspb_debug_pm.size() != _num_surfaces || _fspb_debug_intnl.size() != _num_models ||
62  mooseError("The debug parameters have the wrong size\n");
63 
64  Moose::err << "plastic multipliers =\n";
65  for (unsigned surface = 0; surface < _num_surfaces; ++surface)
66  Moose::err << _fspb_debug_pm[surface] << "\n";
67 
68  Moose::err << "internal parameters =\n";
69  for (unsigned model = 0; model < _num_models; ++model)
70  Moose::err << _fspb_debug_intnl[model] << "\n";
71 
72  Moose::err << "finite-differencing parameter for stress-changes:\n"
73  << _fspb_debug_stress_change << "\n";
74  Moose::err << "finite-differencing parameter(s) for plastic-multiplier(s):\n";
75  for (unsigned surface = 0; surface < _num_surfaces; ++surface)
76  Moose::err << _fspb_debug_pm_change[surface] << "\n";
77  Moose::err << "finite-differencing parameter(s) for internal-parameter(s):\n";
78  for (unsigned model = 0; model < _num_models; ++model)
79  Moose::err << _fspb_debug_intnl_change[model] << "\n";
80 
81  Moose::err << std::flush;
82 }
void mooseError(Args &&... args)
void print(std::ostream &stm=Moose::out) const
Real _fspb_debug_stress_change
Debug finite-differencing parameter for the stress.
std::vector< Real > _fspb_debug_intnl_change
Debug finite-differencing parameters for the internal parameters.
std::vector< Real > _fspb_debug_intnl
Debug the Jacobian entires at these internal parameters.
unsigned int _num_surfaces
Number of surfaces within the plastic models.
RankTwoTensor _fspb_debug_stress
Debug the Jacobian entries at this stress.
const PertinentGeochemicalSystem model(database, {"H2O", "H+", "HCO3-", "O2(aq)", "Ca++", ">(s)FeOH", "radius_neg1", "radius_neg1.5"}, {"Calcite"}, {}, {"Calcite_asdf"}, {"CH4(aq)"}, {">(s)FeOCa+"}, "O2(aq)", "e-")
std::vector< Real > _fspb_debug_pm_change
Debug finite-differencing parameters for the plastic multipliers.
unsigned int _num_models
Number of plastic models for this material.
std::vector< Real > _fspb_debug_pm
Debug the Jacobian entires at these plastic multipliers.

◆ returnMapAll()

bool MultiPlasticityRawComponentAssembler::returnMapAll ( const RankTwoTensor trial_stress,
const std::vector< Real > &  intnl_old,
const RankFourTensor E_ijkl,
Real  ep_plastic_tolerance,
RankTwoTensor stress,
std::vector< Real > &  intnl,
std::vector< Real > &  pm,
std::vector< Real > &  cumulative_pm,
RankTwoTensor delta_dp,
std::vector< Real > &  yf,
unsigned &  num_successful_plastic_returns,
unsigned &  custom_model 
)
protectedinherited

Performs a returnMap for each plastic model using their inbuilt returnMap functions.

Performs a returnMap for each plastic model.

This may be used to quickly ascertain whether a (trial_stress, intnl_old) configuration is admissible, or whether a single model's customized returnMap function can provide a solution to the return-map problem, or whether a full Newton-Raphson approach such as implemented in ComputeMultiPlasticityStress is needed.

There are three cases mentioned below: (A) The (trial_stress, intnl_old) configuration is admissible according to all plastic models (B) The (trial_stress, intnl_old) configuration is inadmissible to exactly one plastic model, and that model can successfully use its customized returnMap function to provide a returned (stress, intnl) configuration, and that configuration is admissible according to all plastic models (C) All other cases. This includes customized returnMap functions failing, or more than one plastic_model being inadmissible, etc

Parameters
trial_stressthe trial stress
intnl_oldthe old values of the internal parameters
E_ijklthe elasticity tensor
ep_plastic_tolerancethe tolerance on the plastic strain
[out]stressis set to trial_stress in case (A) or (C), and the returned value of stress in case (B).
[out]intnlis set to intnl_old in case (A) or (C), and the returned value of intnl in case (B)
[out]pmZero in case (A) or (C), otherwise the plastic multipliers needed to bring about the returnMap in case (B)
[in/out]cumulative_pm cumulative plastic multipliers, updated in case (B), otherwise left untouched
[out]delta_dpis unchanged in case (A) or (C), and is set to the change in plastic strain in case(B)
[out]yfwill contain the yield function values at (stress, intnl)
[out]num_successful_plastic_returnswill be 0 for (A) and (C), and 1 for (B)
Returns
true in case (A) and (B), and false in case (C)

If all models actually signal "elastic" by returning true from their returnMap, and by returning model_plastically_active=0, then yf will contain the yield function values num_successful_plastic_returns will be zero intnl = intnl_old delta_dp will be unchanged from its input value stress will be set to trial_stress pm will be zero cumulative_pm will be unchanged return value will be true num_successful_plastic_returns = 0

If only one model signals "plastically active" by returning true from its returnMap, and by returning model_plastically_active=1, then yf will contain the yield function values num_successful_plastic_returns will be one intnl will be set by the returnMap algorithm delta_dp will be set by the returnMap algorithm stress will be set by the returnMap algorithm pm will be nonzero for the single model, and zero for other models cumulative_pm will be updated return value will be true num_successful_plastic_returns = 1

If >1 model signals "plastically active" or if >=1 model's returnMap fails, then yf will contain the yield function values num_successful_plastic_returns will be set appropriately intnl = intnl_old delta_dp will be unchanged from its input value stress will be set to trial_stress pm will be zero cumulative_pm will be unchanged return value will be true if all returnMap functions returned true, otherwise it will be false num_successful_plastic_returns is set appropriately

Definition at line 597 of file MultiPlasticityRawComponentAssembler.C.

Referenced by ComputeMultiPlasticityStress::quickStep().

609 {
610  mooseAssert(intnl_old.size() == _num_models,
611  "returnMapAll: Incorrect size of internal parameters");
612  mooseAssert(intnl.size() == _num_models, "returnMapAll: Incorrect size of internal parameters");
613  mooseAssert(pm.size() == _num_surfaces, "returnMapAll: Incorrect size of pm");
614 
615  num_successful_plastic_returns = 0;
616  yf.resize(0);
617  pm.assign(_num_surfaces, 0.0);
618 
619  RankTwoTensor returned_stress; // each model will give a returned_stress. if only one model is
620  // plastically active, i set stress=returned_stress, so as to
621  // record this returned value
622  std::vector<Real> model_f;
623  RankTwoTensor model_delta_dp;
624  std::vector<Real> model_pm;
625  bool trial_stress_inadmissible;
626  bool successful_return = true;
627  unsigned the_single_plastic_model = 0;
628  bool using_custom_return_map = true;
629 
630  // run through all the plastic models, performing their
631  // returnMap algorithms.
632  // If one finds (trial_stress, intnl) inadmissible and
633  // successfully returns, break from the loop to evaluate
634  // all the others at that returned stress
635  for (unsigned model = 0; model < _num_models; ++model)
636  {
637  if (using_custom_return_map)
638  {
639  model_pm.assign(_f[model]->numberSurfaces(), 0.0);
640  bool model_returned = _f[model]->returnMap(trial_stress,
641  intnl_old[model],
642  E_ijkl,
643  ep_plastic_tolerance,
644  returned_stress,
645  intnl[model],
646  model_pm,
647  model_delta_dp,
648  model_f,
649  trial_stress_inadmissible);
650  if (!trial_stress_inadmissible)
651  {
652  // in the elastic zone: record the yield-function values (returned_stress, intnl, model_pm
653  // and model_delta_dp are undefined)
654  for (unsigned model_surface = 0; model_surface < _f[model]->numberSurfaces();
655  ++model_surface)
656  yf.push_back(model_f[model_surface]);
657  }
658  else if (trial_stress_inadmissible && !model_returned)
659  {
660  // in the plastic zone, and the customized returnMap failed
661  // for some reason (or wasn't implemented). The coder
662  // should have correctly returned model_f(trial_stress, intnl_old)
663  // so record them
664  // (returned_stress, intnl, model_pm and model_delta_dp are undefined)
665  for (unsigned model_surface = 0; model_surface < _f[model]->numberSurfaces();
666  ++model_surface)
667  yf.push_back(model_f[model_surface]);
668  // now there's almost zero point in using the custom
669  // returnMap functions
670  using_custom_return_map = false;
671  successful_return = false;
672  }
673  else
674  {
675  // in the plastic zone, and the customized returnMap
676  // succeeded.
677  // record the first returned_stress and delta_dp if everything is going OK
678  // as they could be the actual answer
679  if (trial_stress_inadmissible)
680  num_successful_plastic_returns++;
681  the_single_plastic_model = model;
682  stress = returned_stress;
683  // note that i break here, and don't push_back
684  // model_f to yf. So now yf contains only the values of
685  // model_f from previous models to the_single_plastic_model
686  // also i don't set delta_dp = model_delta_dp yet, because
687  // i might find problems later on
688  // also, don't increment cumulative_pm for the same reason
689 
690  break;
691  }
692  }
693  else
694  {
695  // not using custom returnMap functions because one
696  // has already failed and that one said trial_stress
697  // was inadmissible. So now calculate the yield functions
698  // at the trial stress
699  _f[model]->yieldFunctionV(trial_stress, intnl_old[model], model_f);
700  for (unsigned model_surface = 0; model_surface < _f[model]->numberSurfaces(); ++model_surface)
701  yf.push_back(model_f[model_surface]);
702  }
703  }
704 
705  if (num_successful_plastic_returns == 0)
706  {
707  // here either all the models were elastic (successful_return=true),
708  // or some were plastic and either the customized returnMap failed
709  // or wasn't implemented (successful_return=false).
710  // In either case, have to set the following:
711  stress = trial_stress;
712  for (unsigned model = 0; model < _num_models; ++model)
713  intnl[model] = intnl_old[model];
714  return successful_return;
715  }
716 
717  // Now we know that num_successful_plastic_returns == 1 and all the other
718  // models (with model number < the_single_plastic_model) must have been
719  // admissible at (trial_stress, intnl). However, all models might
720  // not be admissible at (trial_stress, intnl), so must check that
721  std::vector<Real> yf_at_returned_stress(0);
722  bool all_admissible = true;
723  for (unsigned model = 0; model < _num_models; ++model)
724  {
725  if (model == the_single_plastic_model)
726  {
727  // no need to spend time calculating the yield function: we know its admissible
728  for (unsigned model_surface = 0; model_surface < _f[model]->numberSurfaces(); ++model_surface)
729  yf_at_returned_stress.push_back(model_f[model_surface]);
730  continue;
731  }
732  _f[model]->yieldFunctionV(stress, intnl_old[model], model_f);
733  for (unsigned model_surface = 0; model_surface < _f[model]->numberSurfaces(); ++model_surface)
734  {
735  if (model_f[model_surface] > _f[model]->_f_tol)
736  // bummer, this model is not admissible at the returned_stress
737  all_admissible = false;
738  yf_at_returned_stress.push_back(model_f[model_surface]);
739  }
740  if (!all_admissible)
741  // no point in continuing computing yield functions
742  break;
743  }
744 
745  if (!all_admissible)
746  {
747  // we tried using the returned value of stress predicted by
748  // the_single_plastic_model, but it wasn't admissible according
749  // to other plastic models. We need to set:
750  stress = trial_stress;
751  for (unsigned model = 0; model < _num_models; ++model)
752  intnl[model] = intnl_old[model];
753  // and calculate the remainder of the yield functions at trial_stress
754  for (unsigned model = the_single_plastic_model; model < _num_models; ++model)
755  {
756  _f[model]->yieldFunctionV(trial_stress, intnl[model], model_f);
757  for (unsigned model_surface = 0; model_surface < _f[model]->numberSurfaces(); ++model_surface)
758  yf.push_back(model_f[model_surface]);
759  }
760  num_successful_plastic_returns = 0;
761  return false;
762  }
763 
764  // So the customized returnMap algorithm can provide a returned
765  // (stress, intnl) configuration, and that is admissible according
766  // to all plastic models
767  yf.resize(0);
768  for (unsigned surface = 0; surface < yf_at_returned_stress.size(); ++surface)
769  yf.push_back(yf_at_returned_stress[surface]);
770  delta_dp = model_delta_dp;
771  for (unsigned model_surface = 0; model_surface < _f[the_single_plastic_model]->numberSurfaces();
772  ++model_surface)
773  {
774  cumulative_pm[_surfaces_given_model[the_single_plastic_model][model_surface]] +=
775  model_pm[model_surface];
776  pm[_surfaces_given_model[the_single_plastic_model][model_surface]] = model_pm[model_surface];
777  }
778  custom_model = the_single_plastic_model;
779  return true;
780 }
std::vector< std::vector< unsigned int > > _surfaces_given_model
_surfaces_given_model[model_number] = vector of surface numbers for this model
unsigned int _num_surfaces
Number of surfaces within the plastic models.
void assign(const TypeTensor< T2 > &)
const PertinentGeochemicalSystem model(database, {"H2O", "H+", "HCO3-", "O2(aq)", "Ca++", ">(s)FeOH", "radius_neg1", "radius_neg1.5"}, {"Calcite"}, {}, {"Calcite_asdf"}, {"CH4(aq)"}, {">(s)FeOCa+"}, "O2(aq)", "e-")
unsigned int _num_models
Number of plastic models for this material.
std::vector< const SolidMechanicsPlasticModel * > _f
User objects that define the yield functions, flow potentials, etc.

◆ validParams()

InputParameters MultiPlasticityDebugger::validParams ( )
static

Definition at line 15 of file MultiPlasticityDebugger.C.

Referenced by ComputeMultiPlasticityStress::validParams().

16 {
18  MooseEnum debug_fspb_type("none crash jacobian jacobian_and_linear_system", "none");
19  params.addParam<MooseEnum>("debug_fspb",
20  debug_fspb_type,
21  "Debug types for use by developers when creating new "
22  "plasticity models, not for general use. 2 = debug Jacobian "
23  "entries, 3 = check the entire Jacobian, and check Ax=b");
24  params.addParam<RealTensorValue>("debug_jac_at_stress",
26  "Debug Jacobian entries at this stress. For use by developers");
27  params.addParam<std::vector<Real>>("debug_jac_at_pm",
28  "Debug Jacobian entries at these plastic multipliers");
29  params.addParam<std::vector<Real>>("debug_jac_at_intnl",
30  "Debug Jacobian entries at these internal parameters");
31  params.addParam<Real>(
32  "debug_stress_change", 1.0, "Debug finite differencing parameter for the stress");
33  params.addParam<std::vector<Real>>(
34  "debug_pm_change", "Debug finite differencing parameters for the plastic multipliers");
35  params.addParam<std::vector<Real>>(
36  "debug_intnl_change", "Debug finite differencing parameters for the internal parameters");
37  return params;
38 }
void addParam(const std::string &name, const std::initializer_list< typename T::value_type > &value, const std::string &doc_string)
TensorValue< Real > RealTensorValue
static InputParameters validParams()
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real

◆ yieldFunction()

void MultiPlasticityRawComponentAssembler::yieldFunction ( const RankTwoTensor stress,
const std::vector< Real > &  intnl,
const std::vector< bool > &  active,
std::vector< Real > &  f 
)
protectedvirtualinherited

The active yield function(s)

Parameters
stressthe stress at which to calculate the yield function
intnlvector of internal parameters
activeset of active constraints - only the active yield functions are put into "f"
[out]fthe yield function (or functions in the case of multisurface plasticity)

Definition at line 96 of file MultiPlasticityRawComponentAssembler.C.

Referenced by ComputeMultiPlasticityStress::buildDumbOrder(), MultiPlasticityLinearSystem::calculateConstraints(), ComputeMultiPlasticityStress::checkAdmissible(), fddyieldFunction_dintnl(), fddyieldFunction_dstress(), and ComputeMultiPlasticityStress::returnMap().

100 {
101  mooseAssert(intnl.size() == _num_models, "Incorrect size of internal parameters");
102  mooseAssert(active.size() == _num_surfaces, "Incorrect size of active");
103 
104  f.resize(0);
105  std::vector<unsigned int> active_surfaces_of_model;
106  std::vector<unsigned int>::iterator active_surface;
107  std::vector<Real> model_f;
108  for (unsigned model = 0; model < _num_models; ++model)
109  {
110  activeModelSurfaces(model, active, active_surfaces_of_model);
111  if (active_surfaces_of_model.size() > 0)
112  {
113  _f[model]->yieldFunctionV(stress, intnl[model], model_f);
114  for (active_surface = active_surfaces_of_model.begin();
115  active_surface != active_surfaces_of_model.end();
116  ++active_surface)
117  f.push_back(model_f[*active_surface]);
118  }
119  }
120 }
void activeModelSurfaces(int model, const std::vector< bool > &active, std::vector< unsigned int > &active_surfaces_of_model)
Returns the internal surface number(s) of the active surfaces of the given model This may be of size=...
unsigned int _num_surfaces
Number of surfaces within the plastic models.
Real f(Real x)
Test function for Brents method.
const PertinentGeochemicalSystem model(database, {"H2O", "H+", "HCO3-", "O2(aq)", "Ca++", ">(s)FeOH", "radius_neg1", "radius_neg1.5"}, {"Calcite"}, {}, {"Calcite_asdf"}, {"CH4(aq)"}, {">(s)FeOCa+"}, "O2(aq)", "e-")
unsigned int _num_models
Number of plastic models for this material.
std::vector< const SolidMechanicsPlasticModel * > _f
User objects that define the yield functions, flow potentials, etc.

Member Data Documentation

◆ _f

std::vector<const SolidMechanicsPlasticModel *> MultiPlasticityRawComponentAssembler::_f
protectedinherited

User objects that define the yield functions, flow potentials, etc.

Definition at line 70 of file MultiPlasticityRawComponentAssembler.h.

Referenced by MultiPlasticityRawComponentAssembler::activeModelSurfaces(), MultiPlasticityRawComponentAssembler::activeSurfaces(), MultiPlasticityRawComponentAssembler::anyActiveSurfaces(), ComputeMultiPlasticityStress::applyKuhnTucker(), MultiPlasticityRawComponentAssembler::buildActiveConstraints(), MultiPlasticityRawComponentAssembler::buildActiveConstraintsJoint(), MultiPlasticityRawComponentAssembler::buildActiveConstraintsRock(), ComputeMultiPlasticityStress::canAddConstraints(), ComputeMultiPlasticityStress::checkAdmissible(), ComputeMultiPlasticityStress::checkKuhnTucker(), MultiPlasticityRawComponentAssembler::dflowPotential_dintnl(), MultiPlasticityRawComponentAssembler::dflowPotential_dstress(), MultiPlasticityRawComponentAssembler::dhardPotential_dintnl(), MultiPlasticityRawComponentAssembler::dhardPotential_dstress(), MultiPlasticityRawComponentAssembler::dyieldFunction_dintnl(), MultiPlasticityRawComponentAssembler::dyieldFunction_dstress(), MultiPlasticityRawComponentAssembler::flowPotential(), MultiPlasticityRawComponentAssembler::hardPotential(), MultiPlasticityLinearSystem::MultiPlasticityLinearSystem(), MultiPlasticityRawComponentAssembler::MultiPlasticityRawComponentAssembler(), ComputeMultiPlasticityStress::quickStep(), ComputeMultiPlasticityStress::residual2(), ComputeMultiPlasticityStress::returnMap(), MultiPlasticityRawComponentAssembler::returnMapAll(), ComputeMultiPlasticityStress::singleStep(), and MultiPlasticityRawComponentAssembler::yieldFunction().

◆ _fspb_debug

MooseEnum MultiPlasticityDebugger::_fspb_debug
protected

none - don't do any debugging crash - currently inactive jacobian - check the jacobian entries jacobian_and_linear_system - check entire jacobian and check that Ax=b

Definition at line 57 of file MultiPlasticityDebugger.h.

Referenced by ComputeMultiPlasticityStress::computeQpStress(), and ComputeMultiPlasticityStress::initQpStatefulProperties().

◆ _fspb_debug_intnl

std::vector<Real> MultiPlasticityDebugger::_fspb_debug_intnl
protected

Debug the Jacobian entires at these internal parameters.

Definition at line 66 of file MultiPlasticityDebugger.h.

Referenced by checkDerivatives(), checkJacobian(), checkSolution(), outputAndCheckDebugParameters(), and ComputeMultiPlasticityStress::plasticStep().

◆ _fspb_debug_intnl_change

std::vector<Real> MultiPlasticityDebugger::_fspb_debug_intnl_change
protected

Debug finite-differencing parameters for the internal parameters.

Definition at line 75 of file MultiPlasticityDebugger.h.

Referenced by fddflowPotential_dintnl(), fddyieldFunction_dintnl(), fdJacobian(), and outputAndCheckDebugParameters().

◆ _fspb_debug_pm

std::vector<Real> MultiPlasticityDebugger::_fspb_debug_pm
protected

Debug the Jacobian entires at these plastic multipliers.

Definition at line 63 of file MultiPlasticityDebugger.h.

Referenced by checkJacobian(), checkSolution(), outputAndCheckDebugParameters(), and ComputeMultiPlasticityStress::plasticStep().

◆ _fspb_debug_pm_change

std::vector<Real> MultiPlasticityDebugger::_fspb_debug_pm_change
protected

Debug finite-differencing parameters for the plastic multipliers.

Definition at line 72 of file MultiPlasticityDebugger.h.

Referenced by fdJacobian(), and outputAndCheckDebugParameters().

◆ _fspb_debug_stress

RankTwoTensor MultiPlasticityDebugger::_fspb_debug_stress
protected

Debug the Jacobian entries at this stress.

Definition at line 60 of file MultiPlasticityDebugger.h.

Referenced by checkDerivatives(), checkJacobian(), checkSolution(), outputAndCheckDebugParameters(), and ComputeMultiPlasticityStress::plasticStep().

◆ _fspb_debug_stress_change

Real MultiPlasticityDebugger::_fspb_debug_stress_change
protected

Debug finite-differencing parameter for the stress.

Definition at line 69 of file MultiPlasticityDebugger.h.

Referenced by fddflowPotential_dstress(), fddyieldFunction_dstress(), fdJacobian(), and outputAndCheckDebugParameters().

◆ _min_f_tol

Real MultiPlasticityLinearSystem::_min_f_tol
protectedinherited

Minimum value of the _f_tol parameters for the Yield Function User Objects.

Definition at line 131 of file MultiPlasticityLinearSystem.h.

Referenced by MultiPlasticityLinearSystem::eliminateLinearDependence(), and MultiPlasticityLinearSystem::MultiPlasticityLinearSystem().

◆ _num_models

unsigned int MultiPlasticityRawComponentAssembler::_num_models
protectedinherited

◆ _num_surfaces

unsigned int MultiPlasticityRawComponentAssembler::_num_surfaces
protectedinherited

Number of surfaces within the plastic models.

For many situations this will be = _num_models since each model will contain just one surface. More generally it is >= _num_models. For instance, Mohr-Coulomb is a single model with 6 surfaces

Definition at line 61 of file MultiPlasticityRawComponentAssembler.h.

Referenced by ComputeMultiPlasticityStress::activeCombinationNumber(), ComputeMultiPlasticityStress::applyKuhnTucker(), MultiPlasticityRawComponentAssembler::buildActiveConstraints(), ComputeMultiPlasticityStress::buildDumbOrder(), MultiPlasticityLinearSystem::calculateConstraints(), MultiPlasticityLinearSystem::calculateJacobian(), MultiPlasticityLinearSystem::calculateRHS(), ComputeMultiPlasticityStress::canAddConstraints(), ComputeMultiPlasticityStress::canIncrementDumb(), ComputeMultiPlasticityStress::changeScheme(), ComputeMultiPlasticityStress::checkAdmissible(), checkDerivatives(), checkJacobian(), ComputeMultiPlasticityStress::checkKuhnTucker(), checkSolution(), ComputeMultiPlasticityStress::ComputeMultiPlasticityStress(), ComputeMultiPlasticityStress::computeQpStress(), ComputeMultiPlasticityStress::consistentTangentOperator(), MultiPlasticityRawComponentAssembler::dflowPotential_dintnl(), MultiPlasticityRawComponentAssembler::dflowPotential_dstress(), MultiPlasticityRawComponentAssembler::dhardPotential_dintnl(), MultiPlasticityRawComponentAssembler::dhardPotential_dstress(), dof_included(), MultiPlasticityRawComponentAssembler::dyieldFunction_dintnl(), MultiPlasticityRawComponentAssembler::dyieldFunction_dstress(), MultiPlasticityLinearSystem::eliminateLinearDependence(), fddflowPotential_dintnl(), fddflowPotential_dstress(), fddyieldFunction_dintnl(), fddyieldFunction_dstress(), fdJacobian(), MultiPlasticityRawComponentAssembler::flowPotential(), MultiPlasticityRawComponentAssembler::hardPotential(), ComputeMultiPlasticityStress::incrementDumb(), ComputeMultiPlasticityStress::initQpStatefulProperties(), MultiPlasticityRawComponentAssembler::MultiPlasticityRawComponentAssembler(), MultiPlasticityLinearSystem::nrStep(), ComputeMultiPlasticityStress::numberActive(), outputAndCheckDebugParameters(), ComputeMultiPlasticityStress::plasticStep(), ComputeMultiPlasticityStress::reinstateLinearDependentConstraints(), ComputeMultiPlasticityStress::residual2(), ComputeMultiPlasticityStress::returnMap(), MultiPlasticityRawComponentAssembler::returnMapAll(), ComputeMultiPlasticityStress::singleStep(), and MultiPlasticityRawComponentAssembler::yieldFunction().

◆ _params

const InputParameters& MultiPlasticityRawComponentAssembler::_params
protectedinherited

◆ _specialIC

MooseEnum MultiPlasticityRawComponentAssembler::_specialIC
protectedinherited

◆ _surfaces_given_model

std::vector<std::vector<unsigned int> > MultiPlasticityRawComponentAssembler::_surfaces_given_model
protectedinherited

◆ _svd_tol

Real MultiPlasticityLinearSystem::_svd_tol
protectedinherited

Tolerance on the minimum ratio of singular values before flow-directions are deemed linearly dependent.

Definition at line 128 of file MultiPlasticityLinearSystem.h.

Referenced by MultiPlasticityLinearSystem::eliminateLinearDependence().


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