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

Revision 3005 - (show annotations)
Thu Apr 22 05:59:31 2010 UTC (9 years, 6 months ago) by gross
File MIME type: text/plain
File size: 6118 byte(s)
```early call of setPreconditioner in FCT solver removed.
```
 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 Paso_Performance pp; 45 Paso_Function *fctfunction=NULL; 46 Paso_ReactiveSolver *rsolver=NULL; 47 const dim_t FAILURES_MAX=100; 48 dim_t i_substeps, Failed=0; 49 double *u_save=NULL; 50 double dt_max=LARGE_POSITIVE_FLOAT, dt2,t=0; 51 err_t errorCode=SOLVER_NO_ERROR; 52 const dim_t n=Paso_SystemMatrix_getTotalNumRows(fctp->transport_matrix); 53 const dim_t blockSize=Paso_TransportProblem_getBlockSize(fctp); 54 55 56 if (dt<=0.) { 57 Paso_setError(VALUE_ERROR,"Paso_TransportProblem_solve: dt must be positive."); 58 } 59 if (blockSize>1) { 60 Paso_setError(VALUE_ERROR,"Paso_TransportProblem_solve: block size >0 is not supported."); 61 } 62 if (Paso_noError()) { 63 dt_max=Paso_TransportProblem_getSafeTimeStepSize(fctp); 64 /* 65 * allocate memory 66 * 67 */ 68 fctfunction=Paso_FCTSolver_Function_alloc(fctp,options); 69 rsolver=Paso_ReactiveSolver_alloc(fctp); 70 u_save=TMPMEMALLOC(n,double); 71 Paso_checkPtr(u_save); 72 } 73 if (Paso_noError()) { 74 /* 75 * let the show begin!!!! 76 * 77 */ 78 i_substeps=0; 79 80 /* while(i_substepsverbose) printf("substep step %d at t=%e (step size= %e)\n",i_substeps+1,t+dt2,dt2); 93 94 /* updates u */ 95 errorCode=Paso_ReactiveSolver_solve(rsolver,fctp,u,dt2/2,q, options, &pp); /* Mu_t=Du+q u(0)=u */ 96 if (errorCode == NO_ERROR) { 97 errorCode=Paso_FCTSolver_solve(fctfunction, u,dt2,options, &pp); /* Mv_t=Lv v(0)=u(dt/2) */ 98 } 99 if (errorCode == NO_ERROR) { 100 errorCode=Paso_ReactiveSolver_solve(rsolver,fctp,u,dt2/2,q, options, &pp); /* Mu_t=Du+q u(dt/2)=v(dt/2) */ 101 } 102 /* 103 * error handeling: 104 */ 105 if (errorCode == SOLVER_NO_ERROR) { 106 Failed=0; 107 i_substeps++; 108 t+=dt2; 109 if (fctp->dt_max < LARGE_POSITIVE_FLOAT) { 110 dt2=MIN3(dt_max,dt2*1.1,dt-t); 111 } else { 112 dt2=MIN(dt2*1.1,dt-t); 113 } 114 Paso_Copy(n,u_save, u); 115 } else if ( (errorCode == SOLVER_MAXITER_REACHED) || (errorCode == SOLVER_DIVERGENCE) ) { 116 /* if FAILURES_MAX failures in a row: give up */ 117 if (Failed > FAILURES_MAX) { 118 Paso_setError(VALUE_ERROR,"Paso_TransportProblem_solve: no convergence after time step reduction."); 119 } else { 120 if (options->verbose) printf("Paso_TransportProblem_solve: no convergence. Time step size is reduced.\n"); 121 dt2*=0.5; 122 Failed++; 123 Paso_Copy(n,u, u_save); 124 } 125 } else if (errorCode == SOLVER_INPUT_ERROR) { 126 Paso_setError(VALUE_ERROR,"Paso_TransportProblem_solve: input error for solver."); 127 } else if (errorCode == SOLVER_MEMORY_ERROR) { 128 Paso_setError(MEMORY_ERROR,"Paso_TransportProblem_solve: memory allication failed."); 129 } else if (errorCode == SOLVER_BREAKDOWN) { 130 Paso_setError(VALUE_ERROR,"Paso_TransportProblem_solve: solver break down."); 131 } else if (errorCode == SOLVER_NEGATIVE_NORM_ERROR) { 132 Paso_setError(VALUE_ERROR,"Paso_TransportProblem_solve: negative norm."); 133 } else { 134 Paso_setError(SYSTEM_ERROR,"Paso_TransportProblem_solve: general error."); 135 } 136 } 137 } 138 /* 139 * clean-up: 140 * 141 */ 142 Paso_FCTSolver_Function_free(fctfunction); 143 Paso_ReactiveSolver_free(rsolver); 144 TMPMEMFREE(u_save); 145 146 } 147 double Paso_TransportProblem_getSafeTimeStepSize(Paso_TransportProblem* fctp) 148 { 149 double dt_max, dt1, dt2; 150 if ( ! fctp->valid_matrices) { 151 /* set row-sum of mass_matrix */ 152 Paso_SystemMatrix_rowSum(fctp->mass_matrix,fctp->lumped_mass_matrix); 153 /* split off row-sum from transport_matrix */ 154 Paso_SystemMatrix_makeZeroRowSums(fctp->transport_matrix,fctp->reactive_matrix); 155 /* get a copy of the main diagonal of the mass matrix */ 156 Paso_SystemMatrix_copyFromMainDiagonal(fctp->mass_matrix,fctp->main_diagonal_mass_matrix); 157 158 159 /* XXXXXXXXX hier geht es weiter XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX */ 160 dt1=Paso_ReactiveSolver_getSafeTimeStepSize(fctp); 161 dt2=Paso_FCTSolver_getSafeTimeStepSize(fctp); 162 dt_max=MIN(dt1,dt2); 163 164 if (dt_max <= 0.) { 165 Paso_setError(TYPE_ERROR,"Paso_TransportProblem_solve: dt must be positive."); 166 } 167 fctp->dt_max=dt_max; 168 fctp->valid_matrices=Paso_noError(); 169 } 170 return fctp->dt_max; 171 }

 ViewVC Help Powered by ViewVC 1.1.26