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 6050 : INSADMaterial::validParams()
23 : {
24 6050 : InputParameters params = Material::validParams();
25 6050 : params.addClassDescription("This is the material class used to compute some of the strong "
26 : "residuals for the INS equations.");
27 12100 : params.addRequiredCoupledVar("velocity", "The velocity");
28 6050 : params.addRequiredCoupledVar(NS::pressure, "The pressure");
29 12100 : params.addParam<MaterialPropertyName>("mu_name", "mu", "The name of the dynamic viscosity");
30 12100 : params.addParam<MaterialPropertyName>("rho_name", "rho", "The name of the density");
31 6050 : return params;
32 0 : }
33 :
34 4731 : INSADMaterial::INSADMaterial(const InputParameters & parameters)
35 : : Material(parameters),
36 4731 : _velocity(adCoupledVectorValue("velocity")),
37 4731 : _grad_velocity(adCoupledVectorGradient("velocity")),
38 4731 : _grad_p(adCoupledGradient(NS::pressure)),
39 9462 : _mu(getADMaterialProperty<Real>("mu_name")),
40 9462 : _rho(getADMaterialProperty<Real>("rho_name")),
41 4731 : _velocity_dot(nullptr),
42 4731 : _mass_strong_residual(declareADProperty<Real>("mass_strong_residual")),
43 4731 : _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 4731 : _td_strong_residual(declareADProperty<RealVectorValue>("td_strong_residual")),
49 4731 : _gravity_strong_residual(declareADProperty<RealVectorValue>("gravity_strong_residual")),
50 4731 : _boussinesq_strong_residual(declareADProperty<RealVectorValue>("boussinesq_strong_residual")),
51 4731 : _coupled_force_strong_residual(
52 4731 : declareADProperty<RealVectorValue>("coupled_force_strong_residual")),
53 4731 : _mesh_velocity(declareADProperty<RealVectorValue>("mesh_velocity")),
54 4731 : _relative_velocity(declareADProperty<RealVectorValue>("relative_velocity")),
55 4731 : _advected_mesh_strong_residual(
56 4731 : declareADProperty<RealVectorValue>("advected_mesh_strong_residual")),
57 : // _mms_function_strong_residual(declareProperty<RealVectorValue>("mms_function_strong_residual")),
58 9462 : _use_displaced_mesh(getParam<bool>("use_displaced_mesh")),
59 4731 : _ad_q_point(_bnd ? _assembly.adQPointsFace() : _assembly.adQPoints()),
60 4731 : _rz_radial_coord(_mesh.getAxisymmetricRadialCoord()),
61 9858 : _rz_axial_coord(_rz_radial_coord == 0 ? 1 : 0)
62 : {
63 9462 : if (!_fe_problem.hasUserObject("ins_ad_object_tracker"))
64 : {
65 1281 : InputParameters tracker_params = INSADObjectTracker::validParams();
66 1281 : tracker_params.addPrivateParam("_moose_app", &_app);
67 :
68 2562 : _fe_problem.addUserObject("INSADObjectTracker", "ins_ad_object_tracker", tracker_params);
69 1281 : }
70 :
71 : // Bypass the UserObjectInterface method because it requires a UserObjectName param which we
72 : // don't need
73 4731 : _object_tracker = &_fe_problem.getUserObject<INSADObjectTracker>("ins_ad_object_tracker");
74 4731 : const_cast<INSADObjectTracker *>(_object_tracker)->addBlockIDs(this->blockIDs());
75 4731 : }
76 :
77 : void
78 4731 : INSADMaterial::resolveOptionalProperties()
79 : {
80 4731 : Material::resolveOptionalProperties();
81 :
82 9528 : for (const auto sub_id : this->blockIDs())
83 4797 : if ((_object_tracker->get<bool>("has_boussinesq", sub_id)))
84 : {
85 924 : _boussinesq_alphas[sub_id] = &getPossiblyConstantGenericMaterialPropertyByName<Real, true>(
86 462 : _object_tracker->get<MaterialPropertyName>("alpha", sub_id), _material_data, 0);
87 924 : _ref_temps[sub_id] = &getPossiblyConstantGenericMaterialPropertyByName<Real, false>(
88 462 : _object_tracker->get<MaterialPropertyName>("ref_temp", sub_id), _material_data, 0);
89 : }
90 4731 : }
91 :
92 : void
93 662442 : INSADMaterial::subdomainSetup()
94 : {
95 662442 : if ((_has_transient = _object_tracker->get<bool>("has_transient", _current_subdomain_id)))
96 169187 : _velocity_dot = &adCoupledVectorDot("velocity");
97 : else
98 493255 : _velocity_dot = nullptr;
99 :
100 662442 : if (auto alpha_it = _boussinesq_alphas.find(_current_subdomain_id);
101 : alpha_it != _boussinesq_alphas.end())
102 : {
103 350113 : _has_boussinesq = true;
104 350113 : _boussinesq_alpha = alpha_it->second;
105 350113 : auto & temp_var = _subproblem.getStandardVariable(
106 350113 : _tid, _object_tracker->get<std::string>("temperature", _current_subdomain_id));
107 350113 : addMooseVariableDependency(&temp_var);
108 350113 : _temperature = &temp_var.adSln();
109 350113 : _ref_temp = libmesh_map_find(_ref_temps, _current_subdomain_id);
110 : }
111 : else
112 : {
113 312329 : _has_boussinesq = false;
114 312329 : _boussinesq_alpha = nullptr;
115 312329 : _temperature = nullptr;
116 312329 : _ref_temp = nullptr;
117 : }
118 :
119 662442 : _has_gravity = _object_tracker->get<bool>("has_gravity", _current_subdomain_id);
120 662442 : if (_has_gravity || _has_boussinesq)
121 350558 : _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 662442 : if ((_has_advected_mesh = _object_tracker->get<bool>("has_advected_mesh", _current_subdomain_id)))
129 : {
130 5364 : auto & disp_x = _subproblem.getStandardVariable(
131 5364 : _tid, _object_tracker->get<VariableName>("disp_x", _current_subdomain_id));
132 5364 : addMooseVariableDependency(&disp_x);
133 5364 : _disp_x_dot = &disp_x.adUDot();
134 5364 : _disp_x_sys_num = disp_x.sys().number();
135 5364 : _disp_x_num = (disp_x.kind() == Moose::VarKindType::VAR_SOLVER) &&
136 5364 : (_disp_x_sys_num == _fe_problem.currentNonlinearSystem().number())
137 5364 : ? disp_x.number()
138 : : libMesh::invalid_uint;
139 10728 : if (_object_tracker->isTrackerParamValid("disp_y", _current_subdomain_id))
140 : {
141 5364 : auto & disp_y = _subproblem.getStandardVariable(
142 5364 : _tid, _object_tracker->get<VariableName>("disp_y", _current_subdomain_id));
143 5364 : addMooseVariableDependency(&disp_y);
144 5364 : _disp_y_dot = &disp_y.adUDot();
145 5364 : _disp_y_sys_num = disp_y.sys().number();
146 5364 : _disp_y_num =
147 : disp_y.kind() == (Moose::VarKindType::VAR_SOLVER &&
148 : (_disp_y_sys_num == _fe_problem.currentNonlinearSystem().number()))
149 5364 : ? 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 10728 : if (_object_tracker->isTrackerParamValid("disp_z", _current_subdomain_id))
159 : {
160 2628 : auto & disp_z = _subproblem.getStandardVariable(
161 2628 : _tid, _object_tracker->get<VariableName>("disp_z", _current_subdomain_id));
162 2628 : addMooseVariableDependency(&disp_z);
163 2628 : _disp_z_dot = &disp_z.adUDot();
164 2628 : _disp_z_sys_num = disp_z.sys().number();
165 2628 : _disp_z_num =
166 : disp_z.kind() == (Moose::VarKindType::VAR_SOLVER &&
167 : (_disp_z_sys_num == _fe_problem.currentNonlinearSystem().number()))
168 2628 : ? disp_z.number()
169 : : libMesh::invalid_uint;
170 : }
171 : else
172 : {
173 2736 : _disp_z_dot = nullptr;
174 2736 : _disp_z_sys_num = libMesh::invalid_uint;
175 2736 : _disp_z_num = libMesh::invalid_uint;
176 : }
177 : }
178 : else
179 657078 : _disp_x_dot = _disp_y_dot = _disp_z_dot = nullptr;
180 :
181 662442 : _viscous_form = static_cast<NS::ViscousForm>(
182 662442 : int(_object_tracker->get<MooseEnum>("viscous_form", _current_subdomain_id)));
183 :
184 662442 : if ((_has_coupled_force = _object_tracker->get<bool>("has_coupled_force", _current_subdomain_id)))
185 : {
186 71443 : _coupled_force_var.clear();
187 71443 : _coupled_force_vector_function.clear();
188 142886 : if (_object_tracker->isTrackerParamValid("coupled_force_var", _current_subdomain_id))
189 : {
190 72778 : const auto & var_names = _object_tracker->get<std::vector<VariableName>>(
191 36389 : "coupled_force_var", _current_subdomain_id);
192 73668 : for (const auto & var_name : var_names)
193 37279 : _coupled_force_var.push_back(&_subproblem.getVectorVariable(_tid, var_name).adSln());
194 : }
195 :
196 142886 : if (_object_tracker->isTrackerParamValid("coupled_force_vector_function",
197 71443 : _current_subdomain_id))
198 : {
199 71888 : const auto & func_names = _object_tracker->get<std::vector<FunctionName>>(
200 35944 : "coupled_force_vector_function", _current_subdomain_id);
201 72778 : for (const auto & func_name : func_names)
202 36834 : _coupled_force_vector_function.push_back(&_fe_problem.getFunction(func_name, _tid));
203 : }
204 : }
205 662442 : }
206 :
207 : void
208 32242651 : INSADMaterial::computeQpProperties()
209 : {
210 32242651 : _mass_strong_residual[_qp] = -NS::divergence(
211 64485302 : _grad_velocity[_qp], _velocity[_qp], _ad_q_point[_qp], _coord_sys, _rz_radial_coord);
212 :
213 64485302 : _advective_strong_residual[_qp] = _rho[_qp] * _grad_velocity[_qp] * _velocity[_qp];
214 32242651 : if (_has_transient)
215 13909142 : _td_strong_residual[_qp] = _rho[_qp] * (*_velocity_dot)[_qp];
216 32242651 : if (_has_gravity)
217 31054272 : _gravity_strong_residual[_qp] = -_rho[_qp] * _gravity_vector;
218 32242651 : if (_has_boussinesq)
219 30925248 : _boussinesq_strong_residual[_qp] = (*_boussinesq_alpha)[_qp] * _gravity_vector * _rho[_qp] *
220 30925248 : ((*_temperature)[_qp] - (*_ref_temp)[_qp]);
221 32242651 : _relative_velocity[_qp] = _velocity[_qp];
222 :
223 32242651 : if (_has_advected_mesh)
224 : {
225 690184 : _mesh_velocity[_qp](0) = (*_disp_x_dot)[_qp];
226 690184 : if (_disp_y_dot)
227 690184 : _mesh_velocity[_qp](1) = (*_disp_y_dot)[_qp];
228 690184 : if (_disp_z_dot)
229 282128 : _mesh_velocity[_qp](2) = (*_disp_z_dot)[_qp];
230 690184 : _relative_velocity[_qp] -= _mesh_velocity[_qp];
231 2070552 : _advected_mesh_strong_residual[_qp] = -_rho[_qp] * _grad_velocity[_qp] * _mesh_velocity[_qp];
232 : }
233 :
234 32242651 : if (_has_coupled_force)
235 : {
236 1668432 : _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 2728440 : for (const auto * var : _coupled_force_var)
241 : {
242 : mooseAssert(var, "null coupled variable in INSADMaterial");
243 1060008 : _coupled_force_strong_residual[_qp] -= (*var)[_qp];
244 : }
245 2663928 : for (const auto * fn : _coupled_force_vector_function)
246 : {
247 : mooseAssert(fn, "null coupled function in INSADMaterial");
248 1990992 : _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 32242651 : }
|
