Contents of /trunk/paso/src/SolverFCT_solve.c

Revision 1375 - (show annotations)
Wed Jan 9 00:15:05 2008 UTC (12 years, 8 months ago) by gross
File MIME type: text/plain
File size: 6892 byte(s)
```bug in interpolation at reduced face elements fixed.
```
 1 /* \$Id:\$ */ 2 3 /******************************************************* 4 * 5 * Copyright 2007 by University of Queensland 6 * 7 * http://esscc.uq.edu.au 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 /* Paso: Flux correction transport solver 17 * 18 * solves Mu_t=Du+Ku+q 19 * 20 * where is D is diffusive (not checked) 21 * - D is symmetric 22 * - row sums are equal to zero. 23 * and K is the advective part. 24 * 25 * u(0) >= 0 26 * 27 * intially fctp->transport_matrix defines the diffusive part 28 * but the matrix is updated by the adevctive part + artificial diffusion 29 * 30 */ 31 /**************************************************************/ 32 33 /* Author: l.gross@uq.edu.au */ 34 35 /**************************************************************/ 36 37 #include "Paso.h" 38 #include "Solver.h" 39 #include "SolverFCT.h" 40 #include "escript/blocktimer.h" 41 #ifdef _OPENMP 42 #include 43 #endif 44 #ifdef PASO_MPI 45 #include 46 #endif 47 48 /*********************************************************************************** 49 50 solve (until convergence) 51 52 (*) [M-theta*dt*L] du = M*u_last + (1-theta)*dt*F(u_last) - M*u + theta*dt F(u) 53 u <- u+du 54 55 with F(u) = L u + f_a(u) (f_a anti-diffusion) 56 and L = D + K + D_K stored in transport_matrix 57 D = Diffusive part (on input stored in transport_matrix) 58 K = flux_matrix 59 D_K = artificial diffusion introduced by K 60 f_a = anti-diffusion introduced by K and u 61 u_last=u of last time step 62 63 For the case theta==0 (explicit case) no iteration is required. One sets 64 65 M*u= M*u_last + dt*F(u_last) 66 67 with F(u) = L u + f_a(u). 68 69 70 For the case theta>0 we set 71 72 A=L-M/(theta*dt) (stored into transport_matrix) 73 F(u)=L u + f_a(u) = M/(theta*dt) u + A u + f_a(u) 74 b(u)= - M/(theta*dt) * u + F(u) = A u + f_a(u) 75 b_last=M/(theta*dt)*u_last + (1-theta)*dt * F(u_last) 76 =M/(theta*dt)*u_last + (1-theta)/theta * [ M/(theta*dt) u_last + A u_last + f_a(u_last) ] 77 =(1/theta**2*dt)*M u_last + (1-theta)/theta*b(u_last) 78 so (*) takes the form 79 80 b_last=(1/theta**2*dt)*M u_last + (1-theta)/theta*b(u_last) 81 while (|du| > tol * |u|) : 82 A*du=b_last + b(u) 83 u <- u-du 84 85 */ 86 87 void Paso_SolverFCT_solve(Paso_FCTransportProblem* fctp, double* u, double dt, double* source, Paso_Options* options) { 88 89 index_t i, j; 90 int n_substeps,n; 91 double dt2=fctp->dt_max, dt2_loc, rtmp,rtmp2,t; 92 dim_t n_rows=Paso_SystemMatrix_getTotalNumRows(fctp->flux_matrix); 93 if (dt<=0.) { 94 Paso_setError(TYPE_ERROR,"Paso_SolverFCT_solve: dt must be positive."); 95 } 96 if (! fctp->valid_matrices) { 97 98 /* extract the row sum of the advective part */ 99 Paso_SystemMatrix_rowSum(fctp->flux_matrix,fctp->row_sum_flux_matrix); 100 /* add the advective part + artificial diffusion to the diffusive part = transport-matrix */ 101 Paso_FCTransportProblem_addAdvectivePart(fctp,1.); 102 /* create a copy of the main diagonal entires of the transport-matrix */ 103 #pragma omp parallel for schedule(static) private(i) 104 for (i=0;itransport_matrix_diagonal[i]= 106 fctp->transport_matrix->mainBlock->val[fctp->main_iptr[i]]; 107 } 108 109 Paso_SystemMatrix_saveMM(fctp->flux_matrix,"flux.mm"); 110 Paso_SystemMatrix_saveMM(fctp->transport_matrix,"trans.mm"); 111 /*=================================================================== * 112 * 113 * calculate time step size: 114 */ 115 dt2=fctp->dt_max; 116 if (fctp->theta < 1.) { 117 dt2=LARGE_POSITIVE_FLOAT; 118 #pragma omp parallel private(dt2_loc) 119 { 120 dt2_loc=LARGE_POSITIVE_FLOAT; 121 #pragma omp for schedule(static) private(i,rtmp,rtmp2) 122 for (i=0;itransport_matrix_diagonal[i]; 124 rtmp2=fctp->lumped_mass_matrix[i]; 125 if ( (rtmp<0 && rtmp2>=0.) || (rtmp>0 && rtmp2<=0.) ) { 126 dt2_loc=MIN(dt2_loc,-rtmp2/rtmp); 127 } 128 } 129 #pragma omp critical 130 { 131 dt2=MIN(dt2,dt2_loc); 132 } 133 } 134 #ifdef PASO_MPI 135 dt2_loc = dt2; 136 MPI_Allreduce(&dt2_loc, &dt2, 1, MPI_DOUBLE, MPI_MIN, fctp->mpi_info->comm); 137 #endif 138 dt2*=1./(1.-fctp->theta); 139 if (fctp->dt_max>0.) dt2=MIN(dt2,fctp->dt_max); 140 } 141 if (dt2 > 0.) { 142 fctp->dt=dt2; 143 } else { 144 fctp->dt=fctp->dt_max; 145 } 146 if (fctp->dt < 0.) { 147 Paso_setError(TYPE_ERROR,"Paso_SolverFCT_solve: dt_max must be positive."); 148 } else { 149 printf("minimum time step size is %e (theta = %e).\n",fctp->dt,fctp->theta); 150 } 151 /* =========================================================== 152 * 153 * 154 */ 155 if (Paso_noError()) { 156 157 } 158 fctp->valid_matrices=Paso_noError(); 159 } 160 161 /* b_last=M*u_last + (1-theta) * F(u_last) */ 162 163 /* decide on substepping */ 164 n_substeps=ceil(dt/fctp->dt); 165 dt2=dt/n_substeps; 166 printf("%d time steps of size is %e (theta = %e, dt_max=%e).\n",n_substeps, dt2,fctp->theta,fctp->dt); 167 168 /* */ 169 Paso_SystemMatrix_allocBuffer(fctp->transport_matrix); 170 171 if (fctp->theta>0) { 172 #pragma omp parallel for schedule(static) private(i) 173 for (i=0;itransport_matrix_diagonal[i]= 175 fctp->transport_matrix->mainBlock->val[fctp->main_iptr[i]]; 176 } 177 } else { 178 /* n_substeps=1; */ 179 /* u= u_last + M^-1*dt*F(u_last) */ 180 t=dt2; 181 n=0; 182 while(nu,u); /* u stores F(u_last)*/ 185 #pragma omp parallel for schedule(static) private(i) 186 for (i=0;iu[i]; 188 fctp->u[i]=fctp->u[i]+dt2*u[i]/fctp->lumped_mass_matrix[i]; 189 printf("%d : %e %e %e\n",i,u[i],fctp->u[i],rtmp); 190 } 191 192 for (i=0;i<21;++i) { 193 for (j=0;j<21;++j) printf("%d->%e ",i*21+j,fctp->u[i*21+j]); 194 printf("\n"); 195 } 196 197 t+=dt2; 198 n++; 199 } 200 } 201 if (Paso_noError()) { 202 #pragma omp parallel for schedule(static) private(i) 203 for (i=0;iu[i]; 205 } 206 } 207 }