Line data Source code
1 : /****************************************************************/
2 : /* DO NOT MODIFY THIS HEADER */
3 : /* BlackBear */
4 : /* */
5 : /* (c) 2017 Battelle Energy Alliance, LLC */
6 : /* ALL RIGHTS RESERVED */
7 : /* */
8 : /* Prepared by Battelle Energy Alliance, LLC */
9 : /* Under Contract No. DE-AC07-05ID14517 */
10 : /* With the U. S. Department of Energy */
11 : /* */
12 : /* See COPYRIGHT for full restrictions */
13 : /****************************************************************/
14 :
15 : #ifdef NEML_ENABLED
16 :
17 : #include "CauchyStressFromNEML.h"
18 :
19 : #include "NEMLStressBase.h"
20 :
21 : registerMooseObject("BlackBearApp", CauchyStressFromNEML);
22 :
23 : InputParameters
24 1174 : CauchyStressFromNEML::validParams()
25 : {
26 1174 : InputParameters params = ComputeLagrangianStressCauchy::validParams();
27 :
28 2348 : params.addRequiredParam<FileName>("database", "Path to NEML XML database.");
29 2348 : params.addRequiredParam<std::string>("model", "Model name in NEML database.");
30 2348 : params.addCoupledVar("temperature", 0.0, "Coupled temperature");
31 :
32 1174 : return params;
33 0 : }
34 :
35 900 : CauchyStressFromNEML::CauchyStressFromNEML(const InputParameters & parameters)
36 : : ComputeLagrangianStressCauchy(parameters),
37 1800 : _fname(getParam<FileName>("database")),
38 1800 : _mname(getParam<std::string>("model")),
39 900 : _temperature(coupledValue("temperature")),
40 900 : _temperature_old(coupledValueOld("temperature")),
41 900 : _history(declareProperty<std::vector<Real>>(_base_name + "history")),
42 1800 : _history_old(getMaterialPropertyOld<std::vector<Real>>(_base_name + "history")),
43 900 : _energy(declareProperty<Real>(_base_name + "energy")),
44 1800 : _energy_old(getMaterialPropertyOld<Real>(_base_name + "energy")),
45 900 : _dissipation(declareProperty<Real>(_base_name + "dissipation")),
46 900 : _dissipation_old(declareProperty<Real>(_base_name + "dissipation_old")),
47 900 : _linear_rotation(declareProperty<RankTwoTensor>(_base_name + "linear_rotation")),
48 1800 : _linear_rotation_old(getMaterialPropertyOld<RankTwoTensor>(_base_name + "linear_rotation")),
49 1800 : _cauchy_stress_old(getMaterialPropertyOld<RankTwoTensor>(_base_name + "cauchy_stress")),
50 1800 : _mechanical_strain(getMaterialProperty<RankTwoTensor>(_base_name + "mechanical_strain")),
51 1800 : _mechanical_strain_old(getMaterialPropertyOld<RankTwoTensor>(_base_name + "mechanical_strain")),
52 900 : _inelastic_strain(declareProperty<RankTwoTensor>(_base_name + "inelastic_strain")),
53 900 : _elastic_strain(declareProperty<RankTwoTensor>(_base_name + "elastic_strain")),
54 2700 : _dissipation_rate(declareProperty<Real>(_base_name + "dissipation_rate"))
55 : {
56 : // Check that the file is readable
57 900 : MooseUtils::checkFileReadable(_fname);
58 :
59 : // Will throw an exception if it doesn't succeed
60 : try
61 : {
62 2700 : _model = neml::parse_xml_unique(_fname, _mname);
63 : }
64 0 : catch (const neml::NEMLError & e)
65 : {
66 0 : paramError("Unable to load NEML model " + _mname + " from file " + _fname);
67 0 : }
68 900 : }
69 :
70 : void
71 280 : CauchyStressFromNEML::reset_state(const std::vector<unsigned int> & indices, unsigned int qp)
72 : {
73 : // Guard from zero, just for MOOSE...
74 280 : if (_model->nstore() == 0)
75 0 : return;
76 :
77 : // Form a NEML history object with the initial state
78 280 : std::vector<Real> init(_model->nstore(), 0);
79 280 : _model->init_store(&init.front());
80 :
81 : // Reset!
82 560 : for (auto i : indices)
83 280 : _history[qp][i] = init[i];
84 : }
85 :
86 : std::vector<unsigned int>
87 14 : CauchyStressFromNEML::provide_indices(const std::vector<std::string> & to_reset)
88 : {
89 : std::vector<unsigned int> indices;
90 :
91 : // Grab the state names...
92 14 : auto names = _model->report_internal_variable_names();
93 :
94 : // Iterate through names resetting each if found,
95 : // raise an error if you don't find it
96 28 : for (auto name : to_reset)
97 : {
98 14 : auto loc = std::find(names.begin(), names.end(), name);
99 14 : if (loc == names.end())
100 0 : mooseError("One of the state variables in the list "
101 : "requested for resetting does not exist "
102 : "in the NEML material model");
103 14 : unsigned int i = loc - names.begin();
104 14 : indices.push_back(i);
105 : }
106 :
107 14 : return indices;
108 14 : }
109 :
110 : void
111 4041681 : CauchyStressFromNEML::computeQpCauchyStress()
112 : {
113 : // Setup all the Mandel notation things we need
114 : double s_np1[6];
115 : double s_n[6];
116 4041681 : NEMLStressBase::rankTwoTensorToNeml(_cauchy_stress_old[_qp], s_n);
117 :
118 : // Strain
119 : double e_np1[6];
120 4041681 : NEMLStressBase::rankTwoTensorToNeml(_mechanical_strain[_qp], e_np1);
121 : double e_n[6];
122 4041681 : NEMLStressBase::rankTwoTensorToNeml(_mechanical_strain_old[_qp], e_n);
123 :
124 : // Vorticity
125 4041681 : RankTwoTensor L;
126 4041681 : if (_large_kinematics)
127 : {
128 1461176 : L = RankTwoTensor::Identity() - _inv_df[_qp];
129 : }
130 : else
131 : {
132 : L.zero();
133 : }
134 4041681 : _linear_rotation[_qp] = _linear_rotation_old[_qp] + (L - L.transpose()) / 2.0;
135 :
136 : double w_np1[3];
137 4041681 : rankTwoSkewToNEML(_linear_rotation[_qp], w_np1);
138 : double w_n[3];
139 4041681 : rankTwoSkewToNEML(_linear_rotation_old[_qp], w_n);
140 :
141 : // Time
142 4041681 : double t_np1 = _t;
143 4041681 : double t_n = _t - _dt;
144 :
145 : // Temperature
146 4041681 : double T_np1 = _temperature[_qp];
147 4041681 : double T_n = _temperature_old[_qp];
148 :
149 : // Internal state
150 : double * h_np1;
151 : const double * h_n;
152 :
153 : // Just to keep MOOSE debug happy
154 4041681 : if (_model->nstore() > 0)
155 : {
156 4041681 : h_np1 = &(_history[_qp][0]);
157 4041681 : h_n = &(_history_old[_qp][0]);
158 : }
159 : else
160 : {
161 : h_np1 = nullptr;
162 : h_n = nullptr;
163 : }
164 :
165 : // Energy
166 : double u_np1;
167 4041681 : double u_n = _energy_old[_qp];
168 :
169 : // Dissipation
170 : double p_np1;
171 4041681 : double p_n = _dissipation_old[_qp];
172 :
173 : // Tangent
174 : double A_np1[36];
175 : double B_np1[18];
176 :
177 : try
178 : {
179 : // Call NEML!
180 4041681 : if (_large_kinematics)
181 : {
182 1461176 : _model->update_ld_inc(e_np1,
183 : e_n,
184 : w_np1,
185 : w_n,
186 : T_np1,
187 : T_n,
188 : t_np1,
189 : t_n,
190 : s_np1,
191 : s_n,
192 : h_np1,
193 : h_n,
194 : A_np1,
195 : B_np1,
196 : u_np1,
197 : u_n,
198 : p_np1,
199 : p_n);
200 : }
201 : else
202 : {
203 2580505 : _model->update_sd(e_np1,
204 : e_n,
205 : T_np1,
206 : T_n,
207 : t_np1,
208 : t_n,
209 : s_np1,
210 : s_n,
211 : h_np1,
212 : h_n,
213 : A_np1,
214 : u_np1,
215 : u_n,
216 : p_np1,
217 : p_n);
218 : std::fill(B_np1, B_np1 + 18, 0.0);
219 : }
220 :
221 : double estrain[6];
222 4041636 : _model->elastic_strains(s_np1, T_np1, h_np1, estrain);
223 :
224 : // Translate back from Mandel notation
225 4041636 : NEMLStressBase::nemlToRankTwoTensor(s_np1, _cauchy_stress[_qp]);
226 4041636 : recombineTangentNEML(A_np1, B_np1, _cauchy_jacobian[_qp]);
227 4041636 : _energy[_qp] = u_np1;
228 4041636 : _dissipation[_qp] = p_np1;
229 4041636 : _dissipation_rate[_qp] = (p_np1 - p_n) / (t_np1 - t_n);
230 :
231 4041636 : NEMLStressBase::nemlToRankTwoTensor(estrain, _elastic_strain[_qp]);
232 4041636 : _inelastic_strain[_qp] = _mechanical_strain[_qp] - _elastic_strain[_qp];
233 : }
234 45 : catch (const neml::NEMLError & e)
235 : {
236 90 : throw MooseException("NEML error: " + e.message());
237 45 : }
238 4041636 : }
239 :
240 : void
241 20984 : CauchyStressFromNEML::initQpStatefulProperties()
242 : {
243 20984 : ComputeLagrangianStressCauchy::initQpStatefulProperties();
244 :
245 20984 : _history[_qp].resize(_model->nstore());
246 : try
247 : {
248 : // This is only needed because MOOSE whines about zero sized vectors
249 : // that are not initialized
250 20984 : if (_history[_qp].size() > 0)
251 20984 : _model->init_store(&_history[_qp].front());
252 : }
253 0 : catch (const neml::NEMLError & e)
254 : {
255 0 : throw MooseException("NEML error: " + e.message());
256 0 : }
257 :
258 20984 : _linear_rotation[_qp].zero();
259 :
260 20984 : _energy[_qp] = 0.0;
261 20984 : _dissipation[_qp] = 0.0;
262 20984 : _dissipation_rate[_qp] = 0.0;
263 20984 : }
264 :
265 : void
266 8083362 : rankTwoSkewToNEML(const RankTwoTensor & in, double * const out)
267 : {
268 8083362 : out[0] = -in(1, 2);
269 8083362 : out[1] = in(0, 2);
270 8083362 : out[2] = -in(0, 1);
271 8083362 : }
272 :
273 : void
274 4041636 : recombineTangentNEML(const double * const Dpart, const double * const Wpart, RankFourTensor & out)
275 : {
276 4041636 : std::vector<double> data(81);
277 4041636 : neml::transform_fourth(Dpart, Wpart, &data[0]);
278 4041636 : out.fillFromInputVector(data, RankFourTensor::FillMethod::general);
279 4041636 : }
280 :
281 : #endif // NEML_ENABLED
|
