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 "SecantSolve.h" 11 : 12 : #include "Executioner.h" 13 : #include "FEProblemBase.h" 14 : #include "NonlinearSystem.h" 15 : #include "AllLocalDofIndicesThread.h" 16 : #include "Console.h" 17 : 18 : InputParameters 19 0 : SecantSolve::validParams() 20 : { 21 0 : InputParameters params = FixedPointSolve::validParams(); 22 : 23 0 : return params; 24 : } 25 : 26 576 : SecantSolve::SecantSolve(Executioner & ex) : FixedPointSolve(ex) 27 : { 28 576 : allocateStorage(true); 29 : 30 576 : _transformed_pps_values.resize(_transformed_pps.size()); 31 676 : for (size_t i = 0; i < _transformed_pps.size(); i++) 32 100 : _transformed_pps_values[i].resize(4); 33 576 : _secondary_transformed_pps_values.resize(_secondary_transformed_pps.size()); 34 644 : for (size_t i = 0; i < _secondary_transformed_pps.size(); i++) 35 68 : _secondary_transformed_pps_values[i].resize(4); 36 576 : } 37 : 38 : void 39 864 : SecantSolve::allocateStorage(const bool primary) 40 : { 41 : TagID fxn_m1_tagid; 42 : TagID xn_m1_tagid; 43 : TagID fxn_m2_tagid; 44 : TagID xn_m2_tagid; 45 864 : if (primary) 46 : { 47 576 : xn_m1_tagid = _problem.addVectorTag("xn_m1", Moose::VECTOR_TAG_SOLUTION); 48 576 : fxn_m1_tagid = _problem.addVectorTag("fxn_m1", Moose::VECTOR_TAG_SOLUTION); 49 576 : xn_m2_tagid = _problem.addVectorTag("xn_m2", Moose::VECTOR_TAG_SOLUTION); 50 576 : fxn_m2_tagid = _problem.addVectorTag("fxn_m2", Moose::VECTOR_TAG_SOLUTION); 51 576 : _xn_m1_tagid = xn_m1_tagid; 52 576 : _fxn_m1_tagid = fxn_m1_tagid; 53 576 : _xn_m2_tagid = xn_m2_tagid; 54 576 : _fxn_m2_tagid = fxn_m2_tagid; 55 : } 56 : else 57 : { 58 288 : xn_m1_tagid = _problem.addVectorTag("secondary_xn_m1", Moose::VECTOR_TAG_SOLUTION); 59 288 : fxn_m1_tagid = _problem.addVectorTag("secondary_fxn_m1", Moose::VECTOR_TAG_SOLUTION); 60 288 : xn_m2_tagid = _problem.addVectorTag("secondary_xn_m2", Moose::VECTOR_TAG_SOLUTION); 61 288 : fxn_m2_tagid = _problem.addVectorTag("secondary_fxn_m2", Moose::VECTOR_TAG_SOLUTION); 62 288 : _secondary_xn_m1_tagid = xn_m1_tagid; 63 288 : _secondary_fxn_m1_tagid = fxn_m1_tagid; 64 288 : _secondary_xn_m2_tagid = xn_m2_tagid; 65 288 : _secondary_fxn_m2_tagid = fxn_m2_tagid; 66 : } 67 : 68 : // TODO: We would only need to store the solution for the degrees of freedom that 69 : // will be transformed, not the entire solution. 70 : // Store solution vectors for the two previous points and their evaluation 71 864 : _solver_sys.addVector(xn_m1_tagid, false, PARALLEL); 72 864 : _solver_sys.addVector(fxn_m1_tagid, false, PARALLEL); 73 864 : _solver_sys.addVector(xn_m2_tagid, false, PARALLEL); 74 864 : _solver_sys.addVector(fxn_m2_tagid, false, PARALLEL); 75 864 : } 76 : 77 : void 78 42852 : SecantSolve::saveVariableValues(const bool primary) 79 : { 80 : TagID fxn_m1_tagid; 81 : TagID xn_m1_tagid; 82 : TagID fxn_m2_tagid; 83 : TagID xn_m2_tagid; 84 42852 : if (primary) 85 : { 86 30592 : fxn_m1_tagid = _fxn_m1_tagid; 87 30592 : xn_m1_tagid = _xn_m1_tagid; 88 30592 : fxn_m2_tagid = _fxn_m2_tagid; 89 30592 : xn_m2_tagid = _xn_m2_tagid; 90 : } 91 : else 92 : { 93 12260 : fxn_m1_tagid = _secondary_fxn_m1_tagid; 94 12260 : xn_m1_tagid = _secondary_xn_m1_tagid; 95 12260 : fxn_m2_tagid = _secondary_fxn_m2_tagid; 96 12260 : xn_m2_tagid = _secondary_xn_m2_tagid; 97 : } 98 : 99 : // Save previous variable values 100 42852 : NumericVector<Number> & solution = _solver_sys.solution(); 101 42852 : NumericVector<Number> & fxn_m1 = _solver_sys.getVector(fxn_m1_tagid); 102 42852 : NumericVector<Number> & xn_m1 = _solver_sys.getVector(xn_m1_tagid); 103 42852 : NumericVector<Number> & fxn_m2 = _solver_sys.getVector(fxn_m2_tagid); 104 42852 : NumericVector<Number> & xn_m2 = _solver_sys.getVector(xn_m2_tagid); 105 : 106 : // Advance one step 107 42852 : xn_m2 = xn_m1; 108 42852 : fxn_m2 = fxn_m1; 109 : 110 : // Before a solve, solution is a sequence term, after a solve, solution is the evaluated term 111 42852 : xn_m1 = solution; 112 42852 : } 113 : 114 : void 115 42852 : SecantSolve::savePostprocessorValues(const bool primary) 116 : { 117 : const std::vector<PostprocessorName> * transformed_pps; 118 : std::vector<std::vector<PostprocessorValue>> * transformed_pps_values; 119 42852 : if (primary) 120 : { 121 30592 : transformed_pps = &_transformed_pps; 122 30592 : transformed_pps_values = &_transformed_pps_values; 123 : } 124 : else 125 : { 126 12260 : transformed_pps = &_secondary_transformed_pps; 127 12260 : transformed_pps_values = &_secondary_transformed_pps_values; 128 : } 129 : 130 : // Save previous postprocessor values 131 49652 : for (size_t i = 0; i < (*transformed_pps).size(); i++) 132 : { 133 : // Advance one step 134 6800 : (*transformed_pps_values)[i][3] = (*transformed_pps_values)[i][1]; 135 6800 : (*transformed_pps_values)[i][2] = (*transformed_pps_values)[i][0]; 136 : 137 : // Save current value 138 : // Primary: this is done before the timestep's solves and before timestep_begin transfers, 139 : // so the value is the result of the previous Secant update (xn_m1) 140 : // Secondary: this is done after the secondary solve, but before timestep_end postprocessors 141 : // are computed, or timestep_end transfers are received. 142 : // This value is the same as before the solve (xn_m1) 143 6800 : (*transformed_pps_values)[i][1] = getPostprocessorValueByName((*transformed_pps)[i]); 144 : } 145 42852 : } 146 : 147 : bool 148 15831 : SecantSolve::useFixedPointAlgorithmUpdateInsteadOfPicard(const bool primary) 149 : { 150 : // Need at least two evaluations to compute the Secant slope 151 15831 : if (primary) 152 9730 : return _fixed_point_it > 1; 153 : else 154 6101 : return _main_fixed_point_it > 1; 155 : } 156 : 157 : void 158 5030 : SecantSolve::transformPostprocessors(const bool primary) 159 : { 160 : Real relaxation_factor; 161 : const std::vector<PostprocessorName> * transformed_pps; 162 : std::vector<std::vector<PostprocessorValue>> * transformed_pps_values; 163 5030 : if (primary) 164 : { 165 2940 : relaxation_factor = _relax_factor; 166 2940 : transformed_pps = &_transformed_pps; 167 2940 : transformed_pps_values = &_transformed_pps_values; 168 : } 169 : else 170 : { 171 2090 : relaxation_factor = _secondary_relaxation_factor; 172 2090 : transformed_pps = &_secondary_transformed_pps; 173 2090 : transformed_pps_values = &_secondary_transformed_pps_values; 174 : } 175 : 176 : // Relax postprocessors for the main application 177 10060 : for (size_t i = 0; i < (*transformed_pps).size(); i++) 178 : { 179 : // Get new postprocessor value 180 5030 : const Real fxn_m1 = getPostprocessorValueByName((*transformed_pps)[i]); 181 5030 : const Real xn_m1 = (*transformed_pps_values)[i][1]; 182 5030 : const Real fxn_m2 = (*transformed_pps_values)[i][2]; 183 5030 : const Real xn_m2 = (*transformed_pps_values)[i][3]; 184 : 185 : // Save fxn_m1, received or computed before the solve 186 5030 : (*transformed_pps_values)[i][0] = fxn_m1; 187 : 188 : // Compute and set relaxed value 189 5030 : Real new_value = fxn_m1; 190 5030 : if (!MooseUtils::absoluteFuzzyEqual(fxn_m1 - xn_m1 - fxn_m2 + xn_m2, 0)) 191 5030 : new_value = xn_m1 - (fxn_m1 - xn_m1) * (xn_m1 - xn_m2) / (fxn_m1 - xn_m1 - fxn_m2 + xn_m2); 192 : 193 : // Relax update if desired 194 5030 : new_value = relaxation_factor * new_value + (1 - relaxation_factor) * xn_m1; 195 : 196 5030 : _problem.setPostprocessorValueByName((*transformed_pps)[i], new_value); 197 : } 198 5030 : } 199 : 200 : void 201 5239 : SecantSolve::transformVariables(const std::set<dof_id_type> & target_dofs, const bool primary) 202 : { 203 : Real relaxation_factor; 204 : TagID fxn_m1_tagid; 205 : TagID xn_m1_tagid; 206 : TagID fxn_m2_tagid; 207 : TagID xn_m2_tagid; 208 5239 : if (primary) 209 : { 210 3688 : relaxation_factor = _relax_factor; 211 3688 : fxn_m1_tagid = _fxn_m1_tagid; 212 3688 : xn_m1_tagid = _xn_m1_tagid; 213 3688 : fxn_m2_tagid = _fxn_m2_tagid; 214 3688 : xn_m2_tagid = _xn_m2_tagid; 215 : } 216 : else 217 : { 218 1551 : relaxation_factor = _secondary_relaxation_factor; 219 1551 : fxn_m1_tagid = _secondary_fxn_m1_tagid; 220 1551 : xn_m1_tagid = _secondary_xn_m1_tagid; 221 1551 : fxn_m2_tagid = _secondary_fxn_m2_tagid; 222 1551 : xn_m2_tagid = _secondary_xn_m2_tagid; 223 : } 224 : 225 5239 : NumericVector<Number> & solution = _solver_sys.solution(); 226 5239 : NumericVector<Number> & xn_m1 = _solver_sys.getVector(xn_m1_tagid); 227 5239 : NumericVector<Number> & fxn_m2 = _solver_sys.getVector(fxn_m2_tagid); 228 5239 : NumericVector<Number> & xn_m2 = _solver_sys.getVector(xn_m2_tagid); 229 : 230 : // Save the most recent evaluation of the coupled problem 231 5239 : NumericVector<Number> & fxn_m1 = _solver_sys.getVector(fxn_m1_tagid); 232 5239 : fxn_m1 = solution; 233 : 234 464555 : for (const auto & dof : target_dofs) 235 : { 236 : // Avoid 0 denominator issue 237 459316 : Real new_value = fxn_m1(dof); 238 459316 : if (!MooseUtils::absoluteFuzzyEqual(solution(dof) - xn_m1(dof) - fxn_m2(dof) + xn_m2(dof), 0)) 239 376662 : new_value = xn_m1(dof) - (solution(dof) - xn_m1(dof)) * (xn_m1(dof) - xn_m2(dof)) / 240 376662 : (solution(dof) - xn_m1(dof) - fxn_m2(dof) + xn_m2(dof)); 241 : 242 : // Relax update 243 459316 : new_value = relaxation_factor * new_value + (1 - relaxation_factor) * xn_m1(dof); 244 : 245 459316 : solution.set(dof, new_value); 246 : } 247 5239 : solution.close(); 248 5239 : _solver_sys.update(); 249 5239 : } 250 : 251 : void 252 14911 : SecantSolve::printFixedPointConvergenceHistory(Real initial_norm, 253 : const std::vector<Real> & timestep_begin_norms, 254 : const std::vector<Real> & timestep_end_norms) const 255 : { 256 14911 : _console << "\n 0 Secant initialization |R| = " 257 14911 : << Console::outputNorm(std::numeric_limits<Real>::max(), initial_norm) << '\n'; 258 : 259 14911 : Real max_norm_old = initial_norm; 260 87037 : for (unsigned int i = 0; i <= _fixed_point_it; ++i) 261 : { 262 72126 : Real max_norm = std::max(timestep_begin_norms[i], timestep_end_norms[i]); 263 72126 : std::stringstream secant_prefix; 264 72126 : if (i < 2) 265 27281 : secant_prefix << " Secant initialization |R| = "; 266 : else 267 44845 : secant_prefix << " Secant step |R| = "; 268 : 269 72126 : _console << std::setw(2) << i + 1 << secant_prefix.str() 270 72126 : << Console::outputNorm(max_norm_old, max_norm) << '\n'; 271 72126 : max_norm_old = max_norm; 272 72126 : } 273 14911 : }