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

Revision 3825 - (show annotations)
Tue Feb 14 06:11:58 2012 UTC (7 years, 2 months ago) by gross
File MIME type: text/plain
File size: 9151 byte(s)
```don't call MPI withing parallel regions
```
 1 2 /******************************************************* 3 * 4 * Copyright (c) 2003-2010 by University of Queensland 5 * Earth Systems Science Computational Center (ESSCC) 6 * http://www.uq.edu.au/esscc 7 * 8 * Primary Business: Queensland, Australia 9 * Licensed under the Open Software License version 3.0 10 * http://www.opensource.org/licenses/osl-3.0.php 11 * 12 *******************************************************/ 13 14 15 /**************************************************************/ 16 17 /* Paso: Transport solver 18 * 19 * solves Mu_t=Ku+q 20 * 21 * using an operator splitting approach (K=L+D where L is row sum zero and D is a diagonal matrix): 22 * 23 * - Mu_t=Du+q u(0)=u (reactive part) 24 * - Mv_t=Lv v(0)=u(dt/2) (transport part uses flux correction (FCT)) 25 * - Mu_t=Du+q u(dt/2)=v(dt/2) (reactive part) 26 * 27 * to return u(dt) 28 * 29 */ 30 /**************************************************************/ 31 32 /* Author: l.gross@uq.edu.au */ 33 34 /**************************************************************/ 35 36 37 38 #include "Transport.h" 39 #include "FCT_Solver.h" 40 #include "ReactiveSolver.h" 41 #include "Solver.h" 42 #include "PasoUtil.h" 43 44 45 void Paso_TransportProblem_solve(Paso_TransportProblem* fctp, double* u, double dt, double* u0, double* q, Paso_Options* options) { 46 const double reduction_after_divergence_factor = 0.5; 47 const dim_t num_failures_max=50; 48 49 Paso_Performance pp; 50 Paso_ReactiveSolver *rsolver=NULL; 51 Paso_FCT_Solver *fctsolver=NULL; 52 53 54 dim_t i_substeps=0, n_substeps=1, num_failures=0; 55 double *u_save=NULL, *u2=NULL; 56 double dt2,t=0, dt3; 57 err_t errorCode=SOLVER_NO_ERROR; 58 const dim_t n=Paso_SystemMatrix_getTotalNumRows(fctp->transport_matrix); 59 const dim_t blockSize=Paso_TransportProblem_getBlockSize(fctp); 60 options->time_step_backtracking_used=FALSE; 61 options->num_iter=0; 62 63 if (dt<=0.) { 64 Esys_setError(VALUE_ERROR,"Paso_TransportProblem_solve: dt must be positive."); 65 } 66 if (blockSize>1) { 67 Esys_setError(VALUE_ERROR,"Paso_TransportProblem_solve: block size >1 is not supported."); 68 } 69 if (options->verbose) { 70 if (options->ode_solver == PASO_BACKWARD_EULER) { 71 printf("Paso_TransportProblem_solve: Backward Euler is used (dt = %e)\n",dt); 72 } else if (options->ode_solver == PASO_LINEAR_CRANK_NICOLSON) { 73 printf("Paso_TransportProblem_solve: linear Crank-Nicolson is used (dt = %e).\n",dt); 74 } else if (options->ode_solver == PASO_CRANK_NICOLSON) { 75 printf("Paso_TransportProblem_solve: Crank-Nicolson is used (dt = %e).\n",dt); 76 } else { 77 Esys_setError(VALUE_ERROR,"Paso_TransportProblem_solve: unknown ODE solver."); 78 } 79 } 80 if (Esys_noError()) { 81 Paso_TransportProblem_getSafeTimeStepSize(fctp); 82 /* 83 * allocate memory 84 * 85 */ 86 fctsolver=Paso_FCT_Solver_alloc(fctp,options); 87 rsolver=Paso_ReactiveSolver_alloc(fctp); 88 u_save=TMPMEMALLOC(n,double); 89 u2=TMPMEMALLOC(n,double); 90 Esys_checkPtr(u_save); 91 Esys_checkPtr(u2); 92 } 93 if (Esys_noError()) { 94 /* 95 * let the show begin!!!! 96 * 97 */ 98 const double dt_R=fctp->dt_max_R; 99 const double dt_T=fctp->dt_max_T; 100 dt2=dt; 101 if (dt_R < LARGE_POSITIVE_FLOAT) dt2 = MIN(dt_R *2, dt); /* as we half the steps size for the RT bit */ 102 if (dt_T < LARGE_POSITIVE_FLOAT) { 103 if ( ( options->ode_solver == PASO_LINEAR_CRANK_NICOLSON) || ( options->ode_solver == PASO_CRANK_NICOLSON ) ) { 104 dt2 = MIN(dt_T, dt); 105 } /* PASO_BACKWARD_EULER does not require a restriction */ 106 } 107 108 num_failures=0; 109 Paso_Copy(n,u,u0); /* copy initial value to return */ 110 111 while( (dt-t)>dt*sqrt(EPSILON) && Esys_noError()) { 112 113 n_substeps=ceil((dt-t)/dt2); 114 if ( n_substeps <= 0) { 115 Esys_setError(VALUE_ERROR,"Paso_TransportProblem_solve: time stepping break down."); 116 } else { 117 dt3=(dt-t)/n_substeps; 118 if (options->verbose) printf("Paso_TransportProblem_solve: number of substeps = %d with dt = %e.\n",n_substeps,dt3); 119 /* initialize the iteration matrix */ 120 Paso_FCT_Solver_initialize(dt3, fctsolver, options, &pp); 121 Paso_ReactiveSolver_initialize(dt3/2, rsolver, options); 122 errorCode = SOLVER_NO_ERROR; 123 124 /* start iteration */ 125 for (i_substeps =0; (i_substepsverbose) printf("Paso_TransportProblem_solve: substep %d of %d at t = %e (dt = %e)\n",i_substeps,n_substeps,t+dt3,dt3); 127 Paso_Copy(n,u_save, u); /* create copy for restart in case of failure */ 128 /* updates u */ 129 errorCode=Paso_ReactiveSolver_solve(rsolver,fctp,u2, u ,q, options, &pp); /* Mu_t=Du+q u(0)=u */ 130 if (errorCode == SOLVER_NO_ERROR) { 131 errorCode=Paso_FCT_Solver_update(fctsolver, u, u2,options, &pp); /* Mv_t=Lv v(0)=u(dt/2) */ 132 133 } 134 if (errorCode == SOLVER_NO_ERROR) { 135 errorCode=Paso_ReactiveSolver_solve(rsolver,fctp,u2, u, q, options, &pp); /* Mu_t=Du+q u(dt/2)=v(dt/2) */ 136 } 137 138 139 /* 140 * error handling: 141 */ 142 if (errorCode == SOLVER_NO_ERROR) { 143 num_failures=0; 144 t+=dt3; 145 Paso_Copy(n,u, u2); 146 } 147 } 148 if (errorCode == SOLVER_NO_ERROR) { 149 } else if ( (errorCode == SOLVER_MAXITER_REACHED) || (errorCode == SOLVER_DIVERGENCE) ) { 150 /* if num_failures_max failures in a row: give up */ 151 if (num_failures >= num_failures_max) { 152 Esys_setError(VALUE_ERROR,"Paso_TransportProblem_solve: No convergence after time step reductions."); 153 } else { 154 options->time_step_backtracking_used=TRUE; 155 if (options->verbose) printf("Paso_TransportProblem_solve: no convergence. Time step size is reduced.\n"); 156 dt2=dt3*reduction_after_divergence_factor; 157 num_failures++; 158 Paso_Copy(n,u, u_save); /* reset initial value */ 159 } 160 } else if (errorCode == SOLVER_INPUT_ERROR) { 161 Esys_setError(VALUE_ERROR,"Paso_TransportProblem_solve: input error for solver."); 162 } else if (errorCode == SOLVER_MEMORY_ERROR) { 163 Esys_setError(MEMORY_ERROR,"Paso_TransportProblem_solve: memory allocation failed."); 164 } else if (errorCode == SOLVER_BREAKDOWN) { 165 Esys_setError(VALUE_ERROR,"Paso_TransportProblem_solve: solver break down."); 166 } else if (errorCode == SOLVER_NEGATIVE_NORM_ERROR) { 167 Esys_setError(VALUE_ERROR,"Paso_TransportProblem_solve: negative norm."); 168 } else { 169 Esys_setError(SYSTEM_ERROR,"Paso_TransportProblem_solve: general error."); 170 } 171 } 172 } 173 } /* end of time loop */ 174 /* 175 * clean-up: 176 * 177 */ 178 179 Paso_FCT_Solver_free(fctsolver); 180 Paso_ReactiveSolver_free(rsolver); 181 TMPMEMFREE(u_save); 182 TMPMEMFREE(u2); 183 } 184 185 186 double Paso_TransportProblem_getSafeTimeStepSize(Paso_TransportProblem* fctp) 187 { 188 double dt_max=0., dt_R=LARGE_POSITIVE_FLOAT, dt_T=LARGE_POSITIVE_FLOAT; 189 index_t fail_loc=0; 190 const dim_t n=Paso_SystemMatrix_getTotalNumRows(fctp->transport_matrix); 191 192 if ( ! fctp->valid_matrices) { 193 /* set row-sum of mass_matrix */ 194 Paso_SystemMatrix_rowSum(fctp->mass_matrix,fctp->lumped_mass_matrix); 195 { 196 /* check for positive entries in lumped_mass_matrix and set negative value for constraints */ 197 index_t fail=0, i; 198 #pragma omp parallel private(i,fail_loc ) 199 { 200 201 fail_loc=0; 202 #pragma omp for schedule(static) 203 for (i=0;ilumped_mass_matrix[i]; 205 if ( m_i > 0 ) { 206 if (fctp->constraint_mask[i] > 0 ) fctp->lumped_mass_matrix[i]=-1.; 207 } else { 208 fail_loc=1; 209 } 210 } 211 #pragma omp critical 212 { 213 fail=MAX(fail, fail_loc); 214 } 215 } 216 #ifdef ESYS_MPI 217 { 218 fail_loc=fail; 219 MPI_Allreduce(&fail_loc, &fail, 1, MPI_INT, MPI_MAX, fctp->mpi_info->comm); 220 } 221 #endif 222 if (fail < 0 ) 223 Esys_setError(VALUE_ERROR, "Paso_FCTSolver_getSafeTimeStepSize: negative mass matrix entries detected."); 224 } 225 /* split off row-sum from transport_matrix */ 226 Paso_SystemMatrix_makeZeroRowSums(fctp->transport_matrix,fctp->reactive_matrix); 227 /* get a copy of the main diagonal of the mass matrix */ 228 Paso_SystemMatrix_copyFromMainDiagonal(fctp->mass_matrix,fctp->main_diagonal_mass_matrix); 229 230 if (Esys_noError()) { 231 dt_R=Paso_ReactiveSolver_getSafeTimeStepSize(fctp); 232 dt_T=Paso_FCT_Solver_getSafeTimeStepSize(fctp); 233 } 234 if (Esys_noError()) { 235 fctp->dt_max_R=dt_R; 236 fctp->dt_max_T=dt_T; 237 fctp->valid_matrices=TRUE; 238 dt_max=MIN( 2 * dt_R,dt_T); 239 } 240 } else { 241 dt_max=MIN(2 * fctp->dt_max_R,fctp->dt_max_T); /* factor 2 as we use operator splitting */ 242 } 243 return dt_max; 244 }

 ViewVC Help Powered by ViewVC 1.1.26