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

Revision 3014 - (show annotations)
Wed Apr 28 04:05:21 2010 UTC (9 years, 9 months ago) by gross
File MIME type: text/plain
File size: 7746 byte(s)
```a more stable handling of divergence in the CN or BE
```
 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 #include "Transport.h" 37 #include "FCTSolver.h" 38 #include "ReactiveSolver.h" 39 #include "Solver.h" 40 #include "PasoUtil.h" 41 42 43 void Paso_TransportProblem_solve(Paso_TransportProblem* fctp, double* u, double dt, double* u0, double* q, Paso_Options* options) { 44 const double reduction_after_divergence_factor = 0.5; 45 const dim_t num_failures_max=50; 46 const double increase_after_convergence_factor = 1.1; /* has not much of an impact if substepping is beeing used */ 47 48 Paso_Performance pp; 49 Paso_Function *fctfunction=NULL; 50 Paso_ReactiveSolver *rsolver=NULL; 51 52 dim_t i_substeps, n_substeps, num_failures=0; 53 double *u_save=NULL; 54 double dt_max=LARGE_POSITIVE_FLOAT, dt2,t=0, dt3; 55 err_t errorCode=SOLVER_NO_ERROR; 56 const dim_t n=Paso_SystemMatrix_getTotalNumRows(fctp->transport_matrix); 57 const dim_t blockSize=Paso_TransportProblem_getBlockSize(fctp); 58 59 60 if (dt<=0.) { 61 Paso_setError(VALUE_ERROR,"Paso_TransportProblem_solve: dt must be positive."); 62 } 63 if (blockSize>1) { 64 Paso_setError(VALUE_ERROR,"Paso_TransportProblem_solve: block size >0 is not supported."); 65 } 66 if (options->verbose) { 67 if (fctp->useBackwardEuler) { 68 printf("Paso_TransportProblem_solve: Backward Euler is used (dt = %e)\n",dt); 69 } else { 70 printf("Paso_TransportProblem_solve: Crank-Nicolson is used (dt = %e).\n",dt); 71 } 72 } 73 if (Paso_noError()) { 74 dt_max=Paso_TransportProblem_getSafeTimeStepSize(fctp); 75 /* 76 * allocate memory 77 * 78 */ 79 fctfunction=Paso_FCTSolver_Function_alloc(fctp,options); 80 rsolver=Paso_ReactiveSolver_alloc(fctp); 81 u_save=TMPMEMALLOC(n,double); 82 Paso_checkPtr(u_save); 83 } 84 if (Paso_noError()) { 85 /* 86 * let the show begin!!!! 87 * 88 */ 89 90 /* while(i_substepsdt*sqrt(EPSILON) && Paso_noError()) { 103 104 n_substeps=ceil((dt-t)/dt2); 105 dt3=(dt-t)/n_substeps; 106 if (options->verbose) printf("Paso_TransportProblem_solve: number of substeps = %d with dt = %e.\n",n_substeps,dt3); 107 i_substeps=0; 108 /* initialize the iteration matrix */ 109 Paso_FCTSolver_Function_initialize(dt3, fctp,options, &pp); 110 Paso_ReactiveSolver_initialize(dt3, fctp, options); 111 /* start iteration */ 112 for (i_substeps =0; (i_substepsverbose) printf("Paso_TransportProblem_solve: substep step %d of %d at t = %e (dt = %e)\n",i_substeps,n_substeps,t+dt3,dt3); 114 /* updates u */ 115 errorCode=Paso_ReactiveSolver_solve(rsolver,fctp,u,dt3/2,q, options, &pp); /* Mu_t=Du+q u(0)=u */ 116 if (errorCode == SOLVER_NO_ERROR) { 117 errorCode=Paso_FCTSolver_solve(fctfunction, u,dt3,options, &pp); /* Mv_t=Lv v(0)=u(dt/2) */ 118 } 119 if (errorCode == SOLVER_NO_ERROR) { 120 errorCode=Paso_ReactiveSolver_solve(rsolver,fctp,u,dt3/2,q, options, &pp); /* Mu_t=Du+q u(dt/2)=v(dt/2) */ 121 } 122 if (errorCode == SOLVER_NO_ERROR) { 123 t+=dt3; 124 Paso_Copy(n,u_save, u); /* copy u to u_save */ 125 } 126 } 127 /* 128 * error handeling: 129 */ 130 if (errorCode == SOLVER_NO_ERROR) { 131 num_failures=0; 132 i_substeps++; 133 if (fctp->dt_max < LARGE_POSITIVE_FLOAT) { 134 fctp->dt_max=MIN3(dt_max, dt3*increase_after_convergence_factor ,fctp->dt_failed/increase_after_convergence_factor); 135 } else { 136 fctp->dt_max=MIN(dt3*increase_after_convergence_factor ,fctp->dt_failed/increase_after_convergence_factor); 137 } 138 } else if ( (errorCode == SOLVER_MAXITER_REACHED) || (errorCode == SOLVER_DIVERGENCE) ) { 139 /* if num_failures_max failures in a row: give up */ 140 fctp->dt_failed=MIN(fctp->dt_failed, dt3); 141 if (num_failures >= num_failures_max) { 142 Paso_setError(VALUE_ERROR,"Paso_TransportProblem_solve: No convergence after time step reductions."); 143 } else { 144 if (options->verbose) printf("Paso_TransportProblem_solve: no convergence. Time step size is reduced.\n"); 145 dt2=dt3*reduction_after_divergence_factor; 146 num_failures++; 147 Paso_Copy(n,u, u_save); /* reset initial value */ 148 errorCode = SOLVER_NO_ERROR; 149 } 150 151 } else if (errorCode == SOLVER_INPUT_ERROR) { 152 Paso_setError(VALUE_ERROR,"Paso_TransportProblem_solve: input error for solver."); 153 } else if (errorCode == SOLVER_MEMORY_ERROR) { 154 Paso_setError(MEMORY_ERROR,"Paso_TransportProblem_solve: memory allication failed."); 155 } else if (errorCode == SOLVER_BREAKDOWN) { 156 Paso_setError(VALUE_ERROR,"Paso_TransportProblem_solve: solver break down."); 157 } else if (errorCode == SOLVER_NEGATIVE_NORM_ERROR) { 158 Paso_setError(VALUE_ERROR,"Paso_TransportProblem_solve: negative norm."); 159 } else { 160 Paso_setError(SYSTEM_ERROR,"Paso_TransportProblem_solve: general error."); 161 } 162 } 163 } 164 /* 165 * clean-up: 166 * 167 */ 168 Paso_FCTSolver_Function_free(fctfunction); 169 Paso_ReactiveSolver_free(rsolver); 170 TMPMEMFREE(u_save); 171 172 } 173 174 175 double Paso_TransportProblem_getSafeTimeStepSize(Paso_TransportProblem* fctp) 176 { 177 double dt_max, dt1, dt2; 178 if ( ! fctp->valid_matrices) { 179 fctp->dt_failed=LARGE_POSITIVE_FLOAT; 180 /* set row-sum of mass_matrix */ 181 Paso_SystemMatrix_rowSum(fctp->mass_matrix,fctp->lumped_mass_matrix); 182 /* split off row-sum from transport_matrix */ 183 Paso_SystemMatrix_makeZeroRowSums(fctp->transport_matrix,fctp->reactive_matrix); 184 /* get a copy of the main diagonal of the mass matrix */ 185 Paso_SystemMatrix_copyFromMainDiagonal(fctp->mass_matrix,fctp->main_diagonal_mass_matrix); 186 187 188 /* XXXXXXXXX hier geht es weiter XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX */ 189 dt1=Paso_ReactiveSolver_getSafeTimeStepSize(fctp); 190 dt2=Paso_FCTSolver_getSafeTimeStepSize(fctp); 191 dt_max=MIN(dt1,dt2); 192 193 if (dt_max <= 0.) { 194 Paso_setError(TYPE_ERROR,"Paso_TransportProblem_solve: dt must be positive."); 195 } 196 fctp->dt_max=dt_max; 197 fctp->valid_matrices=Paso_noError(); 198 } 199 return fctp->dt_max; 200 }