/[escript]/trunk/paso/src/FluxLimiter.c
ViewVC logotype

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3793 - (show annotations)
Wed Feb 1 07:39:43 2012 UTC (7 years, 8 months ago) by gross
File MIME type: text/plain
File size: 10499 byte(s)
new implementation of FCT solver with some modifications to the python interface
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: CFT: fluxlimiter for a transport problem
18 *
19 */
20 /**************************************************************/
21
22 /* Author: l.gross@uq.edu.au */
23
24 /**************************************************************/
25
26 #include "FluxLimiter.h"
27
28
29 Paso_FCT_FluxLimiter* Paso_FCT_FluxLimiter_alloc(Paso_TransportProblem *fctp)
30 {
31 Paso_FCT_FluxLimiter* out=NULL;
32 const dim_t n=Paso_SystemMatrix_getTotalNumRows(fctp->transport_matrix);
33 const dim_t blockSize=Paso_TransportProblem_getBlockSize(fctp);
34
35 out=MEMALLOC(1, Paso_FCT_FluxLimiter);
36 if (! Esys_checkPtr(out)) {
37
38 out->mpi_info = Esys_MPIInfo_getReference(fctp->mpi_info);
39 out->u_tilde=MEMALLOC(n,double);
40 Esys_checkPtr(out->u_tilde);
41
42 out->MQ=MEMALLOC(2*n,double);
43 Esys_checkPtr(out->MQ);
44
45 out->R=MEMALLOC(2*n,double);
46 Esys_checkPtr(out->R);
47
48 out->MQ_coupler=Paso_Coupler_alloc(Paso_TransportProblem_borrowConnector(fctp),2*blockSize);
49 out->R_coupler=Paso_Coupler_alloc(Paso_TransportProblem_borrowConnector(fctp),2*blockSize);
50 out->u_tilde_coupler=Paso_Coupler_alloc(Paso_TransportProblem_borrowConnector(fctp),blockSize);
51 out->antidiffusive_fluxes=Paso_SystemMatrix_alloc(fctp->transport_matrix->type,
52 fctp->transport_matrix->pattern,
53 fctp->transport_matrix->row_block_size,
54 fctp->transport_matrix->col_block_size, TRUE);
55 out->borrowed_lumped_mass_matrix=fctp->lumped_mass_matrix;
56
57 }
58 if (Esys_noError()) {
59 return out;
60 } else {
61 Paso_FCT_FluxLimiter_free(out);
62 return NULL;
63 }
64 }
65
66 void Paso_FCT_FluxLimiter_free(Paso_FCT_FluxLimiter * in)
67 {
68 if (in!=NULL) {
69 Esys_MPIInfo_free(in->mpi_info);
70 Paso_SystemMatrix_free(in->antidiffusive_fluxes);
71 MEMFREE(in->u_tilde);
72 MEMFREE(in->MQ);
73 MEMFREE(in->R);
74 Paso_Coupler_free(in->MQ_coupler);
75 Paso_Coupler_free(in->R_coupler);
76 Paso_Coupler_free(in->u_tilde_coupler);
77 MEMFREE(in);
78 }
79 }
80
81 /* sets the predictor u_tilda from Mu_tilda by solving M_C * u_tilda = Mu_tilda
82 * and calculates the limiters QP and QN: : */
83 void Paso_FCT_FluxLimiter_setU_tilda(Paso_FCT_FluxLimiter* flux_limiter, const double *Mu_tilda) {
84
85 const dim_t n=Paso_FCT_FluxLimiter_getTotalNumRows(flux_limiter);
86 double * remote_u_tilde =NULL;
87 const Paso_SystemMatrixPattern *pattern = Paso_FCT_FluxLimiter_getFluxPattern(flux_limiter);
88 const double *lumped_mass_matrix = flux_limiter->borrowed_lumped_mass_matrix;
89 index_t i, iptr_ij;
90
91 #pragma omp parallel private(i, iptr_ij)
92 {
93 #pragma omp for schedule(static)
94 for (i = 0; i < n; ++i) {
95 const register double m=lumped_mass_matrix[i];
96 flux_limiter->u_tilde[i]=Mu_tilda[i]/m;
97 }
98
99 /* distribute u_tilde: */
100 #pragma omp master
101 {
102 Paso_Coupler_startCollect(flux_limiter->u_tilde_coupler,flux_limiter->u_tilde);
103 }
104 /*
105 * calculate MQ_P[i] = lumped_mass_matrix[i] * max_{j} (\tilde{u}[j]- \tilde{u}[i] )
106 * MQ_N[i] = lumped_mass_matrix[i] * min_{j} (\tilde{u}[j]- \tilde{u}[i] )
107 *
108 */
109 /* first we calculate the min and max of u_tilda in the main block
110 QP, QN are used to hold the result */
111
112 #pragma omp for schedule(static)
113 for (i = 0; i < n; ++i) {
114 register double u_min_i=LARGE_POSITIVE_FLOAT;
115 register double u_max_i=-LARGE_POSITIVE_FLOAT;
116 #pragma ivdep
117 for (iptr_ij=(pattern->mainPattern->ptr[i]);iptr_ij<pattern->mainPattern->ptr[i+1]; ++iptr_ij) {
118 const register index_t j=pattern->mainPattern->index[iptr_ij];
119 const register double u_j=flux_limiter->u_tilde[j];
120 u_min_i=MIN(u_min_i,u_j);
121 u_max_i=MAX(u_max_i,u_j);
122 }
123 flux_limiter->MQ[2*i] = u_min_i;
124 flux_limiter->MQ[2*i+1] = u_max_i;
125 }
126 /* complete distribute u_tilde: */
127 #pragma omp master
128 {
129 Paso_Coupler_finishCollect(flux_limiter->u_tilde_coupler);
130 remote_u_tilde=Paso_Coupler_borrowRemoteData(flux_limiter->u_tilde_coupler);
131 }
132 /* now we look at the couple matrix and set the final value for QP, QN */
133 #pragma omp for schedule(static)
134 for (i = 0; i < n; ++i) {
135 const register double u_i = flux_limiter->u_tilde[i];
136 register double u_min_i = flux_limiter->MQ[2*i];
137 register double u_max_i = flux_limiter->MQ[2*i+1];
138 #pragma ivdep
139 for (iptr_ij=(pattern->col_couplePattern->ptr[i]);iptr_ij<pattern->col_couplePattern->ptr[i+1]; ++iptr_ij) {
140 const register index_t j = pattern->col_couplePattern->index[iptr_ij];
141 const register double u_j = remote_u_tilde[j];
142 u_min_i=MIN(u_min_i,u_j);
143 u_max_i=MAX(u_max_i,u_j);
144 }
145 flux_limiter->MQ[2*i ] = ( u_min_i-u_i ) * lumped_mass_matrix[i]; /* M_C*Q_min*/
146 flux_limiter->MQ[2*i+1] = ( u_max_i-u_i ) * lumped_mass_matrix[i] ; /* M_C*Q_max */
147 }
148 } /* parallel region */
149
150
151
152 Paso_Coupler_startCollect(flux_limiter->MQ_coupler,flux_limiter->MQ);
153 Paso_Coupler_finishCollect(flux_limiter->MQ_coupler);
154 }
155
156 /* starts to update a vector (not given yet) from the antidiffusion fluxes in flux_limiter->antidiffusive_fluxes
157 (needes u_tilde and Q */
158 void Paso_FCT_FluxLimiter_addLimitedFluxes_Start(Paso_FCT_FluxLimiter* flux_limiter) {
159
160 const dim_t n=Paso_FCT_FluxLimiter_getTotalNumRows(flux_limiter);
161 const Paso_SystemMatrixPattern *pattern = Paso_FCT_FluxLimiter_getFluxPattern(flux_limiter);
162 const double* u_tilde = flux_limiter->u_tilde;
163 const double* remote_u_tilde=Paso_Coupler_borrowRemoteData(flux_limiter->u_tilde_coupler);
164 Paso_SystemMatrix * adf=flux_limiter->antidiffusive_fluxes;
165 /* const double *lumped_mass_matrix = flux_limiter->borrowed_lumped_mass_matrix; */
166 register index_t iptr_ij;
167 dim_t i;
168
169 #pragma omp parallel private(i, iptr_ij)
170 {
171 #pragma omp for schedule(static)
172 for (i = 0; i < n; ++i) {
173
174 const register double u_tilde_i=u_tilde[i];
175 register double P_P_i=0.;
176 register double P_N_i=0.;
177 const register double MQ_min=flux_limiter->MQ[2*i];
178 const register double MQ_max=flux_limiter->MQ[2*i+1];
179 register double R_N_i =1;
180 register double R_P_i =1;
181 #pragma ivdep
182 for (iptr_ij=(pattern->mainPattern->ptr[i]);iptr_ij<pattern->mainPattern->ptr[i+1]; ++iptr_ij) {
183 const register index_t j=pattern->mainPattern->index[iptr_ij];
184 if (i != j ) {
185 const register double f_ij=adf->mainBlock->val[iptr_ij];
186 const register double u_tilde_j=u_tilde[j];
187 /* pre-limiter */
188 if (f_ij * (u_tilde_j-u_tilde_i) >= 0) {
189 adf->mainBlock->val[iptr_ij]=0;
190 } else {
191 if (f_ij <=0) {
192 P_N_i+=f_ij;
193 } else {
194 P_P_i+=f_ij;
195 }
196 }
197 }
198 }
199
200 /* now the couple matrix: */
201 #pragma ivdep
202 for (iptr_ij=(pattern->col_couplePattern->ptr[i]);iptr_ij<pattern->col_couplePattern->ptr[i+1]; ++iptr_ij) {
203 const register index_t j=pattern->col_couplePattern->index[iptr_ij];
204 const register double f_ij=adf->col_coupleBlock->val[iptr_ij];
205 const register double u_tilde_j=remote_u_tilde[j];
206 /* pre-limiter */
207 if (f_ij * (u_tilde_j-u_tilde_i) >= 0) {
208 adf->col_coupleBlock->val[iptr_ij]=0;
209 } else {
210 if (f_ij <=0) {
211 P_N_i+=f_ij;
212 } else {
213 P_P_i+=f_ij;
214 }
215 }
216 }
217 /* finally the R+ and R- are calculated */
218
219 if (P_N_i<0) R_N_i=MIN(1,MQ_min/P_N_i);
220 if (P_P_i>0) R_P_i=MIN(1,MQ_max/P_P_i);
221 flux_limiter->R[2*i] = R_N_i;
222 flux_limiter->R[2*i+1] = R_P_i;
223 }
224 } /* end parallel region */
225
226 /* now we kick off the distrribution of the R's */
227 Paso_Coupler_startCollect(flux_limiter->R_coupler,flux_limiter->R);
228 }
229
230
231 /* completes the exchange of the R factors and adds the weighted antdiffusion fluxxes to the residual b */
232 void Paso_FCT_FluxLimiter_addLimitedFluxes_Complete(Paso_FCT_FluxLimiter* flux_limiter, double* b)
233 {
234 const dim_t n=Paso_FCT_FluxLimiter_getTotalNumRows(flux_limiter);
235 const Paso_SystemMatrixPattern *pattern = Paso_FCT_FluxLimiter_getFluxPattern(flux_limiter);
236 const Paso_SystemMatrix * adf=flux_limiter->antidiffusive_fluxes;
237 double *R = flux_limiter->R;
238 double *remote_R=NULL;
239 register dim_t i;
240 register index_t iptr_ij;
241
242 remote_R=Paso_Coupler_finishCollect(flux_limiter->R_coupler);
243
244 #pragma omp parallel for schedule(static) private(i, iptr_ij)
245 for (i = 0; i < n; ++i) {
246
247 const register double R_N_i=R[2*i];
248 const register double R_P_i=R[2*i+1];
249 register double f_i=b[i];
250
251 #pragma ivdep
252 for (iptr_ij=(pattern->mainPattern->ptr[i]);iptr_ij<pattern->mainPattern->ptr[i+1]; ++iptr_ij) {
253 const register index_t j = pattern->mainPattern->index[iptr_ij];
254 const register double f_ij = adf->mainBlock->val[iptr_ij];
255 const register double R_P_j = R[2*j+1];
256 const register double R_N_j = R[2*j];
257 const register double rtmp=((f_ij >=0 ) ? MIN(R_P_i,R_N_j) : MIN(R_N_i,R_P_j) ) ;
258 f_i+=f_ij*rtmp;
259 }
260 #pragma ivdep
261 for (iptr_ij=(pattern->col_couplePattern->ptr[i]);iptr_ij<pattern->col_couplePattern->ptr[i+1]; ++iptr_ij) {
262 const register index_t j = pattern->col_couplePattern->index[iptr_ij];
263 const register double f_ij = adf->col_coupleBlock->val[iptr_ij];
264 const register double R_P_j = remote_R[2*j+1];
265 const register double R_N_j = remote_R[2*j];
266 const register double rtmp=(f_ij >=0 ) ? MIN(R_P_i, R_N_j) : MIN(R_N_i, R_P_j) ;
267 f_i+=f_ij*rtmp;
268 }
269 b[i]=f_i;
270 } /* and of i loop */
271 }

  ViewVC Help
Powered by ViewVC 1.1.26