Line data Source code
1 : //* This file is part of the MOOSE framework
2 : //* https://mooseframework.inl.gov
3 : //*
4 : //* All rights reserved, see COPYRIGHT for full restrictions
5 : //* https://github.com/idaholab/moose/blob/master/COPYRIGHT
6 : //*
7 : //* Licensed under LGPL 2.1, please see LICENSE for details
8 : //* https://www.gnu.org/licenses/lgpl-2.1.html
9 :
10 : #include "INSADMaterial.h"
11 : #include "Function.h"
12 : #include "Assembly.h"
13 : #include "INSADObjectTracker.h"
14 : #include "FEProblemBase.h"
15 : #include "NonlinearSystemBase.h"
16 : #include "NS.h"
17 : #include "NavierStokesMethods.h"
18 :
19 : registerMooseObject("NavierStokesApp", INSADMaterial);
20 :
21 : InputParameters
22 2978 : INSADMaterial::validParams()
23 : {
24 2978 : InputParameters params = Material::validParams();
25 2978 : params.addClassDescription("This is the material class used to compute some of the strong "
26 : "residuals for the INS equations.");
27 5956 : params.addRequiredCoupledVar("velocity", "The velocity");
28 2978 : params.addRequiredCoupledVar(NS::pressure, "The pressure");
29 5956 : params.addParam<MaterialPropertyName>("mu_name", "mu", "The name of the dynamic viscosity");
30 5956 : params.addParam<MaterialPropertyName>("rho_name", "rho", "The name of the density");
31 2978 : return params;
32 0 : }
33 :
34 2307 : INSADMaterial::INSADMaterial(const InputParameters & parameters)
35 : : Material(parameters),
36 2307 : _velocity(adCoupledVectorValue("velocity")),
37 2307 : _grad_velocity(adCoupledVectorGradient("velocity")),
38 2307 : _grad_p(adCoupledGradient(NS::pressure)),
39 4614 : _mu(getADMaterialProperty<Real>("mu_name")),
40 4614 : _rho(getADMaterialProperty<Real>("rho_name")),
41 2307 : _velocity_dot(nullptr),
42 2307 : _mass_strong_residual(declareADProperty<Real>("mass_strong_residual")),
43 2307 : _advective_strong_residual(declareADProperty<RealVectorValue>("advective_strong_residual")),
44 : // We have to declare the below strong residuals for integrity check purposes even though we may
45 : // not compute them. This may incur some unnecessary cost for a non-sparse derivative container
46 : // since when the properties are resized the entire non-sparse derivative containers will be
47 : // initialized to zero
48 2307 : _td_strong_residual(declareADProperty<RealVectorValue>("td_strong_residual")),
49 2307 : _gravity_strong_residual(declareADProperty<RealVectorValue>("gravity_strong_residual")),
50 2307 : _boussinesq_strong_residual(declareADProperty<RealVectorValue>("boussinesq_strong_residual")),
51 2307 : _coupled_force_strong_residual(
52 2307 : declareADProperty<RealVectorValue>("coupled_force_strong_residual")),
53 2307 : _mesh_velocity(declareADProperty<RealVectorValue>("mesh_velocity")),
54 2307 : _relative_velocity(declareADProperty<RealVectorValue>("relative_velocity")),
55 2307 : _advected_mesh_strong_residual(
56 2307 : declareADProperty<RealVectorValue>("advected_mesh_strong_residual")),
57 : // _mms_function_strong_residual(declareProperty<RealVectorValue>("mms_function_strong_residual")),
58 4614 : _use_displaced_mesh(getParam<bool>("use_displaced_mesh")),
59 2307 : _ad_q_point(_bnd ? _assembly.adQPointsFace() : _assembly.adQPoints()),
60 2307 : _rz_radial_coord(_mesh.getAxisymmetricRadialCoord()),
61 4794 : _rz_axial_coord(_rz_radial_coord == 0 ? 1 : 0)
62 : {
63 4614 : if (!_fe_problem.hasUserObject("ins_ad_object_tracker"))
64 : {
65 653 : InputParameters tracker_params = INSADObjectTracker::validParams();
66 653 : tracker_params.addPrivateParam(MooseBase::app_param, &_app);
67 :
68 1306 : _fe_problem.addUserObject("INSADObjectTracker", "ins_ad_object_tracker", tracker_params);
69 653 : }
70 :
71 : // Bypass the UserObjectInterface method because it requires a UserObjectName param which we
72 : // don't need
73 2307 : _object_tracker = &_fe_problem.getUserObject<INSADObjectTracker>("ins_ad_object_tracker");
74 2307 : const_cast<INSADObjectTracker *>(_object_tracker)->addBlockIDs(this->blockIDs());
75 2307 : }
76 :
77 : void
78 2307 : INSADMaterial::resolveOptionalProperties()
79 : {
80 2307 : Material::resolveOptionalProperties();
81 :
82 4644 : for (const auto sub_id : this->blockIDs())
83 2337 : if ((_object_tracker->get<bool>("has_boussinesq", sub_id)))
84 : {
85 564 : _boussinesq_alphas[sub_id] = &getPossiblyConstantGenericMaterialPropertyByName<Real, true>(
86 282 : _object_tracker->get<MaterialPropertyName>("alpha", sub_id), _material_data, 0);
87 564 : _ref_temps[sub_id] = &getPossiblyConstantGenericMaterialPropertyByName<Real, false>(
88 282 : _object_tracker->get<MaterialPropertyName>("ref_temp", sub_id), _material_data, 0);
89 : }
90 2307 : }
91 :
92 : void
93 371745 : INSADMaterial::subdomainSetup()
94 : {
95 371745 : if ((_has_transient = _object_tracker->get<bool>("has_transient", _current_subdomain_id)))
96 103567 : _velocity_dot = &adCoupledVectorDot("velocity");
97 : else
98 268178 : _velocity_dot = nullptr;
99 :
100 371745 : if (auto alpha_it = _boussinesq_alphas.find(_current_subdomain_id);
101 : alpha_it != _boussinesq_alphas.end())
102 : {
103 191404 : _has_boussinesq = true;
104 191404 : _boussinesq_alpha = alpha_it->second;
105 191404 : auto & temp_var = _subproblem.getStandardVariable(
106 191404 : _tid, _object_tracker->get<std::string>("temperature", _current_subdomain_id));
107 191404 : addMooseVariableDependency(&temp_var);
108 191404 : _temperature = &temp_var.adSln();
109 191404 : _ref_temp = libmesh_map_find(_ref_temps, _current_subdomain_id);
110 : }
111 : else
112 : {
113 180341 : _has_boussinesq = false;
114 180341 : _boussinesq_alpha = nullptr;
115 180341 : _temperature = nullptr;
116 180341 : _ref_temp = nullptr;
117 : }
118 :
119 371745 : _has_gravity = _object_tracker->get<bool>("has_gravity", _current_subdomain_id);
120 371745 : if (_has_gravity || _has_boussinesq)
121 191608 : _gravity_vector = _object_tracker->get<RealVectorValue>("gravity", _current_subdomain_id);
122 : else
123 : _gravity_vector = 0;
124 :
125 : // Setup data for Arbitrary Lagrangian Eulerian (ALE) simulations in which the simulation domain
126 : // is displacing. We will need to subtract the mesh velocity from the velocity solution in order
127 : // to get the correct material velocity for the momentum convection term.
128 371745 : if ((_has_advected_mesh = _object_tracker->get<bool>("has_advected_mesh", _current_subdomain_id)))
129 : {
130 2542 : auto & disp_x = _subproblem.getStandardVariable(
131 2542 : _tid, _object_tracker->get<VariableName>("disp_x", _current_subdomain_id));
132 2542 : addMooseVariableDependency(&disp_x);
133 2542 : _disp_x_dot = &disp_x.adUDot();
134 2542 : _disp_x_sys_num = disp_x.sys().number();
135 2542 : _disp_x_num = (disp_x.kind() == Moose::VarKindType::VAR_SOLVER) &&
136 2542 : (_disp_x_sys_num == _fe_problem.currentNonlinearSystem().number())
137 2542 : ? disp_x.number()
138 : : libMesh::invalid_uint;
139 5084 : if (_object_tracker->isTrackerParamValid("disp_y", _current_subdomain_id))
140 : {
141 2542 : auto & disp_y = _subproblem.getStandardVariable(
142 2542 : _tid, _object_tracker->get<VariableName>("disp_y", _current_subdomain_id));
143 2542 : addMooseVariableDependency(&disp_y);
144 2542 : _disp_y_dot = &disp_y.adUDot();
145 2542 : _disp_y_sys_num = disp_y.sys().number();
146 2542 : _disp_y_num =
147 : disp_y.kind() == (Moose::VarKindType::VAR_SOLVER &&
148 : (_disp_y_sys_num == _fe_problem.currentNonlinearSystem().number()))
149 2542 : ? disp_y.number()
150 : : libMesh::invalid_uint;
151 : }
152 : else
153 : {
154 0 : _disp_y_dot = nullptr;
155 0 : _disp_y_sys_num = libMesh::invalid_uint;
156 0 : _disp_y_num = libMesh::invalid_uint;
157 : }
158 5084 : if (_object_tracker->isTrackerParamValid("disp_z", _current_subdomain_id))
159 : {
160 1244 : auto & disp_z = _subproblem.getStandardVariable(
161 1244 : _tid, _object_tracker->get<VariableName>("disp_z", _current_subdomain_id));
162 1244 : addMooseVariableDependency(&disp_z);
163 1244 : _disp_z_dot = &disp_z.adUDot();
164 1244 : _disp_z_sys_num = disp_z.sys().number();
165 1244 : _disp_z_num =
166 : disp_z.kind() == (Moose::VarKindType::VAR_SOLVER &&
167 : (_disp_z_sys_num == _fe_problem.currentNonlinearSystem().number()))
168 1244 : ? disp_z.number()
169 : : libMesh::invalid_uint;
170 : }
171 : else
172 : {
173 1298 : _disp_z_dot = nullptr;
174 1298 : _disp_z_sys_num = libMesh::invalid_uint;
175 1298 : _disp_z_num = libMesh::invalid_uint;
176 : }
177 : }
178 : else
179 369203 : _disp_x_dot = _disp_y_dot = _disp_z_dot = nullptr;
180 :
181 371745 : _viscous_form = static_cast<NS::ViscousForm>(
182 371745 : int(_object_tracker->get<MooseEnum>("viscous_form", _current_subdomain_id)));
183 :
184 371745 : if ((_has_coupled_force = _object_tracker->get<bool>("has_coupled_force", _current_subdomain_id)))
185 : {
186 37740 : _coupled_force_var.clear();
187 37740 : _coupled_force_vector_function.clear();
188 75480 : if (_object_tracker->isTrackerParamValid("coupled_force_var", _current_subdomain_id))
189 : {
190 38352 : const auto & var_names = _object_tracker->get<std::vector<VariableName>>(
191 19176 : "coupled_force_var", _current_subdomain_id);
192 38760 : for (const auto & var_name : var_names)
193 19584 : _coupled_force_var.push_back(&_subproblem.getVectorVariable(_tid, var_name).adSln());
194 : }
195 :
196 75480 : if (_object_tracker->isTrackerParamValid("coupled_force_vector_function",
197 37740 : _current_subdomain_id))
198 : {
199 37944 : const auto & func_names = _object_tracker->get<std::vector<FunctionName>>(
200 18972 : "coupled_force_vector_function", _current_subdomain_id);
201 38352 : for (const auto & func_name : func_names)
202 19380 : _coupled_force_vector_function.push_back(&_fe_problem.getFunction(func_name, _tid));
203 : }
204 : }
205 371745 : }
206 :
207 : void
208 24585542 : INSADMaterial::computeQpProperties()
209 : {
210 24585542 : _mass_strong_residual[_qp] = -NS::divergence(
211 49171084 : _grad_velocity[_qp], _velocity[_qp], _ad_q_point[_qp], _coord_sys, _rz_radial_coord);
212 :
213 49171084 : _advective_strong_residual[_qp] = _rho[_qp] * _grad_velocity[_qp] * _velocity[_qp];
214 24585542 : if (_has_transient)
215 9857884 : _td_strong_residual[_qp] = _rho[_qp] * (*_velocity_dot)[_qp];
216 24585542 : if (_has_gravity)
217 24525984 : _gravity_strong_residual[_qp] = -_rho[_qp] * _gravity_vector;
218 24585542 : if (_has_boussinesq)
219 24439968 : _boussinesq_strong_residual[_qp] = (*_boussinesq_alpha)[_qp] * _gravity_vector * _rho[_qp] *
220 24439968 : ((*_temperature)[_qp] - (*_ref_temp)[_qp]);
221 24585542 : _relative_velocity[_qp] = _velocity[_qp];
222 :
223 24585542 : if (_has_advected_mesh)
224 : {
225 463018 : _mesh_velocity[_qp](0) = (*_disp_x_dot)[_qp];
226 463018 : if (_disp_y_dot)
227 463018 : _mesh_velocity[_qp](1) = (*_disp_y_dot)[_qp];
228 463018 : if (_disp_z_dot)
229 188648 : _mesh_velocity[_qp](2) = (*_disp_z_dot)[_qp];
230 : _relative_velocity[_qp] -= _mesh_velocity[_qp];
231 1389054 : _advected_mesh_strong_residual[_qp] = -_rho[_qp] * _grad_velocity[_qp] * _mesh_velocity[_qp];
232 : }
233 :
234 24585542 : if (_has_coupled_force)
235 : {
236 1210688 : _coupled_force_strong_residual[_qp] = 0;
237 : mooseAssert(!(_coupled_force_var.empty() && _coupled_force_vector_function.empty()),
238 : "Either the coupled force var or the coupled force vector function must be "
239 : "non-empty in 'INSADMaterial'");
240 1966560 : for (const auto * var : _coupled_force_var)
241 : {
242 : mooseAssert(var, "null coupled variable in INSADMaterial");
243 755872 : _coupled_force_strong_residual[_qp] -= (*var)[_qp];
244 : }
245 1923552 : for (const auto * fn : _coupled_force_vector_function)
246 : {
247 : mooseAssert(fn, "null coupled function in INSADMaterial");
248 1425728 : _coupled_force_strong_residual[_qp] -= fn->vectorValue(_t, _q_point[_qp]);
249 : }
250 : }
251 :
252 : // // Future Addition
253 : // _mms_function_strong_residual[_qp] = -RealVectorValue(_x_vel_fn.value(_t, _q_point[_qp]),
254 : // _y_vel_fn.value(_t, _q_point[_qp]),
255 : // _z_vel_fn.value(_t, _q_point[_qp]));
256 24585542 : }
|
