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 425644 : INSADMaterial::subdomainSetup()
94 : {
95 425644 : if ((_has_transient = _object_tracker->get<bool>("has_transient", _current_subdomain_id)))
96 113147 : _velocity_dot = &adCoupledVectorDot("velocity");
97 : else
98 312497 : _velocity_dot = nullptr;
99 :
100 425644 : if (auto alpha_it = _boussinesq_alphas.find(_current_subdomain_id);
101 : alpha_it != _boussinesq_alphas.end())
102 : {
103 225781 : _has_boussinesq = true;
104 225781 : _boussinesq_alpha = alpha_it->second;
105 225781 : auto & temp_var = _subproblem.getStandardVariable(
106 225781 : _tid, _object_tracker->get<std::string>("temperature", _current_subdomain_id));
107 225781 : addMooseVariableDependency(&temp_var);
108 225781 : _temperature = &temp_var.adSln();
109 225781 : _ref_temp = libmesh_map_find(_ref_temps, _current_subdomain_id);
110 : }
111 : else
112 : {
113 199863 : _has_boussinesq = false;
114 199863 : _boussinesq_alpha = nullptr;
115 199863 : _temperature = nullptr;
116 199863 : _ref_temp = nullptr;
117 : }
118 :
119 425644 : _has_gravity = _object_tracker->get<bool>("has_gravity", _current_subdomain_id);
120 425644 : if (_has_gravity || _has_boussinesq)
121 226010 : _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 425644 : if ((_has_advected_mesh = _object_tracker->get<bool>("has_advected_mesh", _current_subdomain_id)))
129 : {
130 2852 : auto & disp_x = _subproblem.getStandardVariable(
131 2852 : _tid, _object_tracker->get<VariableName>("disp_x", _current_subdomain_id));
132 2852 : addMooseVariableDependency(&disp_x);
133 2852 : _disp_x_dot = &disp_x.adUDot();
134 2852 : _disp_x_sys_num = disp_x.sys().number();
135 2852 : _disp_x_num = (disp_x.kind() == Moose::VarKindType::VAR_SOLVER) &&
136 2852 : (_disp_x_sys_num == _fe_problem.currentNonlinearSystem().number())
137 2852 : ? disp_x.number()
138 : : libMesh::invalid_uint;
139 5704 : if (_object_tracker->isTrackerParamValid("disp_y", _current_subdomain_id))
140 : {
141 2852 : auto & disp_y = _subproblem.getStandardVariable(
142 2852 : _tid, _object_tracker->get<VariableName>("disp_y", _current_subdomain_id));
143 2852 : addMooseVariableDependency(&disp_y);
144 2852 : _disp_y_dot = &disp_y.adUDot();
145 2852 : _disp_y_sys_num = disp_y.sys().number();
146 2852 : _disp_y_num =
147 : disp_y.kind() == (Moose::VarKindType::VAR_SOLVER &&
148 : (_disp_y_sys_num == _fe_problem.currentNonlinearSystem().number()))
149 2852 : ? 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 5704 : if (_object_tracker->isTrackerParamValid("disp_z", _current_subdomain_id))
159 : {
160 1396 : auto & disp_z = _subproblem.getStandardVariable(
161 1396 : _tid, _object_tracker->get<VariableName>("disp_z", _current_subdomain_id));
162 1396 : addMooseVariableDependency(&disp_z);
163 1396 : _disp_z_dot = &disp_z.adUDot();
164 1396 : _disp_z_sys_num = disp_z.sys().number();
165 1396 : _disp_z_num =
166 : disp_z.kind() == (Moose::VarKindType::VAR_SOLVER &&
167 : (_disp_z_sys_num == _fe_problem.currentNonlinearSystem().number()))
168 1396 : ? disp_z.number()
169 : : libMesh::invalid_uint;
170 : }
171 : else
172 : {
173 1456 : _disp_z_dot = nullptr;
174 1456 : _disp_z_sys_num = libMesh::invalid_uint;
175 1456 : _disp_z_num = libMesh::invalid_uint;
176 : }
177 : }
178 : else
179 422792 : _disp_x_dot = _disp_y_dot = _disp_z_dot = nullptr;
180 :
181 425644 : _viscous_form = static_cast<NS::ViscousForm>(
182 425644 : int(_object_tracker->get<MooseEnum>("viscous_form", _current_subdomain_id)));
183 :
184 425644 : if ((_has_coupled_force = _object_tracker->get<bool>("has_coupled_force", _current_subdomain_id)))
185 : {
186 42459 : _coupled_force_var.clear();
187 42459 : _coupled_force_vector_function.clear();
188 84918 : if (_object_tracker->isTrackerParamValid("coupled_force_var", _current_subdomain_id))
189 : {
190 43146 : const auto & var_names = _object_tracker->get<std::vector<VariableName>>(
191 21573 : "coupled_force_var", _current_subdomain_id);
192 43604 : for (const auto & var_name : var_names)
193 22031 : _coupled_force_var.push_back(&_subproblem.getVectorVariable(_tid, var_name).adSln());
194 : }
195 :
196 84918 : if (_object_tracker->isTrackerParamValid("coupled_force_vector_function",
197 42459 : _current_subdomain_id))
198 : {
199 42688 : const auto & func_names = _object_tracker->get<std::vector<FunctionName>>(
200 21344 : "coupled_force_vector_function", _current_subdomain_id);
201 43146 : for (const auto & func_name : func_names)
202 21802 : _coupled_force_vector_function.push_back(&_fe_problem.getFunction(func_name, _tid));
203 : }
204 : }
205 425644 : }
206 :
207 : void
208 24590042 : INSADMaterial::computeQpProperties()
209 : {
210 24590042 : _mass_strong_residual[_qp] = -NS::divergence(
211 49180084 : _grad_velocity[_qp], _velocity[_qp], _ad_q_point[_qp], _coord_sys, _rz_radial_coord);
212 :
213 49180084 : _advective_strong_residual[_qp] = _rho[_qp] * _grad_velocity[_qp] * _velocity[_qp];
214 24590042 : if (_has_transient)
215 9857884 : _td_strong_residual[_qp] = _rho[_qp] * (*_velocity_dot)[_qp];
216 24590042 : if (_has_gravity)
217 24525984 : _gravity_strong_residual[_qp] = -_rho[_qp] * _gravity_vector;
218 24590042 : if (_has_boussinesq)
219 24439968 : _boussinesq_strong_residual[_qp] = (*_boussinesq_alpha)[_qp] * _gravity_vector * _rho[_qp] *
220 24439968 : ((*_temperature)[_qp] - (*_ref_temp)[_qp]);
221 24590042 : _relative_velocity[_qp] = _velocity[_qp];
222 :
223 24590042 : 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 24590042 : 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 24590042 : }
|
