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

Revision 4873 - (show annotations)
Wed Apr 16 06:38:51 2014 UTC (5 years, 5 months ago) by caltinay
File size: 9764 byte(s)
whitespace only changes.

 1 2 /***************************************************************************** 3 * 4 * Copyright (c) 2003-2014 by 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 namespace paso { 47 48 void TransportProblem::solve(double* u, double dt, double* u0, double* q, 49 Options* options) 50 { 51 const double reduction_after_divergence_factor = 0.5; 52 const dim_t num_failures_max=50; 53 54 Performance pp; 55 ReactiveSolver* rsolver=NULL; 56 FCT_Solver* fctsolver=NULL; 57 58 dim_t i_substeps=0, n_substeps=1, num_failures=0; 59 double *u_save=NULL, *u2=NULL; 60 double dt2,t=0, dt3; 61 err_t errorCode=SOLVER_NO_ERROR; 62 const dim_t n = transport_matrix->getTotalNumRows(); 63 options->time_step_backtracking_used = false; 64 options->num_iter=0; 65 66 if (dt <= 0.) { 67 Esys_setError(VALUE_ERROR, "TransportProblem::solve: dt must be positive."); 68 } else if (getBlockSize() > 1) { 69 Esys_setError(VALUE_ERROR, "TransportProblem::solve: block size >1 " 70 "is not supported."); 71 } 72 if (options->verbose) { 73 if (options->ode_solver == PASO_BACKWARD_EULER) { 74 printf("TransportProblem::solve: Backward Euler is used (dt = %e)\n",dt); 75 } else if (options->ode_solver == PASO_LINEAR_CRANK_NICOLSON) { 76 printf("TransportProblem::solve: linear Crank-Nicolson is used (dt = %e).\n",dt); 77 } else if (options->ode_solver == PASO_CRANK_NICOLSON) { 78 printf("TransportProblem::solve: Crank-Nicolson is used (dt = %e).\n",dt); 79 } else { 80 Esys_setError(VALUE_ERROR, "TransportProblem::solve: unknown ODE solver."); 81 } 82 } 83 if (Esys_noError()) { 84 getSafeTimeStepSize(); 85 // allocate memory 86 fctsolver = new FCT_Solver(shared_from_this(), options); 87 rsolver = new ReactiveSolver(shared_from_this()); 88 u_save = new double[n]; 89 u2 = new double[n]; 90 } 91 if (Esys_noError()) { 92 // let the show begin!!!! 93 const double dt_R = dt_max_R; 94 const double dt_T = dt_max_T; 95 dt2 = dt; 96 if (dt_R < LARGE_POSITIVE_FLOAT) 97 dt2 = MIN(dt_R*2, dt); // as we half the step size for the RT bit 98 if (dt_T < LARGE_POSITIVE_FLOAT) { 99 if (options->ode_solver == PASO_LINEAR_CRANK_NICOLSON || options->ode_solver == PASO_CRANK_NICOLSON) { 100 dt2 = MIN(dt_T, dt); 101 } // PASO_BACKWARD_EULER does not require a restriction 102 } 103 104 num_failures = 0; 105 util::copy(n, u, u0); // copy initial value to return 106 107 while((dt-t) > dt*sqrt(EPSILON) && Esys_noError()) { 108 n_substeps = ceil((dt-t)/dt2); 109 if (n_substeps <= 0) { 110 Esys_setError(VALUE_ERROR, "TransportProblem::solve: time stepping break down."); 111 } else { 112 dt3 = (dt-t)/n_substeps; 113 if (options->verbose) 114 printf("TransportProblem::solve: number of substeps = %d " 115 "with dt = %e.\n", n_substeps, dt3); 116 // initialize the iteration matrix 117 fctsolver->initialize(dt3, options, &pp); 118 rsolver->initialize(dt3/2, options); 119 errorCode = SOLVER_NO_ERROR; 120 121 // start iteration 122 for (i_substeps=0; i_substepsverbose) 126 printf("TransportProblem::solve: substep %d of %d at t = %e (dt = %e)\n",i_substeps,n_substeps,t+dt3,dt3); 127 128 // create copy for restart in case of failure 129 util::copy(n, u_save, u); 130 // update u 131 132 // Mu_t=Du+q u(0)=u 133 errorCode = rsolver->solve(u2, u, q, options, &pp); 134 135 // Mv_t=Lv v(0)=u(dt/2) 136 if (errorCode == SOLVER_NO_ERROR) { 137 errorCode = fctsolver->update(u, u2, options, &pp); 138 139 } 140 // Mu_t=Du+q u(dt/2)=v(dt/2) 141 if (errorCode == SOLVER_NO_ERROR) { 142 errorCode = rsolver->solve(u2, u, q, options, &pp); 143 } 144 145 if (errorCode == SOLVER_NO_ERROR) { 146 num_failures = 0; 147 t += dt3; 148 util::copy(n, u, u2); 149 } 150 } 151 if (errorCode == SOLVER_MAXITER_REACHED || errorCode == SOLVER_DIVERGENCE) { 152 // if num_failures_max failures in a row: give up 153 if (num_failures >= num_failures_max) { 154 Esys_setError(VALUE_ERROR, "TransportProblem::solve: " 155 "No convergence after time step reductions."); 156 } else { 157 options->time_step_backtracking_used = true; 158 if (options->verbose) 159 printf("TransportProblem::solve: no convergence. " 160 "Time step size is reduced.\n"); 161 dt2 = dt3*reduction_after_divergence_factor; 162 num_failures++; 163 util::copy(n, u, u_save); // reset initial value 164 } 165 } else if (errorCode == SOLVER_INPUT_ERROR) { 166 Esys_setError(VALUE_ERROR, "TransportProblem::solve: " 167 "input error for solver."); 168 } else if (errorCode == SOLVER_MEMORY_ERROR) { 169 Esys_setError(MEMORY_ERROR, "TransportProblem::solve: " 170 "memory allocation failed."); 171 } else if (errorCode == SOLVER_BREAKDOWN) { 172 Esys_setError(VALUE_ERROR, "TransportProblem::solve: " 173 "solver break down."); 174 } else if (errorCode == SOLVER_NEGATIVE_NORM_ERROR) { 175 Esys_setError(VALUE_ERROR, "TransportProblem::solve: " 176 "negative norm."); 177 } else if (errorCode != SOLVER_NO_ERROR) { 178 Esys_setError(SYSTEM_ERROR, "TransportProblem::solve: " 179 "general error."); 180 } 181 } 182 } 183 } // end of time loop 184 185 delete fctsolver; 186 delete rsolver; 187 delete[] u_save; 188 delete[] u2; 189 } 190 191 double TransportProblem::getSafeTimeStepSize() 192 { 193 double dt_max=0.; 194 const dim_t n = transport_matrix->getTotalNumRows(); 195 196 if (!valid_matrices) { 197 // set row-sum of mass_matrix 198 mass_matrix->rowSum(lumped_mass_matrix); 199 // check for positive entries in lumped_mass_matrix and set 200 // negative value for constraints 201 int fail = 0; 202 #pragma omp parallel 203 { 204 index_t fail_loc = 0; 205 #pragma omp for 206 for (index_t i=0; i 0) { 209 if (constraint_mask[i] > 0) 210 lumped_mass_matrix[i]=-1.; 211 } else { 212 fail_loc = 1; 213 } 214 } 215 #pragma omp critical 216 { 217 fail = MAX(fail, fail_loc); 218 } 219 } 220 #ifdef ESYS_MPI 221 int fail_loc = fail; 222 MPI_Allreduce(&fail_loc, &fail, 1, MPI_INT, MPI_MAX, mpi_info->comm); 223 #endif 224 if (fail > 0) 225 Esys_setError(VALUE_ERROR, "TransportProblem::getSafeTimeStepSize: " 226 "negative mass matrix entries detected."); 227 // split off row-sum from transport_matrix 228 transport_matrix->makeZeroRowSums(reactive_matrix); 229 // get a copy of the main diagonal of the mass matrix 230 mass_matrix->copyFromMainDiagonal(main_diagonal_mass_matrix); 231 232 if (Esys_noError()) { 233 const double dt_R = ReactiveSolver::getSafeTimeStepSize(shared_from_this()); 234 const double dt_T = FCT_Solver::getSafeTimeStepSize(shared_from_this()); 235 dt_max_R = dt_R; 236 dt_max_T = dt_T; 237 valid_matrices = true; 238 dt_max = MIN(2*dt_R, dt_T); 239 } 240 } else { 241 // factor 2 as we use operator splitting 242 dt_max = 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/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 /trunk/paso/src/Transport_solve.cpp:4257-4344 /trunk/ripley/test/python/paso/src/Transport_solve.cpp:3480-3515