# Contents of /trunk/paso/src/Transport_solve.cpp

Revision 6000 - (show annotations)
Tue Mar 1 00:24:43 2016 UTC (3 years, 1 month ago) by caltinay
File size: 9226 byte(s)
```a few more include rearrangements.

```
 1 2 /***************************************************************************** 3 * 4 * Copyright (c) 2003-2016 by The University of Queensland 5 6 * 7 * Primary Business: Queensland, Australia 8 * Licensed under the Open Software License version 3.0 9 10 * 11 * Development until 2012 by Earth Systems Science Computational Center (ESSCC) 12 * Development 2012-2013 by School of Earth Sciences 13 * Development from 2014 by Centre for Geoscience Computing (GeoComp) 14 * 15 *****************************************************************************/ 16 17 18 /****************************************************************************/ 19 20 /* Paso: Transport solver 21 * 22 * solves Mu_t=Ku+q 23 * 24 * using an operator splitting approach (K=L+D where L is row sum zero and D 25 * is a diagonal matrix): 26 * 27 * - Mu_t=Du+q u(0)=u (reactive part) 28 * - Mv_t=Lv v(0)=u(dt/2) (transport part uses flux correction (FCT)) 29 * - Mu_t=Du+q u(dt/2)=v(dt/2) (reactive part) 30 * 31 * to return u(dt) 32 * 33 */ 34 /****************************************************************************/ 35 36 /* Author: l.gross@uq.edu.au */ 37 38 /****************************************************************************/ 39 40 #include "Transport.h" 41 #include "FCT_Solver.h" 42 #include "ReactiveSolver.h" 43 #include "Solver.h" 44 #include "PasoUtil.h" 45 46 #include 47 48 namespace paso { 49 50 void TransportProblem::solve(double* u, double dt, double* u0, double* q, 51 Options* options) 52 { 53 const real_t EPSILON = escript::DataTypes::real_t_eps(); 54 const real_t LARGE_POSITIVE_FLOAT = escript::DataTypes::real_t_max(); 55 const double reduction_after_divergence_factor = 0.5; 56 const dim_t num_failures_max=50; 57 58 Performance pp; 59 ReactiveSolver* rsolver=NULL; 60 FCT_Solver* fctsolver=NULL; 61 62 dim_t i_substeps=0, n_substeps=1, num_failures=0; 63 double *u_save=NULL, *u2=NULL; 64 double dt2,t=0, dt3; 65 SolverResult errorCode=NoError; 66 const dim_t n = transport_matrix->getTotalNumRows(); 67 options->time_step_backtracking_used = false; 68 options->num_iter=0; 69 70 if (dt <= 0.) { 71 throw PasoException("TransportProblem::solve: dt must be positive."); 72 } else if (getBlockSize() > 1) { 73 throw PasoException("TransportProblem::solve: block size >1 " 74 "is not supported."); 75 } 76 if (options->verbose) { 77 if (options->ode_solver == PASO_BACKWARD_EULER) { 78 printf("TransportProblem::solve: Backward Euler is used (dt = %e)\n",dt); 79 } else if (options->ode_solver == PASO_LINEAR_CRANK_NICOLSON) { 80 printf("TransportProblem::solve: linear Crank-Nicolson is used (dt = %e).\n",dt); 81 } else if (options->ode_solver == PASO_CRANK_NICOLSON) { 82 printf("TransportProblem::solve: Crank-Nicolson is used (dt = %e).\n",dt); 83 } else { 84 throw PasoException("TransportProblem::solve: unknown ODE solver."); 85 } 86 } 87 getSafeTimeStepSize(); 88 // allocate memory 89 fctsolver = new FCT_Solver(shared_from_this(), options); 90 rsolver = new ReactiveSolver(shared_from_this()); 91 u_save = new double[n]; 92 u2 = new double[n]; 93 94 // let the show begin!!!! 95 const double dt_R = dt_max_R; 96 const double dt_T = dt_max_T; 97 dt2 = dt; 98 if (dt_R < LARGE_POSITIVE_FLOAT) 99 dt2 = std::min(dt_R*2, dt); // as we half the step size for the RT bit 100 if (dt_T < LARGE_POSITIVE_FLOAT) { 101 if (options->ode_solver == PASO_LINEAR_CRANK_NICOLSON || options->ode_solver == PASO_CRANK_NICOLSON) { 102 dt2 = std::min(dt_T, dt); 103 } // PASO_BACKWARD_EULER does not require a restriction 104 } 105 106 num_failures = 0; 107 util::copy(n, u, u0); // copy initial value to return 108 109 while((dt-t) > dt*sqrt(EPSILON)) { 110 n_substeps = ceil((dt-t)/dt2); 111 if (n_substeps <= 0) { 112 throw PasoException("TransportProblem::solve: time stepping break down."); 113 } else { 114 dt3 = (dt-t)/n_substeps; 115 if (options->verbose) { 116 std::cout << "TransportProblem::solve: number of substeps = " 117 << n_substeps << " with dt = " << dt3 << "." 118 << std::endl; 119 } 120 // initialize the iteration matrix 121 fctsolver->initialize(dt3, options, &pp); 122 rsolver->initialize(dt3/2, options); 123 errorCode = NoError; 124 125 // start iteration 126 for (i_substeps=0; i_substepsverbose) { 129 std::cout << "TransportProblem::solve: substep " 130 << i_substeps << " of " << n_substeps << " at t = " 131 << (t+dt3) << " (dt = " << dt3 << ")" << std::endl; 132 } 133 134 // create copy for restart in case of failure 135 util::copy(n, u_save, u); 136 // update u 137 138 // Mu_t=Du+q u(0)=u 139 errorCode = rsolver->solve(u2, u, q, options, &pp); 140 141 // Mv_t=Lv v(0)=u(dt/2) 142 if (errorCode == NoError) { 143 errorCode = fctsolver->update(u, u2, options, &pp); 144 145 } 146 // Mu_t=Du+q u(dt/2)=v(dt/2) 147 if (errorCode == NoError) { 148 errorCode = rsolver->solve(u2, u, q, options, &pp); 149 } 150 151 if (errorCode == NoError) { 152 num_failures = 0; 153 t += dt3; 154 util::copy(n, u, u2); 155 } 156 } 157 if (errorCode == MaxIterReached || errorCode == Divergence) { 158 // if num_failures_max failures in a row: give up 159 if (num_failures >= num_failures_max) { 160 throw PasoException("TransportProblem::solve: " 161 "No convergence after time step reductions."); 162 } else { 163 options->time_step_backtracking_used = true; 164 if (options->verbose) { 165 std::cout << "TransportProblem::solve: " 166 << "no convergence. Time step size is reduced." 167 << std::endl; 168 } 169 dt2 = dt3*reduction_after_divergence_factor; 170 num_failures++; 171 util::copy(n, u, u_save); // reset initial value 172 } 173 } else if (errorCode == InputError) { 174 throw PasoException("TransportProblem::solve: input error for solver."); 175 } else if (errorCode == MemoryError) { 176 throw PasoException("TransportProblem::solve: memory allocation failed."); 177 } else if (errorCode == Breakdown) { 178 throw PasoException("TransportProblem::solve: solver break down."); 179 } else if (errorCode == NegativeNormError) { 180 throw PasoException("TransportProblem::solve: negative norm."); 181 } else if (errorCode != NoError) { 182 throw PasoException("TransportProblem::solve: general error."); 183 } 184 } 185 } // end of time loop 186 187 delete fctsolver; 188 delete rsolver; 189 delete[] u_save; 190 delete[] u2; 191 } 192 193 double TransportProblem::getSafeTimeStepSize() const 194 { 195 double dt_max=0.; 196 const dim_t n = transport_matrix->getTotalNumRows(); 197 198 if (!valid_matrices) { 199 // set row-sum of mass_matrix 200 mass_matrix->rowSum(lumped_mass_matrix); 201 // check for positive entries in lumped_mass_matrix and set 202 // negative value for constraints 203 int fail = 0; 204 #pragma omp parallel 205 { 206 index_t fail_loc = 0; 207 #pragma omp for 208 for (index_t i=0; i 0) { 211 if (constraint_mask[i] > 0) 212 lumped_mass_matrix[i]=-1.; 213 } else { 214 fail_loc = 1; 215 } 216 } 217 #pragma omp critical 218 { 219 fail = std::max(fail, fail_loc); 220 } 221 } 222 #ifdef ESYS_MPI 223 int fail_loc = fail; 224 MPI_Allreduce(&fail_loc, &fail, 1, MPI_INT, MPI_MAX, mpi_info->comm); 225 #endif 226 if (fail > 0) 227 throw PasoException("TransportProblem::getSafeTimeStepSize: " 228 "negative mass matrix entries detected."); 229 // split off row-sum from transport_matrix 230 transport_matrix->makeZeroRowSums(reactive_matrix); 231 // get a copy of the main diagonal of the mass matrix 232 mass_matrix->copyFromMainDiagonal(main_diagonal_mass_matrix); 233 234 const double dt_R = ReactiveSolver::getSafeTimeStepSize(shared_from_this()); 235 const double dt_T = FCT_Solver::getSafeTimeStepSize(shared_from_this()); 236 dt_max_R = dt_R; 237 dt_max_T = dt_T; 238 valid_matrices = true; 239 dt_max = std::min(2*dt_R, dt_T); 240 } else { 241 // factor 2 as we use operator splitting 242 dt_max = std::min(2*dt_max_R, dt_max_T); 243 } 244 return dt_max; 245 } 246 247 } // namespace paso 248

## Properties

Name Value
svn:mergeinfo /branches/4.0fordebian/paso/src/Transport_solve.cpp:5567-5588 /branches/amg_from_3530/paso/src/Transport_solve.cpp:3531-3826 /branches/lapack2681/paso/src/Transport_solve.cpp:2682-2741 /branches/pasowrap/paso/src/Transport_solve.cpp:3661-3674 /branches/py3_attempt2/paso/src/Transport_solve.cpp:3871-3891 /branches/restext/paso/src/Transport_solve.cpp:2610-2624 /branches/ripleygmg_from_3668/paso/src/Transport_solve.cpp:3669-3791 /branches/stage3.0/paso/src/Transport_solve.cpp:2569-2590 /branches/symbolic_from_3470/paso/src/Transport_solve.cpp:3471-3974 /branches/symbolic_from_3470/ripley/test/python/paso/src/Transport_solve.cpp:3517-3974 /release/3.0/paso/src/Transport_solve.cpp:2591-2601 /release/4.0/paso/src/Transport_solve.cpp:5380-5406 /trunk/paso/src/Transport_solve.cpp:4257-4344 /trunk/ripley/test/python/paso/src/Transport_solve.cpp:3480-3515