From 36d097cc4d9e84143412e748c8f631092cbd73d6 Mon Sep 17 00:00:00 2001 From: Guillaume Giudicelli Date: Fri, 10 Jan 2025 14:23:13 -0700 Subject: [PATCH 01/18] Simplify a few loops in MooseVariableData refs #29696 --- framework/src/variables/MooseVariableData.C | 26 +++++++++++---------- 1 file changed, 14 insertions(+), 12 deletions(-) diff --git a/framework/src/variables/MooseVariableData.C b/framework/src/variables/MooseVariableData.C index 4a83bb5a90a4..34a7ae2cdb73 100644 --- a/framework/src/variables/MooseVariableData.C +++ b/framework/src/variables/MooseVariableData.C @@ -1306,13 +1306,13 @@ MooseVariableData::setDofValue(const OutputData & value, unsigned in _has_dof_values = true; auto & u = _vector_tag_u[_solution_tag]; - for (unsigned int qp = 0; qp < u.size(); qp++) - { - u[qp] = (*_phi)[0][qp] * dof_values[0]; - - for (unsigned int i = 1; i < dof_values.size(); i++) + const auto nqps = u.size(); + const auto ndofs = dof_values.size(); + for (unsigned int qp = 0; qp < nqps; qp++) + u[qp] *= 0.; + for (unsigned int qp = 0; qp < nqps; qp++) + for (unsigned int i = 0; i < ndofs; i++) u[qp] += (*_phi)[i][qp] * dof_values[i]; - } } template @@ -1326,12 +1326,13 @@ MooseVariableData::setDofValues(const DenseVector & valu _has_dof_values = true; auto & u = _vector_tag_u[_solution_tag]; - for (unsigned int qp = 0; qp < u.size(); qp++) - { - u[qp] = (*_phi)[0][qp] * dof_values[0]; - for (unsigned int i = 1; i < dof_values.size(); i++) + const auto nqps = u.size(); + const auto ndofs = dof_values.size(); + for (unsigned int qp = 0; qp < nqps; qp++) + u[qp] *= 0.; + for (unsigned int qp = 0; qp < nqps; qp++) + for (unsigned int i = 0; i < ndofs; i++) u[qp] += (*_phi)[i][qp] * dof_values[i]; - } } template @@ -1669,10 +1670,11 @@ MooseVariableData::computeIncrementAtNode( else { unsigned int n = 0; + const auto n_dof_indices = _dof_indices.size(); for (unsigned int j = 0; j < _count; j++) { _increment[0](j) = increment_vec(_dof_indices[0] + n); - n += _dof_indices.size(); + n += n_dof_indices; } } } From ef99ff152fab8b82c51884a3fa00bf4eb9d25fff Mon Sep 17 00:00:00 2001 From: Guillaume Giudicelli Date: Fri, 10 Jan 2025 15:12:30 -0700 Subject: [PATCH 02/18] more memory locality and less if statements refs #29696 --- framework/src/variables/MooseVariableData.C | 363 ++++++++++---------- 1 file changed, 172 insertions(+), 191 deletions(-) diff --git a/framework/src/variables/MooseVariableData.C b/framework/src/variables/MooseVariableData.C index 34a7ae2cdb73..e6234789b582 100644 --- a/framework/src/variables/MooseVariableData.C +++ b/framework/src/variables/MooseVariableData.C @@ -428,240 +428,221 @@ MooseVariableData::computeValues() unsigned int nqp = _current_qrule->n_points(); auto && active_coupleable_matrix_tags = _subproblem.getActiveFEVariableCoupleableMatrixTags(_tid); - for (auto tag : _required_vector_tags) - { - if (_need_vector_tag_u[tag]) - _vector_tag_u[tag].resize(nqp); - if (_need_vector_tag_grad[tag]) - _vector_tag_grad[tag].resize(nqp); - } - - for (auto tag : active_coupleable_matrix_tags) - if (_need_matrix_tag_u[tag]) - _matrix_tag_u[tag].resize(nqp); - - if (_need_second) - _second_u.resize(nqp); + bool second_required = + _need_second || _need_second_old || _need_second_older || _need_second_previous_nl; + bool curl_required = _need_curl || _need_curl_old; + bool div_required = _need_div || _need_div_old; + if (second_required) + mooseAssert(_current_second_phi, + "We're requiring a second calculation but have not set a second shape function!"); + if (curl_required) + mooseAssert(_current_curl_phi, + "We're requiring a curl calculation but have not set a curl shape function!"); + if (div_required) + mooseAssert(_current_div_phi, + "We're requiring a divergence calculation but have not set a div shape function!"); + + // Curl if (_need_curl) + { _curl_u.resize(nqp); + for (unsigned int i = 0; i < nqp; ++i) + _curl_u[i] = 0; + for (unsigned int i = 0; i < num_dofs; i++) + for (unsigned int qp = 0; qp < nqp; qp++) + _curl_u[qp] += (*_current_curl_phi)[i][qp] * _vector_tags_dof_u[_solution_tag][i]; + } + if (_need_curl_old) + { + _curl_u_old.resize(nqp); + for (unsigned int i = 0; i < nqp; ++i) + _curl_u_old[i] = 0; + for (unsigned int i = 0; i < num_dofs; i++) + for (unsigned int qp = 0; qp < nqp; qp++) + _curl_u_old[qp] += (*_current_curl_phi)[i][qp] * _vector_tags_dof_u[_old_solution_tag][i]; + } + // Div if (_need_div) + { _div_u.resize(nqp); + for (unsigned int i = 0; i < nqp; ++i) + _div_u[i] = 0; + for (unsigned int i = 0; i < num_dofs; i++) + for (unsigned int qp = 0; qp < nqp; qp++) + _div_u[qp] += (*_current_div_phi)[i][qp] * _vector_tags_dof_u[_solution_tag][i]; + } + if (is_transient && _need_div_old) + { + _div_u_old.resize(nqp); + for (unsigned int i = 0; i < nqp; ++i) + _div_u_old[i] = 0; + for (unsigned int i = 0; i < num_dofs; i++) + for (unsigned int qp = 0; qp < nqp; qp++) + _div_u_old[qp] += (*_current_div_phi)[i][qp] * _vector_tags_dof_u[_old_solution_tag][i]; + } + // Second + if (_need_second) + { + _second_u.resize(nqp); + for (unsigned int i = 0; i < nqp; ++i) + _second_u[i] = 0; + for (unsigned int i = 0; i < num_dofs; i++) + for (unsigned int qp = 0; qp < nqp; qp++) + _second_u[qp].add_scaled((*_current_second_phi)[i][qp], + _vector_tags_dof_u[_solution_tag][i]); + } if (_need_second_previous_nl) + { _second_u_previous_nl.resize(nqp); - + for (unsigned int i = 0; i < nqp; ++i) + _second_u_previous_nl[i] = 0; + for (unsigned int i = 0; i < num_dofs; i++) + for (unsigned int qp = 0; qp < nqp; qp++) + _second_u_previous_nl[qp].add_scaled((*_current_second_phi)[i][qp], + _vector_tags_dof_u[_previous_nl_solution_tag][i]); + } if (is_transient) { - if (_need_u_dot) - _u_dot.resize(nqp); - - if (_need_u_dotdot) - _u_dotdot.resize(nqp); - - if (_need_u_dot_old) - _u_dot_old.resize(nqp); - - if (_need_u_dotdot_old) - _u_dotdot_old.resize(nqp); - - if (_need_du_dot_du) - _du_dot_du.resize(nqp); - - if (_need_du_dotdot_du) - _du_dotdot_du.resize(nqp); - - if (_need_grad_dot) - _grad_u_dot.resize(nqp); - - if (_need_grad_dotdot) - _grad_u_dotdot.resize(nqp); - if (_need_second_old) + { _second_u_old.resize(nqp); - - if (_need_curl_old) - _curl_u_old.resize(nqp); - - if (_need_div_old) - _div_u_old.resize(nqp); + for (unsigned int i = 0; i < nqp; ++i) + _second_u_old[i] = 0; + for (unsigned int i = 0; i < num_dofs; i++) + for (unsigned int qp = 0; qp < nqp; qp++) + _second_u_old[qp].add_scaled((*_current_second_phi)[i][qp], + _vector_tags_dof_u[_old_solution_tag][i]); + } if (_need_second_older) + { _second_u_older.resize(nqp); + for (unsigned int i = 0; i < nqp; ++i) + _second_u_older[i] = 0; + for (unsigned int i = 0; i < num_dofs; i++) + for (unsigned int qp = 0; qp < nqp; qp++) + _second_u_older[qp].add_scaled((*_current_second_phi)[i][qp], + _vector_tags_dof_u[_older_solution_tag][i]); + } } - for (unsigned int i = 0; i < nqp; ++i) + for (auto tag : _required_vector_tags) { - for (auto tag : _required_vector_tags) + if (_need_vector_tag_u[tag] && _sys.hasVector(tag) && _sys.getVector(tag).closed()) { - if (_need_vector_tag_u[tag]) + _vector_tag_u[tag].resize(nqp); + for (unsigned int i = 0; i < nqp; ++i) _vector_tag_u[tag][i] = 0; - if (_need_vector_tag_grad[tag]) + for (unsigned int i = 0; i < num_dofs; i++) + for (unsigned int qp = 0; qp < nqp; qp++) + _vector_tag_u[tag][qp] += (*_current_phi)[i][qp] * _vector_tags_dof_u[tag][i]; + } + if (_need_vector_tag_grad[tag] && _sys.hasVector(tag) && _sys.getVector(tag).closed()) + { + _vector_tag_grad[tag].resize(nqp); + for (unsigned int i = 0; i < nqp; ++i) _vector_tag_grad[tag][i] = 0; + for (unsigned int i = 0; i < num_dofs; i++) + for (unsigned int qp = 0; qp < nqp; qp++) + _vector_tag_grad[tag][qp].add_scaled((*_current_grad_phi)[i][qp], + _vector_tags_dof_u[tag][i]); } + } - for (auto tag : active_coupleable_matrix_tags) - if (_need_matrix_tag_u[tag]) + for (auto tag : active_coupleable_matrix_tags) + if (_need_matrix_tag_u[tag]) + { + _matrix_tag_u[tag].resize(nqp); + for (unsigned int i = 0; i < nqp; ++i) _matrix_tag_u[tag][i] = 0; + for (unsigned int i = 0; i < num_dofs; i++) + for (unsigned int qp = 0; qp < nqp; qp++) + _matrix_tag_u[tag][qp] += (*_current_phi)[i][qp] * _matrix_tags_dof_u[tag][i]; + } - if (_need_second) - _second_u[i] = 0; - - if (_need_curl) - _curl_u[i] = 0; - - if (_need_div) - _div_u[i] = 0; - - if (_need_second_previous_nl) - _second_u_previous_nl[i] = 0; - - if (is_transient) + if (is_transient) + { + if (_need_u_dot) { - if (_need_u_dot) + _u_dot.resize(nqp); + for (unsigned int i = 0; i < nqp; ++i) _u_dot[i] = 0; + for (unsigned int i = 0; i < num_dofs; i++) + for (unsigned int qp = 0; qp < nqp; qp++) + _u_dot[qp] += (*_current_phi)[i][qp] * _dof_values_dot[i]; + } - if (_need_u_dotdot) + if (_need_u_dotdot) + { + _u_dotdot.resize(nqp); + for (unsigned int i = 0; i < nqp; ++i) _u_dotdot[i] = 0; + for (unsigned int i = 0; i < num_dofs; i++) + for (unsigned int qp = 0; qp < nqp; qp++) + _u_dotdot[qp] += (*_current_phi)[i][qp] * _dof_values_dotdot[i]; + } - if (_need_u_dot_old) + if (_need_u_dot_old) + { + _u_dot_old.resize(nqp); + for (unsigned int i = 0; i < nqp; ++i) _u_dot_old[i] = 0; + for (unsigned int i = 0; i < num_dofs; i++) + for (unsigned int qp = 0; qp < nqp; qp++) + _u_dot_old[qp] += (*_current_phi)[i][qp] * _dof_values_dot_old[i]; + } - if (_need_u_dotdot_old) + if (_need_u_dotdot_old) + { + _u_dotdot_old.resize(nqp); + for (unsigned int i = 0; i < nqp; ++i) _u_dotdot_old[i] = 0; + for (unsigned int i = 0; i < num_dofs; i++) + for (unsigned int qp = 0; qp < nqp; qp++) + _u_dotdot_old[qp] += (*_current_phi)[i][qp] * _dof_values_dotdot_old[i]; + } - if (_need_du_dot_du) + if (_need_du_dot_du) + { + _du_dot_du.resize(nqp); + for (unsigned int i = 0; i < nqp; ++i) _du_dot_du[i] = 0; + for (unsigned int i = 0; i < num_dofs; i++) + for (unsigned int qp = 0; qp < nqp; qp++) + _du_dot_du[qp] = _dof_du_dot_du[i]; + } - if (_need_du_dotdot_du) + if (_need_du_dotdot_du) + { + _du_dotdot_du.resize(nqp); + for (unsigned int i = 0; i < nqp; ++i) _du_dotdot_du[i] = 0; + for (unsigned int i = 0; i < num_dofs; i++) + for (unsigned int qp = 0; qp < nqp; qp++) + _du_dotdot_du[qp] = _dof_du_dotdot_du[i]; + } - if (_need_grad_dot) + if (_need_grad_dot) + { + _grad_u_dot.resize(nqp); + for (unsigned int i = 0; i < nqp; ++i) _grad_u_dot[i] = 0; - - if (_need_grad_dotdot) - _grad_u_dotdot[i] = 0; - - if (_need_second_old) - _second_u_old[i] = 0; - - if (_need_second_older) - _second_u_older[i] = 0; - - if (_need_curl_old) - _curl_u_old[i] = 0; - - if (_need_div_old) - _div_u_old[i] = 0; + for (unsigned int i = 0; i < num_dofs; i++) + for (unsigned int qp = 0; qp < nqp; qp++) + _grad_u_dot[qp].add_scaled((*_current_grad_phi)[i][qp], _dof_values_dot[i]); } - } - - bool second_required = - _need_second || _need_second_old || _need_second_older || _need_second_previous_nl; - bool curl_required = _need_curl || _need_curl_old; - bool div_required = _need_div || _need_div_old; - for (unsigned int i = 0; i < num_dofs; i++) - { - for (unsigned int qp = 0; qp < nqp; qp++) + if (_need_grad_dotdot) { - const OutputType phi_local = (*_current_phi)[i][qp]; - const typename OutputTools::OutputGradient dphi_qp = (*_current_grad_phi)[i][qp]; - - if (is_transient) - { - if (_need_u_dot) - _u_dot[qp] += phi_local * _dof_values_dot[i]; - - if (_need_u_dotdot) - _u_dotdot[qp] += phi_local * _dof_values_dotdot[i]; - - if (_need_u_dot_old) - _u_dot_old[qp] += phi_local * _dof_values_dot_old[i]; - - if (_need_u_dotdot_old) - _u_dotdot_old[qp] += phi_local * _dof_values_dotdot_old[i]; - - if (_need_grad_dot) - _grad_u_dot[qp].add_scaled(dphi_qp, _dof_values_dot[i]); - - if (_need_grad_dotdot) - _grad_u_dotdot[qp].add_scaled(dphi_qp, _dof_values_dotdot[i]); - - if (_need_du_dot_du) - _du_dot_du[qp] = _dof_du_dot_du[i]; - - if (_need_du_dotdot_du) - _du_dotdot_du[qp] = _dof_du_dotdot_du[i]; - } - - if (second_required) - { - mooseAssert( - _current_second_phi, - "We're requiring a second calculation but have not set a second shape function!"); - const typename OutputTools::OutputSecond d2phi_local = - (*_current_second_phi)[i][qp]; - - if (_need_second) - _second_u[qp].add_scaled(d2phi_local, _vector_tags_dof_u[_solution_tag][i]); - - if (_need_second_previous_nl) - _second_u_previous_nl[qp].add_scaled(d2phi_local, - _vector_tags_dof_u[_previous_nl_solution_tag][i]); - - if (is_transient) - { - if (_need_second_old) - _second_u_old[qp].add_scaled(d2phi_local, _vector_tags_dof_u[_old_solution_tag][i]); - - if (_need_second_older) - _second_u_older[qp].add_scaled(d2phi_local, _vector_tags_dof_u[_older_solution_tag][i]); - } - } - - if (curl_required) - { - mooseAssert(_current_curl_phi, - "We're requiring a curl calculation but have not set a curl shape function!"); - const OutputType curl_phi_local = (*_current_curl_phi)[i][qp]; - - if (_need_curl) - _curl_u[qp] += curl_phi_local * _vector_tags_dof_u[_solution_tag][i]; - - if (is_transient && _need_curl_old) - _curl_u_old[qp] += curl_phi_local * _vector_tags_dof_u[_old_solution_tag][i]; - } - - if (div_required) - { - mooseAssert( - _current_div_phi, - "We're requiring a divergence calculation but have not set a div shape function!"); - const OutputShapeDivergence div_phi_local = (*_current_div_phi)[i][qp]; - - if (_need_div) - _div_u[qp] += div_phi_local * _vector_tags_dof_u[_solution_tag][i]; - - if (is_transient && _need_div_old) - _div_u_old[qp] += div_phi_local * _vector_tags_dof_u[_old_solution_tag][i]; - } - - for (auto tag : _required_vector_tags) - { - if (_sys.hasVector(tag) && _sys.getVector(tag).closed()) - { - if (_need_vector_tag_u[tag]) - _vector_tag_u[tag][qp] += phi_local * _vector_tags_dof_u[tag][i]; - if (_need_vector_tag_grad[tag]) - _vector_tag_grad[tag][qp].add_scaled(dphi_qp, _vector_tags_dof_u[tag][i]); - } - } - - for (auto tag : active_coupleable_matrix_tags) - if (_need_matrix_tag_u[tag]) - _matrix_tag_u[tag][qp] += phi_local * _matrix_tags_dof_u[tag][i]; + _grad_u_dotdot.resize(nqp); + for (unsigned int i = 0; i < nqp; ++i) + _grad_u_dotdot[i] = 0; + for (unsigned int i = 0; i < num_dofs; i++) + for (unsigned int qp = 0; qp < nqp; qp++) + _grad_u_dotdot[qp].add_scaled((*_current_grad_phi)[i][qp], _dof_values_dotdot[i]); } } From 824d0e8f490253836318befff87fedff4b1ddeb3 Mon Sep 17 00:00:00 2001 From: Guillaume Giudicelli Date: Sun, 26 Jan 2025 22:29:54 -0700 Subject: [PATCH 03/18] Avoid implicit conversions which prevent vectorization And remove a few harmless others Avoid a few size computations --- framework/src/variables/MooseVariableData.C | 88 +++++++++---------- framework/src/variables/MooseVariableDataFV.C | 12 +-- 2 files changed, 48 insertions(+), 52 deletions(-) diff --git a/framework/src/variables/MooseVariableData.C b/framework/src/variables/MooseVariableData.C index e6234789b582..dfac73e4cb53 100644 --- a/framework/src/variables/MooseVariableData.C +++ b/framework/src/variables/MooseVariableData.C @@ -428,27 +428,21 @@ MooseVariableData::computeValues() unsigned int nqp = _current_qrule->n_points(); auto && active_coupleable_matrix_tags = _subproblem.getActiveFEVariableCoupleableMatrixTags(_tid); - bool second_required = - _need_second || _need_second_old || _need_second_older || _need_second_previous_nl; - bool curl_required = _need_curl || _need_curl_old; - bool div_required = _need_div || _need_div_old; - - if (second_required) - mooseAssert(_current_second_phi, - "We're requiring a second calculation but have not set a second shape function!"); - if (curl_required) - mooseAssert(_current_curl_phi, - "We're requiring a curl calculation but have not set a curl shape function!"); - if (div_required) - mooseAssert(_current_div_phi, - "We're requiring a divergence calculation but have not set a div shape function!"); + mooseAssert( + !(_need_second || _need_second_old || _need_second_older || _need_second_previous_nl) || + _current_second_phi, + "We're requiring a second calculation but have not set a second shape function!"); + mooseAssert(!(_need_curl || _need_curl_old) || _current_curl_phi, + "We're requiring a curl calculation but have not set a curl shape function!"); + mooseAssert(!(_need_div || _need_div_old) || _current_div_phi, + "We're requiring a divergence calculation but have not set a div shape function!"); // Curl if (_need_curl) { _curl_u.resize(nqp); for (unsigned int i = 0; i < nqp; ++i) - _curl_u[i] = 0; + _curl_u[i] = 0.; for (unsigned int i = 0; i < num_dofs; i++) for (unsigned int qp = 0; qp < nqp; qp++) _curl_u[qp] += (*_current_curl_phi)[i][qp] * _vector_tags_dof_u[_solution_tag][i]; @@ -457,7 +451,7 @@ MooseVariableData::computeValues() { _curl_u_old.resize(nqp); for (unsigned int i = 0; i < nqp; ++i) - _curl_u_old[i] = 0; + _curl_u_old[i] = 0.; for (unsigned int i = 0; i < num_dofs; i++) for (unsigned int qp = 0; qp < nqp; qp++) _curl_u_old[qp] += (*_current_curl_phi)[i][qp] * _vector_tags_dof_u[_old_solution_tag][i]; @@ -468,7 +462,7 @@ MooseVariableData::computeValues() { _div_u.resize(nqp); for (unsigned int i = 0; i < nqp; ++i) - _div_u[i] = 0; + _div_u[i] = 0.; for (unsigned int i = 0; i < num_dofs; i++) for (unsigned int qp = 0; qp < nqp; qp++) _div_u[qp] += (*_current_div_phi)[i][qp] * _vector_tags_dof_u[_solution_tag][i]; @@ -477,7 +471,7 @@ MooseVariableData::computeValues() { _div_u_old.resize(nqp); for (unsigned int i = 0; i < nqp; ++i) - _div_u_old[i] = 0; + _div_u_old[i] = 0.; for (unsigned int i = 0; i < num_dofs; i++) for (unsigned int qp = 0; qp < nqp; qp++) _div_u_old[qp] += (*_current_div_phi)[i][qp] * _vector_tags_dof_u[_old_solution_tag][i]; @@ -488,7 +482,7 @@ MooseVariableData::computeValues() { _second_u.resize(nqp); for (unsigned int i = 0; i < nqp; ++i) - _second_u[i] = 0; + _second_u[i] = 0.; for (unsigned int i = 0; i < num_dofs; i++) for (unsigned int qp = 0; qp < nqp; qp++) _second_u[qp].add_scaled((*_current_second_phi)[i][qp], @@ -510,7 +504,7 @@ MooseVariableData::computeValues() { _second_u_old.resize(nqp); for (unsigned int i = 0; i < nqp; ++i) - _second_u_old[i] = 0; + _second_u_old[i] = 0.; for (unsigned int i = 0; i < num_dofs; i++) for (unsigned int qp = 0; qp < nqp; qp++) _second_u_old[qp].add_scaled((*_current_second_phi)[i][qp], @@ -521,7 +515,7 @@ MooseVariableData::computeValues() { _second_u_older.resize(nqp); for (unsigned int i = 0; i < nqp; ++i) - _second_u_older[i] = 0; + _second_u_older[i] = 0.; for (unsigned int i = 0; i < num_dofs; i++) for (unsigned int qp = 0; qp < nqp; qp++) _second_u_older[qp].add_scaled((*_current_second_phi)[i][qp], @@ -535,7 +529,7 @@ MooseVariableData::computeValues() { _vector_tag_u[tag].resize(nqp); for (unsigned int i = 0; i < nqp; ++i) - _vector_tag_u[tag][i] = 0; + _vector_tag_u[tag][i] = 0.; for (unsigned int i = 0; i < num_dofs; i++) for (unsigned int qp = 0; qp < nqp; qp++) _vector_tag_u[tag][qp] += (*_current_phi)[i][qp] * _vector_tags_dof_u[tag][i]; @@ -544,7 +538,7 @@ MooseVariableData::computeValues() { _vector_tag_grad[tag].resize(nqp); for (unsigned int i = 0; i < nqp; ++i) - _vector_tag_grad[tag][i] = 0; + _vector_tag_grad[tag][i] = 0.; for (unsigned int i = 0; i < num_dofs; i++) for (unsigned int qp = 0; qp < nqp; qp++) _vector_tag_grad[tag][qp].add_scaled((*_current_grad_phi)[i][qp], @@ -557,7 +551,7 @@ MooseVariableData::computeValues() { _matrix_tag_u[tag].resize(nqp); for (unsigned int i = 0; i < nqp; ++i) - _matrix_tag_u[tag][i] = 0; + _matrix_tag_u[tag][i] = 0.; for (unsigned int i = 0; i < num_dofs; i++) for (unsigned int qp = 0; qp < nqp; qp++) _matrix_tag_u[tag][qp] += (*_current_phi)[i][qp] * _matrix_tags_dof_u[tag][i]; @@ -569,7 +563,7 @@ MooseVariableData::computeValues() { _u_dot.resize(nqp); for (unsigned int i = 0; i < nqp; ++i) - _u_dot[i] = 0; + _u_dot[i] = 0.; for (unsigned int i = 0; i < num_dofs; i++) for (unsigned int qp = 0; qp < nqp; qp++) _u_dot[qp] += (*_current_phi)[i][qp] * _dof_values_dot[i]; @@ -579,7 +573,7 @@ MooseVariableData::computeValues() { _u_dotdot.resize(nqp); for (unsigned int i = 0; i < nqp; ++i) - _u_dotdot[i] = 0; + _u_dotdot[i] = 0.; for (unsigned int i = 0; i < num_dofs; i++) for (unsigned int qp = 0; qp < nqp; qp++) _u_dotdot[qp] += (*_current_phi)[i][qp] * _dof_values_dotdot[i]; @@ -589,7 +583,7 @@ MooseVariableData::computeValues() { _u_dot_old.resize(nqp); for (unsigned int i = 0; i < nqp; ++i) - _u_dot_old[i] = 0; + _u_dot_old[i] = 0.; for (unsigned int i = 0; i < num_dofs; i++) for (unsigned int qp = 0; qp < nqp; qp++) _u_dot_old[qp] += (*_current_phi)[i][qp] * _dof_values_dot_old[i]; @@ -599,7 +593,7 @@ MooseVariableData::computeValues() { _u_dotdot_old.resize(nqp); for (unsigned int i = 0; i < nqp; ++i) - _u_dotdot_old[i] = 0; + _u_dotdot_old[i] = 0.; for (unsigned int i = 0; i < num_dofs; i++) for (unsigned int qp = 0; qp < nqp; qp++) _u_dotdot_old[qp] += (*_current_phi)[i][qp] * _dof_values_dotdot_old[i]; @@ -609,7 +603,7 @@ MooseVariableData::computeValues() { _du_dot_du.resize(nqp); for (unsigned int i = 0; i < nqp; ++i) - _du_dot_du[i] = 0; + _du_dot_du[i] = 0.; for (unsigned int i = 0; i < num_dofs; i++) for (unsigned int qp = 0; qp < nqp; qp++) _du_dot_du[qp] = _dof_du_dot_du[i]; @@ -619,7 +613,7 @@ MooseVariableData::computeValues() { _du_dotdot_du.resize(nqp); for (unsigned int i = 0; i < nqp; ++i) - _du_dotdot_du[i] = 0; + _du_dotdot_du[i] = 0.; for (unsigned int i = 0; i < num_dofs; i++) for (unsigned int qp = 0; qp < nqp; qp++) _du_dotdot_du[qp] = _dof_du_dotdot_du[i]; @@ -629,7 +623,7 @@ MooseVariableData::computeValues() { _grad_u_dot.resize(nqp); for (unsigned int i = 0; i < nqp; ++i) - _grad_u_dot[i] = 0; + _grad_u_dot[i] = 0.; for (unsigned int i = 0; i < num_dofs; i++) for (unsigned int qp = 0; qp < nqp; qp++) _grad_u_dot[qp].add_scaled((*_current_grad_phi)[i][qp], _dof_values_dot[i]); @@ -639,7 +633,7 @@ MooseVariableData::computeValues() { _grad_u_dotdot.resize(nqp); for (unsigned int i = 0; i < nqp; ++i) - _grad_u_dotdot[i] = 0; + _grad_u_dotdot[i] = 0.; for (unsigned int i = 0; i < num_dofs; i++) for (unsigned int qp = 0; qp < nqp; qp++) _grad_u_dotdot[qp].add_scaled((*_current_grad_phi)[i][qp], _dof_values_dotdot[i]); @@ -794,10 +788,10 @@ MooseVariableData::computeValues() _u_dotdot_old[i].setZero(_count); if (_need_du_dot_du) - _du_dot_du[i] = 0; + _du_dot_du[i] = 0.; if (_need_du_dotdot_du) - _du_dotdot_du[i] = 0; + _du_dotdot_du[i] = 0.; if (_need_grad_dot) _grad_u_dot[i].setZero(_count, LIBMESH_DIM); @@ -1001,10 +995,10 @@ MooseVariableData::computeMonomialValues() } const dof_id_type & idx = _dof_indices[0]; - Real u_dot = 0; - Real u_dotdot = 0; - Real u_dot_old = 0; - Real u_dotdot_old = 0; + Real u_dot = 0.; + Real u_dotdot = 0.; + Real u_dot_old = 0.; + Real u_dotdot_old = 0.; const Real & du_dot_du = _sys.duDotDu(_var_num); const Real & du_dotdot_du = _sys.duDotDotDu(); @@ -1142,7 +1136,7 @@ MooseVariableData::computeAD(const unsigned int num_dofs, const unsi { // Have to do this because upon construction this won't initialize any of the derivatives // (because DualNumber::do_derivatives is false at that time). - _ad_zero = 0; + _ad_zero = 0.; _ad_dof_values.resize(num_dofs); if (_need_ad_u) @@ -1489,8 +1483,9 @@ MooseVariableData::addSolution(NumericVector & sol, unsigned int p = 0; for (unsigned int j = 0; j < _count; ++j) { - unsigned int inc = (isNodal() ? j : j * _dof_indices.size()); - for (unsigned int i = 0; i < _dof_indices.size(); ++i) + const auto num_dofs = _dof_indices.size(); + unsigned int inc = (isNodal() ? j : j * num_dofs); + for (unsigned int i = 0; i < num_dofs; ++i) sol.add(_dof_indices[i] + inc, v(p++)); } } @@ -1581,7 +1576,7 @@ MooseVariableData::computeIncrementAtQps(const NumericVector unsigned int num_dofs = _dof_indices.size(); for (unsigned int qp = 0; qp < nqp; qp++) { - _increment[qp] = 0; + _increment[qp] = 0.; for (unsigned int i = 0; i < num_dofs; i++) _increment[qp] += (*_phi)[i][qp] * increment_vec(_dof_indices[i]); } @@ -1801,7 +1796,7 @@ MooseVariableData::fetchADNodalValues() } // Executing something with a time derivative at initial should not put a NaN else if (_time_integrator && !_time_integrator->dt()) - _ad_dofs_dot[i] = 0; + _ad_dofs_dot[i] = 0.; else mooseError("AD nodal time derivatives not implemented for variables without a time " "integrator (auxiliary variables)"); @@ -1894,11 +1889,12 @@ MooseVariableData::reinitAux() fetchDoFValues(); + const auto num_dofs = _dof_indices.size(); for (auto & dof_u : _vector_tags_dof_u) - dof_u.resize(_dof_indices.size()); + dof_u.resize(num_dofs); for (auto & dof_u : _matrix_tags_dof_u) - dof_u.resize(_dof_indices.size()); + dof_u.resize(num_dofs); _has_dof_indices = true; } @@ -1927,7 +1923,7 @@ MooseVariableData::reinitNodes(const std::vector & node } } - if (_dof_indices.size() > 0) + if (!_dof_indices.empty()) _has_dof_indices = true; else _has_dof_indices = false; diff --git a/framework/src/variables/MooseVariableDataFV.C b/framework/src/variables/MooseVariableDataFV.C index b2bd81a85cd2..2bef995e264c 100644 --- a/framework/src/variables/MooseVariableDataFV.C +++ b/framework/src/variables/MooseVariableDataFV.C @@ -468,28 +468,28 @@ MooseVariableDataFV::computeValues() if (second_required) { if (_need_second) - _second_u[qp] = 0; + _second_u[qp] = 0.; if (_need_second_previous_nl) - _second_u_previous_nl[qp] = 0; + _second_u_previous_nl[qp] = 0.; if (is_transient) { if (_need_second_old) - _second_u_old[qp] = 0; + _second_u_old[qp] = 0.; if (_need_second_older) - _second_u_older[qp] = 0; + _second_u_older[qp] = 0.; } } if (curl_required) { if (_need_curl) - _curl_u[qp] = 0; + _curl_u[qp] = 0.; if (is_transient && _need_curl_old) - _curl_u_old[qp] = 0; + _curl_u_old[qp] = 0.; } for (auto tag : _required_vector_tags) From 2f7663d770e982cdf92fec743923202643a4bf20 Mon Sep 17 00:00:00 2001 From: Guillaume Giudicelli Date: Sun, 26 Jan 2025 22:53:50 -0700 Subject: [PATCH 04/18] Use same technique on computeAD for regular variables --- framework/src/variables/MooseVariableData.C | 196 ++++++++++---------- 1 file changed, 101 insertions(+), 95 deletions(-) diff --git a/framework/src/variables/MooseVariableData.C b/framework/src/variables/MooseVariableData.C index dfac73e4cb53..679a0e8fc562 100644 --- a/framework/src/variables/MooseVariableData.C +++ b/framework/src/variables/MooseVariableData.C @@ -81,6 +81,12 @@ MooseVariableData::MooseVariableData(const MooseVariableField void MooseVariableData::computeAD(const unsigned int num_dofs, const unsigned int nqp) { - // Have to do this because upon construction this won't initialize any of the derivatives - // (because DualNumber::do_derivatives is false at that time). - _ad_zero = 0.; + const bool do_derivatives = Moose::doDerivatives(_subproblem, _sys); _ad_dof_values.resize(num_dofs); - if (_need_ad_u) - _ad_u.resize(nqp); - - if (_need_ad_grad_u) - _ad_grad_u.resize(nqp); - - if (_need_ad_second_u) - _ad_second_u.resize(nqp); + for (unsigned int i = 0; i < num_dofs; i++) + _ad_dof_values[i] = (*_sys.currentSolution())(_dof_indices[i]); + // NOTE! You have to do this AFTER setting the value! + if (do_derivatives) + for (unsigned int i = 0; i < num_dofs; i++) + Moose::derivInsert(_ad_dof_values[i].derivatives(), _dof_indices[i], 1.); - if (_need_ad_u_dot) + if (_need_ad_u) { - _ad_dofs_dot.resize(num_dofs); - _ad_u_dot.resize(nqp); - } - if (_need_ad_grad_u_dot) - _ad_grad_u_dot.resize(nqp); + _ad_u.resize(nqp); + for (unsigned int qp = 0; qp < nqp; qp++) + _ad_u[qp] = _ad_zero; - if (_need_ad_u_dotdot) - { - _ad_dofs_dotdot.resize(num_dofs); - _ad_u_dotdot.resize(nqp); + for (unsigned int i = 0; i < num_dofs; i++) + for (unsigned int qp = 0; qp < nqp; qp++) + _ad_u[qp] += _ad_dof_values[i] * (*_current_phi)[i][qp]; } - const bool do_derivatives = Moose::doDerivatives(_subproblem, _sys); - - for (unsigned int qp = 0; qp < nqp; qp++) + if (_need_ad_grad_u) { - if (_need_ad_u) - _ad_u[qp] = _ad_zero; - - if (_need_ad_grad_u) + _ad_grad_u.resize(nqp); + for (unsigned int qp = 0; qp < nqp; qp++) _ad_grad_u[qp] = _ad_zero; - if (_need_ad_second_u) - _ad_second_u[qp] = _ad_zero; - - if (_need_ad_u_dot) - _ad_u_dot[qp] = _ad_zero; - - if (_need_ad_u_dotdot) - _ad_u_dotdot[qp] = _ad_zero; - - if (_need_ad_grad_u_dot) - _ad_grad_u_dot[qp] = _ad_zero; + // The latter check here is for handling the fact that we have not yet implemented + // calculation of ad_grad_phi for neighbor and neighbor-face, so if we are in that + // situation we need to default to using the non-ad grad_phi + if (_displaced && _current_ad_grad_phi) + for (unsigned int i = 0; i < num_dofs; i++) + for (unsigned int qp = 0; qp < nqp; qp++) + _ad_grad_u[qp] += _ad_dof_values[i] * (*_current_ad_grad_phi)[i][qp]; + else + for (unsigned int i = 0; i < num_dofs; i++) + for (unsigned int qp = 0; qp < nqp; qp++) + _ad_grad_u[qp] += _ad_dof_values[i] * (*_current_grad_phi)[i][qp]; } - for (unsigned int i = 0; i < num_dofs; i++) + if (_need_ad_second_u) { - _ad_dof_values[i] = (*_sys.currentSolution())(_dof_indices[i]); - - // NOTE! You have to do this AFTER setting the value! - if (do_derivatives) - Moose::derivInsert(_ad_dof_values[i].derivatives(), _dof_indices[i], 1.); + _ad_second_u.resize(nqp); + for (unsigned int qp = 0; qp < nqp; qp++) + _ad_second_u[qp] = _ad_zero; - if (_need_ad_u_dot && _time_integrator && _time_integrator->dt()) - { - _ad_dofs_dot[i] = _ad_dof_values[i]; - _time_integrator->computeADTimeDerivatives(_ad_dofs_dot[i], - _dof_indices[i], - _need_ad_u_dotdot ? _ad_dofs_dotdot[i] - : _ad_real_dummy); - } + for (unsigned int i = 0; i < num_dofs; i++) + for (unsigned int qp = 0; qp < nqp; qp++) + // Note that this will not carry any derivatives with respect to displacements because + // those calculations have not yet been implemented in Assembly + _ad_second_u[qp] += _ad_dof_values[i] * (*_current_second_phi)[i][qp]; } - // Now build up the solution at each quadrature point: - for (unsigned int i = 0; i < num_dofs; i++) + bool is_transient = _subproblem.isTransient(); + if (is_transient) { - for (unsigned int qp = 0; qp < nqp; qp++) + if (_need_ad_u_dot) { - if (_need_ad_u) - _ad_u[qp] += _ad_dof_values[i] * (*_current_phi)[i][qp]; + _ad_dofs_dot.resize(num_dofs); + if (_need_ad_u_dotdot) + _ad_dofs_dotdot.resize(num_dofs); + _ad_u_dot.resize(nqp); + for (unsigned int qp = 0; qp < nqp; qp++) + _ad_u_dot[qp] = _ad_zero; - if (_need_ad_grad_u) + if (_time_integrator && _time_integrator->dt()) { - // The latter check here is for handling the fact that we have not yet implemented - // calculation of ad_grad_phi for neighbor and neighbor-face, so if we are in that - // situation we need to default to using the non-ad grad_phi - if (_displaced && _current_ad_grad_phi) - _ad_grad_u[qp] += _ad_dof_values[i] * (*_current_ad_grad_phi)[i][qp]; - else - _ad_grad_u[qp] += _ad_dof_values[i] * (*_current_grad_phi)[i][qp]; - } - - if (_need_ad_second_u) - // Note that this will not carry any derivatives with respect to displacements because - // those calculations have not yet been implemented in Assembly - _ad_second_u[qp] += _ad_dof_values[i] * (*_current_second_phi)[i][qp]; + for (unsigned int i = 0; i < num_dofs; i++) + _ad_dofs_dot[i] = _ad_dof_values[i]; + for (unsigned int i = 0; i < num_dofs; i++) + _time_integrator->computeADTimeDerivatives(_ad_dofs_dot[i], + _dof_indices[i], + _need_ad_u_dotdot ? _ad_dofs_dotdot[i] + : _ad_real_dummy); - if (_need_ad_u_dot && _time_integrator && _time_integrator->dt()) + for (unsigned int i = 0; i < num_dofs; i++) + for (unsigned int qp = 0; qp < nqp; qp++) + _ad_u_dot[qp] += (*_current_phi)[i][qp] * _ad_dofs_dot[i]; + } + else if (!_time_integrator) { - _ad_u_dot[qp] += (*_current_phi)[i][qp] * _ad_dofs_dot[i]; - if (_need_ad_u_dotdot) - _ad_u_dotdot[qp] += (*_current_phi)[i][qp] * _ad_dofs_dotdot[i]; + for (unsigned int i = 0; i < num_dofs; i++) + _ad_dofs_dot[i] = _dof_values_dot[i]; + for (unsigned int qp = 0; qp < nqp; qp++) + _ad_u_dot[qp] = _u_dot[qp]; } + } + + if (_need_ad_u_dotdot) + { + _ad_u_dotdot.resize(nqp); + for (unsigned int qp = 0; qp < nqp; qp++) + _ad_u_dotdot[qp] = _ad_zero; + + if (_time_integrator && _time_integrator->dt()) + for (unsigned int i = 0; i < num_dofs; i++) + for (unsigned int qp = 0; qp < nqp; qp++) + _ad_u_dotdot[qp] += (*_current_phi)[i][qp] * _ad_dofs_dotdot[i]; + else if (!_time_integrator) + for (unsigned int qp = 0; qp < nqp; qp++) + _ad_u_dotdot[qp] = _u_dotdot[qp]; + } - if (_need_ad_grad_u_dot && _time_integrator && _time_integrator->dt()) + if (_need_ad_grad_u_dot) + { + _ad_grad_u_dot.resize(nqp); + for (unsigned int qp = 0; qp < nqp; qp++) + _ad_grad_u_dot[qp] = _ad_zero; + + if (_time_integrator && _time_integrator->dt()) { // The latter check here is for handling the fact that we have not yet implemented // calculation of ad_grad_phi for neighbor and neighbor-face, so if we are in that // situation we need to default to using the non-ad grad_phi if (_displaced && _current_ad_grad_phi) - _ad_grad_u_dot[qp] += _ad_dofs_dot[i] * (*_current_ad_grad_phi)[i][qp]; + for (unsigned int i = 0; i < num_dofs; i++) + for (unsigned int qp = 0; qp < nqp; qp++) + _ad_grad_u_dot[qp] += _ad_dofs_dot[i] * (*_current_ad_grad_phi)[i][qp]; else - _ad_grad_u_dot[qp] += _ad_dofs_dot[i] * (*_current_grad_phi)[i][qp]; + for (unsigned int i = 0; i < num_dofs; i++) + for (unsigned int qp = 0; qp < nqp; qp++) + _ad_grad_u_dot[qp] += _ad_dofs_dot[i] * (*_current_grad_phi)[i][qp]; } + else if (!_time_integrator) + for (unsigned int qp = 0; qp < nqp; qp++) + _ad_grad_u_dot[qp] = _grad_u_dot[qp]; } } - - if (_need_ad_u_dot && !_time_integrator) - { - for (MooseIndex(nqp) qp = 0; qp < nqp; ++qp) - { - _ad_u_dot[qp] = _u_dot[qp]; - if (_need_ad_u_dotdot) - _ad_u_dotdot[qp] = _u_dotdot[qp]; - } - for (unsigned int i = 0; i < num_dofs; i++) - _ad_dofs_dot[i] = _dof_values_dot[i]; - } - - if (_need_ad_grad_u_dot && !_time_integrator) - for (MooseIndex(nqp) qp = 0; qp < nqp; ++qp) - _ad_grad_u_dot[qp] = _grad_u_dot[qp]; } template <> From dd74a6c7f9b2673fb35d083d82f37d7089bbea0d Mon Sep 17 00:00:00 2001 From: Guillaume Giudicelli Date: Sun, 26 Jan 2025 23:42:51 -0700 Subject: [PATCH 05/18] Remove unused attributes --- framework/include/variables/MooseVariableData.h | 8 ++++---- framework/include/variables/MooseVariableDataFV.h | 8 ++++---- 2 files changed, 8 insertions(+), 8 deletions(-) diff --git a/framework/include/variables/MooseVariableData.h b/framework/include/variables/MooseVariableData.h index 1bab817da8b0..a1f976673a43 100644 --- a/framework/include/variables/MooseVariableData.h +++ b/framework/include/variables/MooseVariableData.h @@ -554,19 +554,19 @@ class MooseVariableData : public MooseVariableDataBase FieldVariableValue _u_dot; /// u_dotdot (second time derivative) - FieldVariableValue _u_dotdot, _u_dotdot_bak; + FieldVariableValue _u_dotdot; /// u_dot_old (time derivative) - FieldVariableValue _u_dot_old, _u_dot_old_bak; + FieldVariableValue _u_dot_old; /// u_dotdot_old (second time derivative) - FieldVariableValue _u_dotdot_old, _u_dotdot_old_bak; + FieldVariableValue _u_dotdot_old; /// derivative of u_dot wrt u VariableValue _du_dot_du; /// derivative of u_dotdot wrt u - VariableValue _du_dotdot_du, _du_dotdot_du_bak; + VariableValue _du_dotdot_du; /// The current qrule. This has to be a reference because the current qrule will be constantly /// changing. If we initialized this to point to one qrule, then in the next calculation we would diff --git a/framework/include/variables/MooseVariableDataFV.h b/framework/include/variables/MooseVariableDataFV.h index a54badf25d1e..126167bcb225 100644 --- a/framework/include/variables/MooseVariableDataFV.h +++ b/framework/include/variables/MooseVariableDataFV.h @@ -373,19 +373,19 @@ class MooseVariableDataFV : public MooseVariableDataBase, public Mes FieldVariableValue _u_dot; /// u_dotdot (second time derivative) - FieldVariableValue _u_dotdot, _u_dotdot_bak; + FieldVariableValue _u_dotdot; /// u_dot_old (time derivative) - FieldVariableValue _u_dot_old, _u_dot_old_bak; + FieldVariableValue _u_dot_old; /// u_dotdot_old (second time derivative) - FieldVariableValue _u_dotdot_old, _u_dotdot_old_bak; + FieldVariableValue _u_dotdot_old; /// derivative of u_dot wrt u VariableValue _du_dot_du; /// derivative of u_dotdot wrt u - VariableValue _du_dotdot_du, _du_dotdot_du_bak; + VariableValue _du_dotdot_du; /// Pointer to time integrator const TimeIntegrator * const _time_integrator; From f2bc7ba9b8ab4e3e7732eabb6389dcef1ffdf93a Mon Sep 17 00:00:00 2001 From: Logan Harbour Date: Tue, 28 Jan 2025 20:44:09 -0700 Subject: [PATCH 06/18] Generalize filling to a lambda and combine RealEigenVector case --- framework/src/variables/MooseVariableData.C | 559 +++++--------------- 1 file changed, 121 insertions(+), 438 deletions(-) diff --git a/framework/src/variables/MooseVariableData.C b/framework/src/variables/MooseVariableData.C index 679a0e8fc562..d902b71644bd 100644 --- a/framework/src/variables/MooseVariableData.C +++ b/framework/src/variables/MooseVariableData.C @@ -430,10 +430,39 @@ MooseVariableData::computeValues() if (num_dofs > 0) fetchDoFValues(); - bool is_transient = _subproblem.isTransient(); - unsigned int nqp = _current_qrule->n_points(); + const bool is_transient = _subproblem.isTransient(); + const unsigned int nqp = _current_qrule->n_points(); auto && active_coupleable_matrix_tags = _subproblem.getActiveFEVariableCoupleableMatrixTags(_tid); + // Map grad_phi using Eigen so that we can perform array operations easier + if constexpr (std::is_same_v) + { + if (_qrule == _current_qrule) + { + _mapped_grad_phi.resize(num_dofs); + for (unsigned int i = 0; i < num_dofs; i++) + { + _mapped_grad_phi[i].resize(nqp, Eigen::Map(nullptr)); + for (unsigned int qp = 0; qp < nqp; qp++) + // Note: this does NOT do any allocation. It is "reconstructing" the object in place + new (&_mapped_grad_phi[i][qp]) + Eigen::Map(const_cast(&(*_current_grad_phi)[i][qp](0))); + } + } + else + { + _mapped_grad_phi_face.resize(num_dofs); + for (unsigned int i = 0; i < num_dofs; i++) + { + _mapped_grad_phi_face[i].resize(nqp, Eigen::Map(nullptr)); + for (unsigned int qp = 0; qp < nqp; qp++) + // Note: this does NOT do any allocation. It is "reconstructing" the object in place + new (&_mapped_grad_phi_face[i][qp]) + Eigen::Map(const_cast(&(*_current_grad_phi)[i][qp](0))); + } + } + } + mooseAssert( !(_need_second || _need_second_old || _need_second_older || _need_second_previous_nl) || _current_second_phi, @@ -443,167 +472,124 @@ MooseVariableData::computeValues() mooseAssert(!(_need_div || _need_div_old) || _current_div_phi, "We're requiring a divergence calculation but have not set a div shape function!"); - // Curl - if (_need_curl) + const auto fill = [this, &nqp, &num_dofs](auto & dest, const auto & phi, const auto & vector) { - _curl_u.resize(nqp); - for (unsigned int i = 0; i < nqp; ++i) - _curl_u[i] = 0.; + constexpr bool is_real = std::is_same_v; + constexpr bool is_real_vector = std::is_same_v; + constexpr bool is_eigen = std::is_same_v; + mooseAssert(is_real || is_real_vector || is_eigen, "Unsupported type"); + + // Output type + using T = typename std::remove_reference_t::value_type; + constexpr bool is_output = std::is_same_v; + constexpr bool is_gradient = std::is_same_v; + constexpr bool is_second = std::is_same_v; + + dest.resize(nqp); + for (unsigned int qp = 0; qp < nqp; ++qp) + { + if constexpr (is_real || is_real_vector) + dest[qp] = 0; + else if constexpr (is_eigen) + { + if constexpr (is_output) + dest[qp].setZero(this->_count); + else if constexpr (is_gradient) + dest[qp].setZero(this->_count, LIBMESH_DIM); + else if constexpr (is_second) + dest[qp].setZero(this->_count, LIBMESH_DIM * LIBMESH_DIM); + else + mooseAssert(false, "Unsupported type"); + } + else + mooseAssert(false, "Unsupported type"); + } + for (unsigned int i = 0; i < num_dofs; i++) for (unsigned int qp = 0; qp < nqp; qp++) - _curl_u[qp] += (*_current_curl_phi)[i][qp] * _vector_tags_dof_u[_solution_tag][i]; - } + { + if constexpr (is_real || is_real_vector || (is_eigen && is_output)) + { + if constexpr (is_output) + dest[qp] += phi[i][qp] * vector[i]; + else if constexpr (is_gradient || is_second) + dest[qp].add_scaled(phi[i][qp], vector[i]); + else + mooseAssert(false, "Unsupported type"); + } + else if constexpr (is_eigen) + { + if constexpr (is_gradient) + { + for (const auto d : make_range(Moose::dim)) + dest[qp].col(d) += phi[i][qp](d) * vector[i]; + } + else if constexpr (is_second) + { + for (unsigned int d = 0, d1 = 0; d1 < LIBMESH_DIM; ++d1) + for (unsigned int d2 = 0; d2 < LIBMESH_DIM; ++d2) + dest[qp].col(d++) += phi[i][qp](d1, d2) * vector[i]; + } + else + mooseAssert(false, "Unsupported type"); + } + } + }; + + // Curl + if (_need_curl) + fill(_curl_u, *_current_curl_phi, _vector_tags_dof_u[_solution_tag]); if (_need_curl_old) - { - _curl_u_old.resize(nqp); - for (unsigned int i = 0; i < nqp; ++i) - _curl_u_old[i] = 0.; - for (unsigned int i = 0; i < num_dofs; i++) - for (unsigned int qp = 0; qp < nqp; qp++) - _curl_u_old[qp] += (*_current_curl_phi)[i][qp] * _vector_tags_dof_u[_old_solution_tag][i]; - } + fill(_curl_u_old, *_current_curl_phi, _vector_tags_dof_u[_old_solution_tag]); // Div if (_need_div) - { - _div_u.resize(nqp); - for (unsigned int i = 0; i < nqp; ++i) - _div_u[i] = 0.; - for (unsigned int i = 0; i < num_dofs; i++) - for (unsigned int qp = 0; qp < nqp; qp++) - _div_u[qp] += (*_current_div_phi)[i][qp] * _vector_tags_dof_u[_solution_tag][i]; - } + fill(_div_u, *_current_div_phi, _vector_tags_dof_u[_solution_tag]); if (is_transient && _need_div_old) - { - _div_u_old.resize(nqp); - for (unsigned int i = 0; i < nqp; ++i) - _div_u_old[i] = 0.; - for (unsigned int i = 0; i < num_dofs; i++) - for (unsigned int qp = 0; qp < nqp; qp++) - _div_u_old[qp] += (*_current_div_phi)[i][qp] * _vector_tags_dof_u[_old_solution_tag][i]; - } + fill(_div_u_old, *_current_div_phi, _vector_tags_dof_u[_old_solution_tag]); // Second if (_need_second) - { - _second_u.resize(nqp); - for (unsigned int i = 0; i < nqp; ++i) - _second_u[i] = 0.; - for (unsigned int i = 0; i < num_dofs; i++) - for (unsigned int qp = 0; qp < nqp; qp++) - _second_u[qp].add_scaled((*_current_second_phi)[i][qp], - _vector_tags_dof_u[_solution_tag][i]); - } + fill(_second_u, *_current_second_phi, _vector_tags_dof_u[_solution_tag]); if (_need_second_previous_nl) - { - _second_u_previous_nl.resize(nqp); - for (unsigned int i = 0; i < nqp; ++i) - _second_u_previous_nl[i] = 0; - for (unsigned int i = 0; i < num_dofs; i++) - for (unsigned int qp = 0; qp < nqp; qp++) - _second_u_previous_nl[qp].add_scaled((*_current_second_phi)[i][qp], - _vector_tags_dof_u[_previous_nl_solution_tag][i]); - } + fill( + _second_u_previous_nl, *_current_second_phi, _vector_tags_dof_u[_previous_nl_solution_tag]); if (is_transient) { if (_need_second_old) - { - _second_u_old.resize(nqp); - for (unsigned int i = 0; i < nqp; ++i) - _second_u_old[i] = 0.; - for (unsigned int i = 0; i < num_dofs; i++) - for (unsigned int qp = 0; qp < nqp; qp++) - _second_u_old[qp].add_scaled((*_current_second_phi)[i][qp], - _vector_tags_dof_u[_old_solution_tag][i]); - } - + fill(_second_u_old, *_current_second_phi, _vector_tags_dof_u[_old_solution_tag]); if (_need_second_older) - { - _second_u_older.resize(nqp); - for (unsigned int i = 0; i < nqp; ++i) - _second_u_older[i] = 0.; - for (unsigned int i = 0; i < num_dofs; i++) - for (unsigned int qp = 0; qp < nqp; qp++) - _second_u_older[qp].add_scaled((*_current_second_phi)[i][qp], - _vector_tags_dof_u[_older_solution_tag][i]); - } + fill(_second_u_older, *_current_second_phi, _vector_tags_dof_u[_older_solution_tag]); } + // Vector tags for (auto tag : _required_vector_tags) { if (_need_vector_tag_u[tag] && _sys.hasVector(tag) && _sys.getVector(tag).closed()) - { - _vector_tag_u[tag].resize(nqp); - for (unsigned int i = 0; i < nqp; ++i) - _vector_tag_u[tag][i] = 0.; - for (unsigned int i = 0; i < num_dofs; i++) - for (unsigned int qp = 0; qp < nqp; qp++) - _vector_tag_u[tag][qp] += (*_current_phi)[i][qp] * _vector_tags_dof_u[tag][i]; - } + fill(_vector_tag_u[tag], *_current_phi, _vector_tags_dof_u[tag]); if (_need_vector_tag_grad[tag] && _sys.hasVector(tag) && _sys.getVector(tag).closed()) - { - _vector_tag_grad[tag].resize(nqp); - for (unsigned int i = 0; i < nqp; ++i) - _vector_tag_grad[tag][i] = 0.; - for (unsigned int i = 0; i < num_dofs; i++) - for (unsigned int qp = 0; qp < nqp; qp++) - _vector_tag_grad[tag][qp].add_scaled((*_current_grad_phi)[i][qp], - _vector_tags_dof_u[tag][i]); - } + fill(_vector_tag_grad[tag], *_current_grad_phi, _vector_tags_dof_u[tag]); } + // Matrix tags for (auto tag : active_coupleable_matrix_tags) if (_need_matrix_tag_u[tag]) - { - _matrix_tag_u[tag].resize(nqp); - for (unsigned int i = 0; i < nqp; ++i) - _matrix_tag_u[tag][i] = 0.; - for (unsigned int i = 0; i < num_dofs; i++) - for (unsigned int qp = 0; qp < nqp; qp++) - _matrix_tag_u[tag][qp] += (*_current_phi)[i][qp] * _matrix_tags_dof_u[tag][i]; - } + fill(_matrix_tag_u[tag], *_current_phi, _matrix_tags_dof_u[tag]); + // Derivatives if (is_transient) { if (_need_u_dot) - { - _u_dot.resize(nqp); - for (unsigned int i = 0; i < nqp; ++i) - _u_dot[i] = 0.; - for (unsigned int i = 0; i < num_dofs; i++) - for (unsigned int qp = 0; qp < nqp; qp++) - _u_dot[qp] += (*_current_phi)[i][qp] * _dof_values_dot[i]; - } + fill(_u_dot, *_current_phi, _dof_values_dot); if (_need_u_dotdot) - { - _u_dotdot.resize(nqp); - for (unsigned int i = 0; i < nqp; ++i) - _u_dotdot[i] = 0.; - for (unsigned int i = 0; i < num_dofs; i++) - for (unsigned int qp = 0; qp < nqp; qp++) - _u_dotdot[qp] += (*_current_phi)[i][qp] * _dof_values_dotdot[i]; - } + fill(_u_dotdot, *_current_phi, _dof_values_dotdot); if (_need_u_dot_old) - { - _u_dot_old.resize(nqp); - for (unsigned int i = 0; i < nqp; ++i) - _u_dot_old[i] = 0.; - for (unsigned int i = 0; i < num_dofs; i++) - for (unsigned int qp = 0; qp < nqp; qp++) - _u_dot_old[qp] += (*_current_phi)[i][qp] * _dof_values_dot_old[i]; - } + fill(_u_dot_old, *_current_phi, _dof_values_dot_old); if (_need_u_dotdot_old) - { - _u_dotdot_old.resize(nqp); - for (unsigned int i = 0; i < nqp; ++i) - _u_dotdot_old[i] = 0.; - for (unsigned int i = 0; i < num_dofs; i++) - for (unsigned int qp = 0; qp < nqp; qp++) - _u_dotdot_old[qp] += (*_current_phi)[i][qp] * _dof_values_dotdot_old[i]; - } + fill(_u_dotdot_old, *_current_phi, _dof_values_dotdot_old); if (_need_du_dot_du) { @@ -625,317 +611,17 @@ MooseVariableData::computeValues() _du_dotdot_du[qp] = _dof_du_dotdot_du[i]; } + // Derivative gradients if (_need_grad_dot) - { - _grad_u_dot.resize(nqp); - for (unsigned int i = 0; i < nqp; ++i) - _grad_u_dot[i] = 0.; - for (unsigned int i = 0; i < num_dofs; i++) - for (unsigned int qp = 0; qp < nqp; qp++) - _grad_u_dot[qp].add_scaled((*_current_grad_phi)[i][qp], _dof_values_dot[i]); - } - + fill(_grad_u_dot, *_current_grad_phi, _dof_values_dot); if (_need_grad_dotdot) - { - _grad_u_dotdot.resize(nqp); - for (unsigned int i = 0; i < nqp; ++i) - _grad_u_dotdot[i] = 0.; - for (unsigned int i = 0; i < num_dofs; i++) - for (unsigned int qp = 0; qp < nqp; qp++) - _grad_u_dotdot[qp].add_scaled((*_current_grad_phi)[i][qp], _dof_values_dotdot[i]); - } + fill(_grad_u_dotdot, *_current_grad_phi, _dof_values_dotdot); } - // Automatic differentiation - if (_need_ad) - computeAD(num_dofs, nqp); -} - -template <> -void -MooseVariableData::computeValues() -{ - unsigned int num_dofs = _dof_indices.size(); - - if (num_dofs > 0) - fetchDoFValues(); - - bool is_transient = _subproblem.isTransient(); - unsigned int nqp = _current_qrule->n_points(); - auto && active_coupleable_matrix_tags = _subproblem.getActiveFEVariableCoupleableMatrixTags(_tid); - - // Map grad_phi using Eigen so that we can perform array operations easier - if (_qrule == _current_qrule) - { - _mapped_grad_phi.resize(num_dofs); - for (unsigned int i = 0; i < num_dofs; i++) - { - _mapped_grad_phi[i].resize(nqp, Eigen::Map(nullptr)); - for (unsigned int qp = 0; qp < nqp; qp++) - // Note: this does NOT do any allocation. It is "reconstructing" the object in place - new (&_mapped_grad_phi[i][qp]) - Eigen::Map(const_cast(&(*_current_grad_phi)[i][qp](0))); - } - } - else - { - _mapped_grad_phi_face.resize(num_dofs); - for (unsigned int i = 0; i < num_dofs; i++) - { - _mapped_grad_phi_face[i].resize(nqp, Eigen::Map(nullptr)); - for (unsigned int qp = 0; qp < nqp; qp++) - // Note: this does NOT do any allocation. It is "reconstructing" the object in place - new (&_mapped_grad_phi_face[i][qp]) - Eigen::Map(const_cast(&(*_current_grad_phi)[i][qp](0))); - } - } - - for (auto tag : _required_vector_tags) - { - if (_need_vector_tag_u[tag]) - _vector_tag_u[tag].resize(nqp); - if (_need_vector_tag_grad[tag]) - _vector_tag_grad[tag].resize(nqp); - } - - for (auto tag : active_coupleable_matrix_tags) - if (_need_matrix_tag_u[tag]) - _matrix_tag_u[tag].resize(nqp); - - if (_need_second) - _second_u.resize(nqp); - - if (_need_curl) - _curl_u.resize(nqp); - - if (_need_div) - _div_u.resize(nqp); - - if (_need_second_previous_nl) - _second_u_previous_nl.resize(nqp); - - if (is_transient) - { - if (_need_u_dot) - _u_dot.resize(nqp); - - if (_need_u_dotdot) - _u_dotdot.resize(nqp); - - if (_need_u_dot_old) - _u_dot_old.resize(nqp); - - if (_need_u_dotdot_old) - _u_dotdot_old.resize(nqp); - - if (_need_du_dot_du) - _du_dot_du.resize(nqp); - - if (_need_du_dotdot_du) - _du_dotdot_du.resize(nqp); - - if (_need_grad_dot) - _grad_u_dot.resize(nqp); - - if (_need_grad_dotdot) - _grad_u_dotdot.resize(nqp); - - if (_need_second_old) - _second_u_old.resize(nqp); - - if (_need_curl_old) - _curl_u_old.resize(nqp); - - if (_need_div_old) - _div_u_old.resize(nqp); - - if (_need_second_older) - _second_u_older.resize(nqp); - } - - for (unsigned int i = 0; i < nqp; ++i) - { - for (auto tag : _required_vector_tags) - { - if (_need_vector_tag_u[tag]) - _vector_tag_u[tag][i].setZero(_count); - if (_need_vector_tag_grad[tag]) - _vector_tag_grad[tag][i].setZero(_count, LIBMESH_DIM); - } - - for (auto tag : active_coupleable_matrix_tags) - if (_need_matrix_tag_u[tag]) - _matrix_tag_u[tag][i].setZero(_count); - - if (_need_second) - _second_u[i].setZero(_count, LIBMESH_DIM * LIBMESH_DIM); - - if (_need_curl) - _curl_u[i].setZero(_count); - - if (_need_div) - _div_u[i].setZero(_count); - - if (_need_second_previous_nl) - _second_u_previous_nl[i].setZero(_count, LIBMESH_DIM * LIBMESH_DIM); - - if (is_transient) - { - if (_need_u_dot) - _u_dot[i].setZero(_count); - - if (_need_u_dotdot) - _u_dotdot[i].setZero(_count); - - if (_need_u_dot_old) - _u_dot_old[i].setZero(_count); - - if (_need_u_dotdot_old) - _u_dotdot_old[i].setZero(_count); - - if (_need_du_dot_du) - _du_dot_du[i] = 0.; - - if (_need_du_dotdot_du) - _du_dotdot_du[i] = 0.; - - if (_need_grad_dot) - _grad_u_dot[i].setZero(_count, LIBMESH_DIM); - - if (_need_grad_dotdot) - _grad_u_dotdot[i].setZero(_count, LIBMESH_DIM); - - if (_need_second_old) - _second_u_old[i].setZero(_count, LIBMESH_DIM * LIBMESH_DIM); - - if (_need_second_older) - _second_u_older[i].setZero(_count, LIBMESH_DIM * LIBMESH_DIM); - - if (_need_curl_old) - _curl_u_old[i].setZero(_count); - - if (_need_div_old) - _div_u_old[i].setZero(_count); - } - } - - bool second_required = - _need_second || _need_second_old || _need_second_older || _need_second_previous_nl; - bool curl_required = _need_curl || _need_curl_old; - bool div_required = _need_div || _need_div_old; - - for (unsigned int i = 0; i < num_dofs; i++) - { - for (unsigned int qp = 0; qp < nqp; qp++) - { - const OutputShape phi_local = (*_current_phi)[i][qp]; - const OutputShapeGradient dphi_qp = (*_current_grad_phi)[i][qp]; - - if (is_transient) - { - if (_need_u_dot) - _u_dot[qp] += phi_local * _dof_values_dot[i]; - - if (_need_u_dotdot) - _u_dotdot[qp] += phi_local * _dof_values_dotdot[i]; - - if (_need_u_dot_old) - _u_dot_old[qp] += phi_local * _dof_values_dot_old[i]; - - if (_need_u_dotdot_old) - _u_dotdot_old[qp] += phi_local * _dof_values_dotdot_old[i]; - - if (_need_grad_dot) - for (const auto d : make_range(Moose::dim)) - _grad_u_dot[qp].col(d) += dphi_qp(d) * _dof_values_dot[i]; - - if (_need_grad_dotdot) - for (const auto d : make_range(Moose::dim)) - _grad_u_dotdot[qp].col(d) += dphi_qp(d) * _dof_values_dotdot[i]; - - if (_need_du_dot_du) - _du_dot_du[qp] = _dof_du_dot_du[i]; - - if (_need_du_dotdot_du) - _du_dotdot_du[qp] = _dof_du_dotdot_du[i]; - } - - if (second_required) - { - mooseAssert( - _current_second_phi, - "We're requiring a second calculation but have not set a second shape function!"); - const RealTensorValue d2phi_local = (*_current_second_phi)[i][qp]; - - if (_need_second) - for (unsigned int d = 0, d1 = 0; d1 < LIBMESH_DIM; ++d1) - for (unsigned int d2 = 0; d2 < LIBMESH_DIM; ++d2) - _second_u[qp].col(d++) += d2phi_local(d1, d2) * _vector_tags_dof_u[_solution_tag][i]; - - if (_need_second_previous_nl) - for (unsigned int d = 0, d1 = 0; d1 < LIBMESH_DIM; ++d1) - for (unsigned int d2 = 0; d2 < LIBMESH_DIM; ++d2) - _second_u_previous_nl[qp].col(d++) += - d2phi_local(d1, d2) * _vector_tags_dof_u[_previous_nl_solution_tag][i]; - - if (is_transient) - { - if (_need_second_old) - for (unsigned int d = 0, d1 = 0; d1 < LIBMESH_DIM; ++d1) - for (unsigned int d2 = 0; d2 < LIBMESH_DIM; ++d2) - _second_u_old[qp].col(d++) += - d2phi_local(d1, d2) * _vector_tags_dof_u[_old_solution_tag][i]; - - if (_need_second_older) - for (unsigned int d = 0, d1 = 0; d1 < LIBMESH_DIM; ++d1) - for (unsigned int d2 = 0; d2 < LIBMESH_DIM; ++d2) - _second_u_older[qp].col(d++) += - d2phi_local(d1, d2) * _vector_tags_dof_u[_older_solution_tag][i]; - } - } - - if (curl_required) - { - mooseAssert(_current_curl_phi, - "We're requiring a curl calculation but have not set a curl shape function!"); - const auto curl_phi_local = (*_current_curl_phi)[i][qp]; - - if (_need_curl) - _curl_u[qp] += curl_phi_local * _vector_tags_dof_u[_solution_tag][i]; - - if (is_transient && _need_curl_old) - _curl_u_old[qp] += curl_phi_local * _vector_tags_dof_u[_old_solution_tag][i]; - } - - if (div_required) - { - mooseAssert(_current_div_phi, - "We're requiring a divergence calculation but have not set a divergence shape " - "function!"); - const auto div_phi_local = (*_current_div_phi)[i][qp]; - - if (_need_div) - _div_u[qp] += div_phi_local * _vector_tags_dof_u[_solution_tag][i]; - - if (is_transient && _need_div_old) - _div_u_old[qp] += div_phi_local * _vector_tags_dof_u[_old_solution_tag][i]; - } - - for (auto tag : _required_vector_tags) - { - if (_need_vector_tag_u[tag]) - _vector_tag_u[tag][qp] += phi_local * _vector_tags_dof_u[tag][i]; - if (_need_vector_tag_grad[tag]) - for (const auto d : make_range(Moose::dim)) - _vector_tag_grad[tag][qp].col(d) += dphi_qp(d) * _vector_tags_dof_u[tag][i]; - } - - for (auto tag : active_coupleable_matrix_tags) - if (_need_matrix_tag_u[tag]) - _matrix_tag_u[tag][qp] += phi_local * _matrix_tags_dof_u[tag][i]; - } - } - // No AD support for array variable yet. + // Automatic differentiation, not for eigen + if constexpr (!std::is_same_v) + if (_need_ad) + computeAD(num_dofs, nqp); } template @@ -1852,10 +1538,7 @@ MooseVariableData::prepare() // FIXME: remove this when the Richard's module is migrated to use the new NodalCoupleable // interface. - if (_dof_indices.size() > 0) - _has_dof_indices = true; - else - _has_dof_indices = false; + _has_dof_indices = _dof_indices.size(); } template From bcf3e2d6827887c0954005fc1eee583e2026bf6c Mon Sep 17 00:00:00 2001 From: Logan Harbour Date: Tue, 28 Jan 2025 21:23:41 -0700 Subject: [PATCH 07/18] Silence unused this --- framework/src/variables/MooseVariableData.C | 3 +++ 1 file changed, 3 insertions(+) diff --git a/framework/src/variables/MooseVariableData.C b/framework/src/variables/MooseVariableData.C index d902b71644bd..823e286a97ac 100644 --- a/framework/src/variables/MooseVariableData.C +++ b/framework/src/variables/MooseVariableData.C @@ -479,6 +479,9 @@ MooseVariableData::computeValues() constexpr bool is_eigen = std::is_same_v; mooseAssert(is_real || is_real_vector || is_eigen, "Unsupported type"); + if constexpr (!is_eigen) + (void)this; + // Output type using T = typename std::remove_reference_t::value_type; constexpr bool is_output = std::is_same_v; From bb047053a183ad16bdcadf10a58387563a1d7361 Mon Sep 17 00:00:00 2001 From: Logan Harbour Date: Wed, 29 Jan 2025 08:30:43 -0700 Subject: [PATCH 08/18] Add fill for divergence --- framework/src/variables/MooseVariableData.C | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/framework/src/variables/MooseVariableData.C b/framework/src/variables/MooseVariableData.C index 823e286a97ac..0e930d6212c4 100644 --- a/framework/src/variables/MooseVariableData.C +++ b/framework/src/variables/MooseVariableData.C @@ -487,6 +487,7 @@ MooseVariableData::computeValues() constexpr bool is_output = std::is_same_v; constexpr bool is_gradient = std::is_same_v; constexpr bool is_second = std::is_same_v; + constexpr bool is_divergence = std::is_same_v; dest.resize(nqp); for (unsigned int qp = 0; qp < nqp; ++qp) @@ -513,7 +514,7 @@ MooseVariableData::computeValues() { if constexpr (is_real || is_real_vector || (is_eigen && is_output)) { - if constexpr (is_output) + if constexpr (is_output || is_divergence) dest[qp] += phi[i][qp] * vector[i]; else if constexpr (is_gradient || is_second) dest[qp].add_scaled(phi[i][qp], vector[i]); From 13ed43d3285935ff17b28fbdb285e661957fd84e Mon Sep 17 00:00:00 2001 From: Logan Harbour Date: Wed, 29 Jan 2025 08:45:35 -0700 Subject: [PATCH 09/18] Make assertions compile time errors, add comments --- framework/src/variables/MooseVariableData.C | 27 ++++++++++++--------- 1 file changed, 16 insertions(+), 11 deletions(-) diff --git a/framework/src/variables/MooseVariableData.C b/framework/src/variables/MooseVariableData.C index 0e930d6212c4..6cbfbdc1fd11 100644 --- a/framework/src/variables/MooseVariableData.C +++ b/framework/src/variables/MooseVariableData.C @@ -472,23 +472,27 @@ MooseVariableData::computeValues() mooseAssert(!(_need_div || _need_div_old) || _current_div_phi, "We're requiring a divergence calculation but have not set a div shape function!"); + // Helper for filling values const auto fill = [this, &nqp, &num_dofs](auto & dest, const auto & phi, const auto & vector) { + // Deduce OutputType constexpr bool is_real = std::is_same_v; constexpr bool is_real_vector = std::is_same_v; constexpr bool is_eigen = std::is_same_v; - mooseAssert(is_real || is_real_vector || is_eigen, "Unsupported type"); + static_assert(is_real || is_real_vector || is_eigen, "Unsupported type"); + // this is only used in the RealEigenVector case to get this->_count if constexpr (!is_eigen) (void)this; - // Output type - using T = typename std::remove_reference_t::value_type; - constexpr bool is_output = std::is_same_v; - constexpr bool is_gradient = std::is_same_v; - constexpr bool is_second = std::is_same_v; - constexpr bool is_divergence = std::is_same_v; + // Deduce type of value within dest MooseArray + using dest_array_type = typename std::remove_reference_t::value_type; + constexpr bool is_output = std::is_same_v; + constexpr bool is_gradient = std::is_same_v; + constexpr bool is_second = std::is_same_v; + constexpr bool is_divergence = std::is_same_v; + // Resize and zero dest.resize(nqp); for (unsigned int qp = 0; qp < nqp; ++qp) { @@ -503,12 +507,13 @@ MooseVariableData::computeValues() else if constexpr (is_second) dest[qp].setZero(this->_count, LIBMESH_DIM * LIBMESH_DIM); else - mooseAssert(false, "Unsupported type"); + static_assert(Moose::always_false, "Unsupported type"); } else - mooseAssert(false, "Unsupported type"); + static_assert(Moose::always_false, "Unsupported type"); } + // Fill for (unsigned int i = 0; i < num_dofs; i++) for (unsigned int qp = 0; qp < nqp; qp++) { @@ -519,7 +524,7 @@ MooseVariableData::computeValues() else if constexpr (is_gradient || is_second) dest[qp].add_scaled(phi[i][qp], vector[i]); else - mooseAssert(false, "Unsupported type"); + static_assert(Moose::always_false, "Unsupported type"); } else if constexpr (is_eigen) { @@ -535,7 +540,7 @@ MooseVariableData::computeValues() dest[qp].col(d++) += phi[i][qp](d1, d2) * vector[i]; } else - mooseAssert(false, "Unsupported type"); + static_assert(Moose::always_false, "Unsupported type"); } } }; From b383c3a8762a14119f4d2bdae2067148c04958f9 Mon Sep 17 00:00:00 2001 From: Logan Harbour Date: Wed, 29 Jan 2025 08:46:55 -0700 Subject: [PATCH 10/18] Remove a branch --- framework/src/variables/MooseVariableData.C | 18 +++++------------- 1 file changed, 5 insertions(+), 13 deletions(-) diff --git a/framework/src/variables/MooseVariableData.C b/framework/src/variables/MooseVariableData.C index 6cbfbdc1fd11..29b338ad5c31 100644 --- a/framework/src/variables/MooseVariableData.C +++ b/framework/src/variables/MooseVariableData.C @@ -563,13 +563,6 @@ MooseVariableData::computeValues() if (_need_second_previous_nl) fill( _second_u_previous_nl, *_current_second_phi, _vector_tags_dof_u[_previous_nl_solution_tag]); - if (is_transient) - { - if (_need_second_old) - fill(_second_u_old, *_current_second_phi, _vector_tags_dof_u[_old_solution_tag]); - if (_need_second_older) - fill(_second_u_older, *_current_second_phi, _vector_tags_dof_u[_older_solution_tag]); - } // Vector tags for (auto tag : _required_vector_tags) @@ -585,18 +578,19 @@ MooseVariableData::computeValues() if (_need_matrix_tag_u[tag]) fill(_matrix_tag_u[tag], *_current_phi, _matrix_tags_dof_u[tag]); - // Derivatives + // Derivatives and old values if (is_transient) { + if (_need_second_old) + fill(_second_u_old, *_current_second_phi, _vector_tags_dof_u[_old_solution_tag]); + if (_need_second_older) + fill(_second_u_older, *_current_second_phi, _vector_tags_dof_u[_older_solution_tag]); if (_need_u_dot) fill(_u_dot, *_current_phi, _dof_values_dot); - if (_need_u_dotdot) fill(_u_dotdot, *_current_phi, _dof_values_dotdot); - if (_need_u_dot_old) fill(_u_dot_old, *_current_phi, _dof_values_dot_old); - if (_need_u_dotdot_old) fill(_u_dotdot_old, *_current_phi, _dof_values_dotdot_old); @@ -609,7 +603,6 @@ MooseVariableData::computeValues() for (unsigned int qp = 0; qp < nqp; qp++) _du_dot_du[qp] = _dof_du_dot_du[i]; } - if (_need_du_dotdot_du) { _du_dotdot_du.resize(nqp); @@ -620,7 +613,6 @@ MooseVariableData::computeValues() _du_dotdot_du[qp] = _dof_du_dotdot_du[i]; } - // Derivative gradients if (_need_grad_dot) fill(_grad_u_dot, *_current_grad_phi, _dof_values_dot); if (_need_grad_dotdot) From abb7b1ee236fa7d6df3c21e32fe9d90aa1842c26 Mon Sep 17 00:00:00 2001 From: Logan Harbour Date: Wed, 29 Jan 2025 08:47:27 -0700 Subject: [PATCH 11/18] Add an extra static assertion --- framework/src/variables/MooseVariableData.C | 2 ++ 1 file changed, 2 insertions(+) diff --git a/framework/src/variables/MooseVariableData.C b/framework/src/variables/MooseVariableData.C index 29b338ad5c31..7f1213a7c0eb 100644 --- a/framework/src/variables/MooseVariableData.C +++ b/framework/src/variables/MooseVariableData.C @@ -542,6 +542,8 @@ MooseVariableData::computeValues() else static_assert(Moose::always_false, "Unsupported type"); } + else + static_assert(Moose::always_false, "Unsupported type"); } }; From 91360aea269cac398cc7356f84aaab0f94698331 Mon Sep 17 00:00:00 2001 From: Logan Harbour Date: Wed, 29 Jan 2025 08:53:26 -0700 Subject: [PATCH 12/18] Correct static assertion --- framework/src/variables/MooseVariableData.C | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/framework/src/variables/MooseVariableData.C b/framework/src/variables/MooseVariableData.C index 7f1213a7c0eb..3f76ba4af1d5 100644 --- a/framework/src/variables/MooseVariableData.C +++ b/framework/src/variables/MooseVariableData.C @@ -507,7 +507,7 @@ MooseVariableData::computeValues() else if constexpr (is_second) dest[qp].setZero(this->_count, LIBMESH_DIM * LIBMESH_DIM); else - static_assert(Moose::always_false, "Unsupported type"); + static_assert(Moose::always_false, "Unsupported type"); } else static_assert(Moose::always_false, "Unsupported type"); @@ -524,7 +524,7 @@ MooseVariableData::computeValues() else if constexpr (is_gradient || is_second) dest[qp].add_scaled(phi[i][qp], vector[i]); else - static_assert(Moose::always_false, "Unsupported type"); + static_assert(Moose::always_false, "Unsupported type"); } else if constexpr (is_eigen) { @@ -540,7 +540,7 @@ MooseVariableData::computeValues() dest[qp].col(d++) += phi[i][qp](d1, d2) * vector[i]; } else - static_assert(Moose::always_false, "Unsupported type"); + static_assert(Moose::always_false, "Unsupported type"); } else static_assert(Moose::always_false, "Unsupported type"); From 8dbabe3f5c46560c67f96c27cd9ac03f93f4132c Mon Sep 17 00:00:00 2001 From: Logan Harbour Date: Wed, 29 Jan 2025 10:27:30 -0700 Subject: [PATCH 13/18] Combine monomial optimization into the same function --- .../include/variables/MooseVariableData.h | 9 + .../variables/MooseVariableConstMonomial.C | 32 +- framework/src/variables/MooseVariableData.C | 283 ++++-------------- 3 files changed, 77 insertions(+), 247 deletions(-) diff --git a/framework/include/variables/MooseVariableData.h b/framework/include/variables/MooseVariableData.h index a1f976673a43..2db0a1cfcba8 100644 --- a/framework/include/variables/MooseVariableData.h +++ b/framework/include/variables/MooseVariableData.h @@ -461,6 +461,15 @@ class MooseVariableData : public MooseVariableDataBase void assignADNodalValue(const ADReal & value, const unsigned int & component); void fetchADNodalValues(); + /** + * Internal method for computeValues() and computeMonomialValues() + * + * Monomial is a template parameter so that we get compile time optimization + * for monomial vs non-monomial + */ + template + void computeValuesInternal(); + const libMesh::FEType & _fe_type; const unsigned int _var_num; diff --git a/framework/src/variables/MooseVariableConstMonomial.C b/framework/src/variables/MooseVariableConstMonomial.C index d797470d3b04..93f6567a23d1 100644 --- a/framework/src/variables/MooseVariableConstMonomial.C +++ b/framework/src/variables/MooseVariableConstMonomial.C @@ -42,50 +42,26 @@ void MooseVariableConstMonomial::computeElemValues() { _element_data->setGeometry(Moose::Volume); - - // We have not optimized AD calculations for const monomials yet, so we fall back on the - // non-optimized routine - if (_element_data->needsAD()) - _element_data->computeValues(); - else - _element_data->computeMonomialValues(); + _element_data->computeMonomialValues(); } void MooseVariableConstMonomial::computeElemValuesFace() { _element_data->setGeometry(Moose::Face); - - // We have not optimized AD calculations for const monomials yet, so we fall back on the - // non-optimized routine - if (_element_data->needsAD()) - _element_data->computeValues(); - else - _element_data->computeMonomialValues(); + _element_data->computeMonomialValues(); } void MooseVariableConstMonomial::computeNeighborValues() { _neighbor_data->setGeometry(Moose::Volume); - - // We have not optimized AD calculations for const monomials yet, so we fall back on the - // non-optimized routine - if (_neighbor_data->needsAD()) - _neighbor_data->computeValues(); - else - _neighbor_data->computeMonomialValues(); + _neighbor_data->computeMonomialValues(); } void MooseVariableConstMonomial::computeNeighborValuesFace() { _neighbor_data->setGeometry(Moose::Face); - - // We have not optimized AD calculations for const monomials yet, so we fall back on the - // non-optimized routine - if (_neighbor_data->needsAD()) - _neighbor_data->computeValues(); - else - _neighbor_data->computeMonomialValues(); + _neighbor_data->computeMonomialValues(); } diff --git a/framework/src/variables/MooseVariableData.C b/framework/src/variables/MooseVariableData.C index 3f76ba4af1d5..c2f67be5385a 100644 --- a/framework/src/variables/MooseVariableData.C +++ b/framework/src/variables/MooseVariableData.C @@ -422,11 +422,11 @@ MooseVariableData::divPhiFace() const } template +template void -MooseVariableData::computeValues() +MooseVariableData::computeValuesInternal() { - unsigned int num_dofs = _dof_indices.size(); - + const unsigned int num_dofs = _dof_indices.size(); if (num_dofs > 0) fetchDoFValues(); @@ -492,10 +492,12 @@ MooseVariableData::computeValues() constexpr bool is_second = std::is_same_v; constexpr bool is_divergence = std::is_same_v; - // Resize and zero - dest.resize(nqp); - for (unsigned int qp = 0; qp < nqp; ++qp) + // Sets a value to zero at a quadrature point + const auto set_zero = [this, &dest](const auto qp) { + if constexpr (!is_eigen) + (void)this; + if constexpr (is_real || is_real_vector) dest[qp] = 0; else if constexpr (is_eigen) @@ -511,40 +513,60 @@ MooseVariableData::computeValues() } else static_assert(Moose::always_false, "Unsupported type"); - } + }; - // Fill - for (unsigned int i = 0; i < num_dofs; i++) - for (unsigned int qp = 0; qp < nqp; qp++) + // Accumulates a value + const auto accumulate = [this, &dest, &phi, &vector](const auto i, const auto qp) + { + if constexpr (is_real || is_real_vector || (is_eigen && is_output)) { - if constexpr (is_real || is_real_vector || (is_eigen && is_output)) + if constexpr (is_output || is_divergence) + dest[qp] += phi[i][qp] * vector[i]; + else if constexpr (is_gradient || is_second) + dest[qp].add_scaled(phi[i][qp], vector[i]); + else + static_assert(Moose::always_false, "Unsupported type"); + } + else if constexpr (is_eigen) + { + if constexpr (is_gradient) { - if constexpr (is_output || is_divergence) - dest[qp] += phi[i][qp] * vector[i]; - else if constexpr (is_gradient || is_second) - dest[qp].add_scaled(phi[i][qp], vector[i]); - else - static_assert(Moose::always_false, "Unsupported type"); + for (const auto d : make_range(Moose::dim)) + dest[qp].col(d) += phi[i][qp](d) * vector[i]; } - else if constexpr (is_eigen) + else if constexpr (is_second) { - if constexpr (is_gradient) - { - for (const auto d : make_range(Moose::dim)) - dest[qp].col(d) += phi[i][qp](d) * vector[i]; - } - else if constexpr (is_second) - { - for (unsigned int d = 0, d1 = 0; d1 < LIBMESH_DIM; ++d1) - for (unsigned int d2 = 0; d2 < LIBMESH_DIM; ++d2) - dest[qp].col(d++) += phi[i][qp](d1, d2) * vector[i]; - } - else - static_assert(Moose::always_false, "Unsupported type"); + for (unsigned int d = 0, d1 = 0; d1 < LIBMESH_DIM; ++d1) + for (unsigned int d2 = 0; d2 < LIBMESH_DIM; ++d2) + dest[qp].col(d++) += phi[i][qp](d1, d2) * vector[i]; } else - static_assert(Moose::always_false, "Unsupported type"); + static_assert(Moose::always_false, "Unsupported type"); } + else + static_assert(Moose::always_false, "Unsupported type"); + }; + + dest.resize(nqp); + + // Monomial case, accumulate dest[0] and set dest[>0] to dest[0] + if constexpr (monomial) + { + mooseAssert(num_dofs == 1, "Should have only one dof"); + set_zero(0); + accumulate(0, 0); + for (unsigned int qp = 0; qp < nqp; ++qp) + dest[qp] = dest[0]; + } + // Non monomial case + else + { + for (unsigned int qp = 0; qp < nqp; ++qp) + set_zero(qp); + for (unsigned int i = 0; i < num_dofs; i++) + for (unsigned int qp = 0; qp < nqp; qp++) + accumulate(i, qp); + } }; // Curl @@ -627,6 +649,13 @@ MooseVariableData::computeValues() computeAD(num_dofs, nqp); } +template +void +MooseVariableData::computeValues() +{ + computeValuesInternal(); +} + template void MooseVariableData::computeMonomialValues() @@ -634,195 +663,11 @@ MooseVariableData::computeMonomialValues() if (_dof_indices.size() == 0) return; + // Monomial optimizations are not appropriate after p-refinement if (_elem->p_level()) - { - // The optimizations in this routine are not appropriate after p-refinement computeValues(); - return; - } - - bool is_transient = _subproblem.isTransient(); - unsigned int nqp = _current_qrule->n_points(); - - if (_need_second) - _second_u.resize(nqp); - - if (_need_second_previous_nl) - _second_u_previous_nl.resize(nqp); - - if (is_transient) - { - if (_need_u_dot) - _u_dot.resize(nqp); - - if (_need_u_dotdot) - _u_dotdot.resize(nqp); - - if (_need_u_dot_old) - _u_dot_old.resize(nqp); - - if (_need_u_dotdot_old) - _u_dotdot_old.resize(nqp); - - if (_need_du_dot_du) - _du_dot_du.resize(nqp); - - if (_need_du_dotdot_du) - _du_dotdot_du.resize(nqp); - - if (_need_second_old) - _second_u_old.resize(nqp); - - if (_need_second_older) - _second_u_older.resize(nqp); - } - - if (is_transient) - { - if (_need_dof_values_dot) - _dof_values_dot.resize(1); - if (_need_dof_values_dotdot) - _dof_values_dotdot.resize(1); - if (_need_dof_values_dot_old) - _dof_values_dot_old.resize(1); - if (_need_dof_values_dotdot_old) - _dof_values_dotdot_old.resize(1); - } - - const dof_id_type & idx = _dof_indices[0]; - Real u_dot = 0.; - Real u_dotdot = 0.; - Real u_dot_old = 0.; - Real u_dotdot_old = 0.; - const Real & du_dot_du = _sys.duDotDu(_var_num); - const Real & du_dotdot_du = _sys.duDotDotDu(); - - if (is_transient) - { - if (_sys.solutionUDot()) - u_dot = (*_sys.solutionUDot())(idx); - if (_sys.solutionUDotDot()) - u_dotdot = (*_sys.solutionUDotDot())(idx); - if (_sys.solutionUDotOld()) - u_dot_old = (*_sys.solutionUDotOld())(idx); - if (_sys.solutionUDotDotOld()) - u_dotdot_old = (*_sys.solutionUDotDotOld())(idx); - - if (_need_dof_values_dot) - _dof_values_dot[0] = u_dot; - - if (_need_dof_values_dotdot) - _dof_values_dotdot[0] = u_dotdot; - } - - auto phi = (*_current_phi)[0][0]; - - if (is_transient) - { - if (_need_u_dot) - _u_dot[0] = phi * u_dot; - - if (_need_u_dotdot) - _u_dotdot[0] = phi * u_dotdot; - - if (_need_u_dot_old) - _u_dot_old[0] = phi * u_dot_old; - - if (_need_u_dotdot_old) - _u_dotdot_old[0] = phi * u_dotdot_old; - - if (_need_du_dot_du) - _du_dot_du[0] = du_dot_du; - - if (_need_du_dotdot_du) - _du_dotdot_du[0] = du_dotdot_du; - } - - for (unsigned qp = 1; qp < nqp; ++qp) - { - if (is_transient) - { - if (_need_u_dot) - _u_dot[qp] = _u_dot[0]; - - if (_need_u_dotdot) - _u_dotdot[qp] = _u_dotdot[0]; - - if (_need_u_dot_old) - _u_dot_old[qp] = _u_dot_old[0]; - - if (_need_u_dotdot_old) - _u_dotdot_old[qp] = _u_dotdot_old[0]; - - if (_need_du_dot_du) - _du_dot_du[qp] = _du_dot_du[0]; - - if (_need_du_dotdot_du) - _du_dotdot_du[qp] = _du_dotdot_du[0]; - } - } - - auto && active_coupleable_matrix_tags = _subproblem.getActiveFEVariableCoupleableMatrixTags(_tid); - - for (auto tag : _required_vector_tags) - { - if (_need_vector_tag_u[tag] || _need_vector_tag_grad[tag] || _need_vector_tag_dof_u[tag]) - if ((_subproblem.vectorTagType(tag) == Moose::VECTOR_TAG_RESIDUAL && - _subproblem.safeAccessTaggedVectors()) || - _subproblem.vectorTagType(tag) == Moose::VECTOR_TAG_SOLUTION) - // tag is defined on problem but may not be used by a system - if (_sys.hasVector(tag) && _sys.getVector(tag).closed()) - { - auto & vec = _sys.getVector(tag); - _vector_tags_dof_u[tag].resize(1); - _vector_tags_dof_u[tag][0] = vec(_dof_indices[0]); - } - - if (_need_vector_tag_u[tag]) - { - _vector_tag_u[tag].resize(nqp); - auto v = phi * _vector_tags_dof_u[tag][0]; - for (unsigned int qp = 0; qp < nqp; ++qp) - _vector_tag_u[tag][qp] = v; - } - if (_need_vector_tag_grad[tag]) - _vector_tag_grad[tag].resize(nqp); - } - - if (_subproblem.safeAccessTaggedMatrices()) - { - auto & active_coupleable_matrix_tags = - _subproblem.getActiveFEVariableCoupleableMatrixTags(_tid); - for (auto tag : active_coupleable_matrix_tags) - { - _matrix_tags_dof_u[tag].resize(1); - if (_need_matrix_tag_dof_u[tag] || _need_matrix_tag_u[tag]) - if (_sys.hasMatrix(tag) && _sys.matrixTagActive(tag) && _sys.getMatrix(tag).closed()) - { - auto & mat = _sys.getMatrix(tag); - { - Threads::spin_mutex::scoped_lock lock(Threads::spin_mtx); - _matrix_tags_dof_u[tag][0] = mat(_dof_indices[0], _dof_indices[0]); - } - } - } - } - for (auto tag : active_coupleable_matrix_tags) - if (_need_matrix_tag_u[tag]) - { - _matrix_tag_u[tag].resize(nqp); - auto v = phi * _matrix_tags_dof_u[tag][0]; - for (unsigned int qp = 0; qp < nqp; ++qp) - _matrix_tag_u[tag][qp] = v; - } -} - -template <> -void -MooseVariableData::computeMonomialValues() -{ - // FIXME: will think of optimization later - computeValues(); + else + computeValuesInternal(); } template From 55c1632c88755b4e35f9dbe62a6f9ca4cf610368 Mon Sep 17 00:00:00 2001 From: Logan Harbour Date: Wed, 29 Jan 2025 14:17:39 -0700 Subject: [PATCH 14/18] Fix ignores, remove unused captures --- framework/src/variables/MooseVariableData.C | 9 ++++++--- 1 file changed, 6 insertions(+), 3 deletions(-) diff --git a/framework/src/variables/MooseVariableData.C b/framework/src/variables/MooseVariableData.C index c2f67be5385a..429b511d9d78 100644 --- a/framework/src/variables/MooseVariableData.C +++ b/framework/src/variables/MooseVariableData.C @@ -475,6 +475,9 @@ MooseVariableData::computeValuesInternal() // Helper for filling values const auto fill = [this, &nqp, &num_dofs](auto & dest, const auto & phi, const auto & vector) { + if constexpr (monomial) + libmesh_ignore(num_dofs); + // Deduce OutputType constexpr bool is_real = std::is_same_v; constexpr bool is_real_vector = std::is_same_v; @@ -483,7 +486,7 @@ MooseVariableData::computeValuesInternal() // this is only used in the RealEigenVector case to get this->_count if constexpr (!is_eigen) - (void)this; + libmesh_ignore(this); // Deduce type of value within dest MooseArray using dest_array_type = typename std::remove_reference_t::value_type; @@ -496,7 +499,7 @@ MooseVariableData::computeValuesInternal() const auto set_zero = [this, &dest](const auto qp) { if constexpr (!is_eigen) - (void)this; + libmesh_ignore(this); if constexpr (is_real || is_real_vector) dest[qp] = 0; @@ -516,7 +519,7 @@ MooseVariableData::computeValuesInternal() }; // Accumulates a value - const auto accumulate = [this, &dest, &phi, &vector](const auto i, const auto qp) + const auto accumulate = [&dest, &phi, &vector](const auto i, const auto qp) { if constexpr (is_real || is_real_vector || (is_eigen && is_output)) { From f1052da95e818910bf6bb89ab658fdb94aeeb5f8 Mon Sep 17 00:00:00 2001 From: Logan Harbour Date: Wed, 29 Jan 2025 14:18:12 -0700 Subject: [PATCH 15/18] Rename vector to dof_values --- framework/src/variables/MooseVariableData.C | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/framework/src/variables/MooseVariableData.C b/framework/src/variables/MooseVariableData.C index 429b511d9d78..282f8ccf5f16 100644 --- a/framework/src/variables/MooseVariableData.C +++ b/framework/src/variables/MooseVariableData.C @@ -473,7 +473,7 @@ MooseVariableData::computeValuesInternal() "We're requiring a divergence calculation but have not set a div shape function!"); // Helper for filling values - const auto fill = [this, &nqp, &num_dofs](auto & dest, const auto & phi, const auto & vector) + const auto fill = [this, &nqp, &num_dofs](auto & dest, const auto & phi, const auto & dof_values) { if constexpr (monomial) libmesh_ignore(num_dofs); @@ -519,14 +519,14 @@ MooseVariableData::computeValuesInternal() }; // Accumulates a value - const auto accumulate = [&dest, &phi, &vector](const auto i, const auto qp) + const auto accumulate = [&dest, &phi, &dof_values](const auto i, const auto qp) { if constexpr (is_real || is_real_vector || (is_eigen && is_output)) { if constexpr (is_output || is_divergence) - dest[qp] += phi[i][qp] * vector[i]; + dest[qp] += phi[i][qp] * dof_values[i]; else if constexpr (is_gradient || is_second) - dest[qp].add_scaled(phi[i][qp], vector[i]); + dest[qp].add_scaled(phi[i][qp], dof_values[i]); else static_assert(Moose::always_false, "Unsupported type"); } @@ -535,13 +535,13 @@ MooseVariableData::computeValuesInternal() if constexpr (is_gradient) { for (const auto d : make_range(Moose::dim)) - dest[qp].col(d) += phi[i][qp](d) * vector[i]; + dest[qp].col(d) += phi[i][qp](d) * dof_values[i]; } else if constexpr (is_second) { for (unsigned int d = 0, d1 = 0; d1 < LIBMESH_DIM; ++d1) for (unsigned int d2 = 0; d2 < LIBMESH_DIM; ++d2) - dest[qp].col(d++) += phi[i][qp](d1, d2) * vector[i]; + dest[qp].col(d++) += phi[i][qp](d1, d2) * dof_values[i]; } else static_assert(Moose::always_false, "Unsupported type"); From bb074dd0054b2adc8e264849eed00bc8af849a73 Mon Sep 17 00:00:00 2001 From: Logan Harbour Date: Wed, 29 Jan 2025 16:29:35 -0700 Subject: [PATCH 16/18] Add extra assertion --- framework/src/variables/MooseVariableData.C | 2 ++ 1 file changed, 2 insertions(+) diff --git a/framework/src/variables/MooseVariableData.C b/framework/src/variables/MooseVariableData.C index 282f8ccf5f16..dafed22e177b 100644 --- a/framework/src/variables/MooseVariableData.C +++ b/framework/src/variables/MooseVariableData.C @@ -494,6 +494,8 @@ MooseVariableData::computeValuesInternal() constexpr bool is_gradient = std::is_same_v; constexpr bool is_second = std::is_same_v; constexpr bool is_divergence = std::is_same_v; + static_assert(is_output || is_gradient || is_second || is_divergence, + "Unsupported destination array type"); // Sets a value to zero at a quadrature point const auto set_zero = [this, &dest](const auto qp) From 9895d9d810c9da6b2953eba22c04c2fcebae4610 Mon Sep 17 00:00:00 2001 From: Logan Harbour Date: Wed, 29 Jan 2025 16:35:26 -0700 Subject: [PATCH 17/18] Expand always_false to allow multiple types, fix clang min bugs --- framework/include/utils/MooseTypes.h | 2 +- framework/src/variables/MooseVariableData.C | 10 +++++----- 2 files changed, 6 insertions(+), 6 deletions(-) diff --git a/framework/include/utils/MooseTypes.h b/framework/include/utils/MooseTypes.h index 2ae52824349b..468949505849 100644 --- a/framework/include/utils/MooseTypes.h +++ b/framework/include/utils/MooseTypes.h @@ -581,7 +581,7 @@ struct IsADType> * error with constexpr-based if conditions. The templating delays the triggering * of the static assertion until the template is instantiated. */ -template +template constexpr std::false_type always_false{}; } // namespace Moose diff --git a/framework/src/variables/MooseVariableData.C b/framework/src/variables/MooseVariableData.C index dafed22e177b..0bdbcf781564 100644 --- a/framework/src/variables/MooseVariableData.C +++ b/framework/src/variables/MooseVariableData.C @@ -514,10 +514,10 @@ MooseVariableData::computeValuesInternal() else if constexpr (is_second) dest[qp].setZero(this->_count, LIBMESH_DIM * LIBMESH_DIM); else - static_assert(Moose::always_false, "Unsupported type"); + static_assert(Moose::always_false, "Unsupported type"); } else - static_assert(Moose::always_false, "Unsupported type"); + static_assert(Moose::always_false, "Unsupported type"); }; // Accumulates a value @@ -530,7 +530,7 @@ MooseVariableData::computeValuesInternal() else if constexpr (is_gradient || is_second) dest[qp].add_scaled(phi[i][qp], dof_values[i]); else - static_assert(Moose::always_false, "Unsupported type"); + static_assert(Moose::always_false, "Unsupported type"); } else if constexpr (is_eigen) { @@ -546,10 +546,10 @@ MooseVariableData::computeValuesInternal() dest[qp].col(d++) += phi[i][qp](d1, d2) * dof_values[i]; } else - static_assert(Moose::always_false, "Unsupported type"); + static_assert(Moose::always_false, "Unsupported type"); } else - static_assert(Moose::always_false, "Unsupported type"); + static_assert(Moose::always_false, "Unsupported type"); }; dest.resize(nqp); From d7b226f00660fcfa5fd390ffc1f503731ae92c83 Mon Sep 17 00:00:00 2001 From: Logan Harbour Date: Wed, 29 Jan 2025 16:40:12 -0700 Subject: [PATCH 18/18] Only set starting with the 1st qp --- framework/src/variables/MooseVariableData.C | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/framework/src/variables/MooseVariableData.C b/framework/src/variables/MooseVariableData.C index 0bdbcf781564..6d4f7c393d62 100644 --- a/framework/src/variables/MooseVariableData.C +++ b/framework/src/variables/MooseVariableData.C @@ -560,7 +560,7 @@ MooseVariableData::computeValuesInternal() mooseAssert(num_dofs == 1, "Should have only one dof"); set_zero(0); accumulate(0, 0); - for (unsigned int qp = 0; qp < nqp; ++qp) + for (unsigned int qp = 1; qp < nqp; ++qp) dest[qp] = dest[0]; } // Non monomial case