From 6bc930ca3bf04054e3bfe3d801f5d7f25e6118da Mon Sep 17 00:00:00 2001 From: Thomas Poulet Date: Thu, 10 Jun 2021 11:47:00 +0800 Subject: [PATCH] #180, renamed RedbackFaultAnisotropicPermAux in more generic RedbackAnisotropicTensorAux Now using functions for longitudinal and transversal components. Those function can take the potential field as parameter --- .../auxkernels/RedbackAnisotropicTensorAux.h | 54 +++++++ .../RedbackFaultAnisotropicPermAux.h | 38 ----- src/auxkernels/RedbackAnisotropicTensorAux.C | 153 ++++++++++++++++++ .../RedbackFaultAnisotropicPermAux.C | 85 ---------- 4 files changed, 207 insertions(+), 123 deletions(-) create mode 100644 include/auxkernels/RedbackAnisotropicTensorAux.h delete mode 100644 include/auxkernels/RedbackFaultAnisotropicPermAux.h create mode 100644 src/auxkernels/RedbackAnisotropicTensorAux.C delete mode 100644 src/auxkernels/RedbackFaultAnisotropicPermAux.C diff --git a/include/auxkernels/RedbackAnisotropicTensorAux.h b/include/auxkernels/RedbackAnisotropicTensorAux.h new file mode 100644 index 00000000..e43534f1 --- /dev/null +++ b/include/auxkernels/RedbackAnisotropicTensorAux.h @@ -0,0 +1,54 @@ +/****************************************************************/ +/* DO NOT MODIFY THIS HEADER */ +/* REDBACK - Rock mEchanics with Dissipative feedBACKs */ +/* */ +/* (c) 2014 CSIRO */ +/* ALL RIGHTS RESERVED */ +/* */ +/* Prepared by CSIRO and UNSW Australia */ +/* */ +/* See COPYRIGHT for full restrictions */ +/****************************************************************/ + +#ifndef REDBACKANISOTROPICTENSORAUX +#define REDBACKANISOTROPICTENSORAUX + +#include "AuxKernel.h" +#include "RankTwoTensor.h" +#include "FunctionParserUtils.h" + +/** + * AuxKernel to compute anisotropic tensor component, defined by + * its longitudinal and transversal permeability components, as well as + * a directionality provided by a normal (potential) field. The longitudinal + * direction follows the isopotential surfaces. + * (Used in practice to set permeability tensor.) + */ + +// Forward declarations +class RedbackAnisotropicTensorAux; + +template <> +InputParameters validParams(); + +class RedbackAnisotropicTensorAux : public AuxKernel, public FunctionParserUtils +{ +public: + RedbackAnisotropicTensorAux(const InputParameters & parameters); + virtual ~RedbackAnisotropicTensorAux() {} + +protected: + virtual Real computeValue(); + const VariableValue & _potential; + const VariableGradient & _n; // normal vector (non-normalised) to iso-potentials + std::string _function_long; // function expression for longitudinal component + std::string _function_tran; // function expression for longitudinal component + SymFunctionPtr _func_F_long; // function parser object + SymFunctionPtr _func_F_tran; // function parser object + usingFunctionParserUtilsMembers(false); + + int _i_row, _i_column; // indices i,j of tensor to set + Real _k_long, _k_tran; // values for longitudinal and transversal components (computed from functions) +}; + +#endif // REDBACKANISOTROPICTENSORAUX diff --git a/include/auxkernels/RedbackFaultAnisotropicPermAux.h b/include/auxkernels/RedbackFaultAnisotropicPermAux.h deleted file mode 100644 index 56a3888a..00000000 --- a/include/auxkernels/RedbackFaultAnisotropicPermAux.h +++ /dev/null @@ -1,38 +0,0 @@ -/****************************************************************/ -/* DO NOT MODIFY THIS HEADER */ -/* REDBACK - Rock mEchanics with Dissipative feedBACKs */ -/* */ -/* (c) 2014 CSIRO */ -/* ALL RIGHTS RESERVED */ -/* */ -/* Prepared by CSIRO and UNSW Australia */ -/* */ -/* See COPYRIGHT for full restrictions */ -/****************************************************************/ - -#ifndef REDBACKFAULTANISOTROPICPERMAUX_H -#define REDBACKFAULTANISOTROPICPERMAUX_H - -#include "AuxKernel.h" -#include "RankTwoTensor.h" - -// Forward declarations -class RedbackFaultAnisotropicPermAux; - -template <> -InputParameters validParams(); - -class RedbackFaultAnisotropicPermAux : public AuxKernel -{ -public: - RedbackFaultAnisotropicPermAux(const InputParameters & parameters); - virtual ~RedbackFaultAnisotropicPermAux() {} - -protected: - virtual Real computeValue(); - const VariableGradient & _n; // normal vector (non-normalised) - Real _k_long, _k_tran; // perm values for longitudinal and transversal components - int _i_row, _i_column; // indices i,j of perm tensor to set -}; - -#endif // REDBACKFAULTANISOTROPICPERMAUX_H diff --git a/src/auxkernels/RedbackAnisotropicTensorAux.C b/src/auxkernels/RedbackAnisotropicTensorAux.C new file mode 100644 index 00000000..8a80846d --- /dev/null +++ b/src/auxkernels/RedbackAnisotropicTensorAux.C @@ -0,0 +1,153 @@ +/****************************************************************/ +/* DO NOT MODIFY THIS HEADER */ +/* REDBACK - Rock mEchanics with Dissipative feedBACKs */ +/* */ +/* (c) 2014 CSIRO */ +/* ALL RIGHTS RESERVED */ +/* */ +/* Prepared by CSIRO and UNSW Australia */ +/* */ +/* See COPYRIGHT for full restrictions */ +/****************************************************************/ + +// AuxKernel to compute anisotropic tensor component, defined by +// its longitudinal and transversal components, as well as +// a directionality provided by a normal (potential) field. The longitudinal +// direction follows the isopotential surfaces. +// (Used in practice to set permeability tensor.) +#include "RedbackAnisotropicTensorAux.h" + +registerMooseObject("RedbackApp", RedbackAnisotropicTensorAux); + +template <> +InputParameters +validParams() +{ + MooseEnum component("x=0 y=1 z=2"); + InputParameters params = AuxKernel::validParams(); + params += FunctionParserUtils::validParams(); + params.addClassDescription("AuxKernel to set anisotropic tensor from " + "longitudinal and transversal components."); + + params.addRequiredCoupledVar("potential_field", "Scalar field from which " + "to infer directions in that block. " + "(longitudinal = along iso-potentials; transversal = normal to iso-potentials)"); + params.addRequiredCustomTypeParam( + "function_long", "FunctionExpression", "function expression defining " + "longitudinal component of tensor. It can use the potential field as argument."); + params.addRequiredCustomTypeParam( + "function_tran", "FunctionExpression", "function expression defining " + "transversal component of tensor. It can use the potential field as argument."); + params.addParam>( + "constant_names", "Vector of constants used in the parsed functions."); + params.addParam>( + "constant_expressions", + "Vector of values for the constants in constant_names (can be an FParser expression)"); + params.addParam("row", component, "The tensor component (row) to assign: x, y or z"); + params.addParam("column", component, "The tensor component (column) to assign: x, y or z"); + + return params; +} + +RedbackAnisotropicTensorAux::RedbackAnisotropicTensorAux(const InputParameters & parameters) + : AuxKernel(parameters), + FunctionParserUtils(parameters), + _potential(coupledValue("potential_field")), + _n(coupledGradient("potential_field")), + _function_long(getParam("function_long")), + _function_tran(getParam("function_tran")), + _i_row(getParam("row")), + _i_column(getParam("column")) +{ + // build variables argument (single variable representing the potential field) + std::string variables; + variables += getVar("potential_field", 0)->name(); + + // base function objects + _func_F_long = std::make_shared(); + _func_F_tran = std::make_shared(); + + // set FParser internal feature flags + setParserFeatureFlags(_func_F_long); + setParserFeatureFlags(_func_F_tran); + + // add the constant expressions + addFParserConstants(_func_F_long, + getParam>("constant_names"), + getParam>("constant_expressions")); + addFParserConstants(_func_F_tran, + getParam>("constant_names"), + getParam>("constant_expressions")); + + // parse functions + if (_func_F_long->Parse(_function_long, variables) >= 0) + mooseError( + "Invalid function\n", _function_long, "\nin RedbackAnisotropicTensorAux ", name(), ".\n", _func_F_long->ErrorMsg()); + if (_func_F_tran->Parse(_function_tran, variables) >= 0) + mooseError( + "Invalid function\n", _function_tran, "\nin RedbackAnisotropicTensorAux ", name(), ".\n", _func_F_tran->ErrorMsg()); + + // optimize + if (!_disable_fpoptimizer) + { + _func_F_long->Optimize(); + _func_F_tran->Optimize(); + } + + // just-in-time compile + if (_enable_jit) + { + _func_F_long->JITCompile(); + _func_F_tran->JITCompile(); + } + + // reserve storage for parameter passing bufefr + _func_params.resize(1); // single argument for the potential field +} + +Real +RedbackAnisotropicTensorAux::computeValue() +{ + RealTensorValue my_tensor, transfer_matrix; + // Set first column of transfer matrix as normalised normal vector + RealVectorValue u(_n[_qp]);//(0), _n[_qp](1), _n[_qp](2)); + u /= u.norm(); + // Second column vector is orthogonal to first one + RealVectorValue v(0, 0, 0); + if (_n[_qp](2) == 0) { + // normal vector in horizontal plane + if (_n[_qp](1) == 0) { + v(1) = 1; + // keep other 2 components zero + } else { + v(0) = _n[_qp](1); + v(1) = -_n[_qp](0); + // keep third component zero + } + } else { + // normal vector has non-zero vertical component + // keep first component zero + v(1) = _n[_qp](2); + v(2) = -_n[_qp](1); + } + v /= v.norm(); + // Third vector as cross product to get orthonormal basis + transfer_matrix = RealTensorValue(u, v, u.cross(v)).transpose(); + + // Set function parameter for potential field and evaluate both functions + _func_params[0] = _potential[_qp]; + _k_long = evaluate(_func_F_long); + _k_tran = evaluate(_func_F_tran); + + // Define tensor in rotated referential (aligned with iso-potentials) + RealTensorValue k_prime; + k_prime.zero(); + k_prime(0, 0) = _k_tran; + k_prime(1, 1) = _k_long; + k_prime(2, 2) = _k_long; + + // Get tensor in original referential + my_tensor = transfer_matrix * (k_prime * transfer_matrix.transpose()); + + return my_tensor(_i_row, _i_column); +} diff --git a/src/auxkernels/RedbackFaultAnisotropicPermAux.C b/src/auxkernels/RedbackFaultAnisotropicPermAux.C deleted file mode 100644 index b29357e7..00000000 --- a/src/auxkernels/RedbackFaultAnisotropicPermAux.C +++ /dev/null @@ -1,85 +0,0 @@ -/****************************************************************/ -/* DO NOT MODIFY THIS HEADER */ -/* REDBACK - Rock mEchanics with Dissipative feedBACKs */ -/* */ -/* (c) 2014 CSIRO */ -/* ALL RIGHTS RESERVED */ -/* */ -/* Prepared by CSIRO and UNSW Australia */ -/* */ -/* See COPYRIGHT for full restrictions */ -/****************************************************************/ - -// AuxKernel to compute anistropic permeability component around a fault -// given by its (non-normalised) normal vectors, as well as -// longitudinal and transversal permeability components to assign -#include "RedbackFaultAnisotropicPermAux.h" - -registerMooseObject("RedbackApp", RedbackFaultAnisotropicPermAux); - -template <> -InputParameters -validParams() -{ - MooseEnum component("x=0 y=1 z=2"); - InputParameters params = validParams(); - params.addRequiredCoupledVar("fault_normal_field", "Scalar field from which to infer vector normal to fault in that block"); - params.addParam("k_long", 1.0, "Permeability, longitudinal component"); - params.addParam("k_tran", 1.0, "Permeability, transversal component"); - params.addParam("row", component, "The permeability tensor component (row) to assign"); - params.addParam("column", component, "The permeability tensor component (column) to assign"); - - return params; -} - -RedbackFaultAnisotropicPermAux::RedbackFaultAnisotropicPermAux(const InputParameters & parameters) - : AuxKernel(parameters), - _n(coupledGradient("fault_normal_field")), - _k_long(getParam("k_long")), - _k_tran(getParam("k_tran")), - _i_row(getParam("row")), - _i_column(getParam("column")) -{ -} - -Real -RedbackFaultAnisotropicPermAux::computeValue() -{ - RealTensorValue perm_tensor, transfer_matrix; - // Set first column of transfer matrix as normalised normal vector - RealVectorValue u(_n[_qp]);//(0), _n[_qp](1), _n[_qp](2)); - u /= u.norm(); - // Second column vector is orthogonal to first one - RealVectorValue v(0, 0, 0); - if (_n[_qp](2) == 0) { - // normal vector in horizontal plane - if (_n[_qp](1) == 0) { - v(1) = 1; - // keep other 2 components zero - } else { - v(0) = _n[_qp](1); - v(1) = -_n[_qp](0); - // keep third component zero - } - } else { - // normal vector has non-zero vertical component - // keep first component zero - v(1) = _n[_qp](2); - v(2) = -_n[_qp](1); - } - v /= v.norm(); - // Third vector as cross product to get orthonormal basis - transfer_matrix = RealTensorValue(u, v, u.cross(v)).transpose(); - - // Define permeability tensor in rotated (fault-aligned) referential - RealTensorValue k_prime; - k_prime.zero(); - k_prime(0, 0) = _k_tran; - k_prime(1, 1) = _k_long; - k_prime(2, 2) = _k_long; - - // Get perm tensor in original referential - perm_tensor = transfer_matrix * (k_prime * transfer_matrix.transpose()); - - return perm_tensor(_i_row, _i_column); -}