/[escript]/branches/symbolic_from_3470/paso/src/FluxLimiter.c
ViewVC logotype

Contents of /branches/symbolic_from_3470/paso/src/FluxLimiter.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3868 - (show annotations)
Thu Mar 15 06:07:08 2012 UTC (7 years, 1 month ago) by caltinay
File MIME type: text/plain
File size: 9747 byte(s)
Update to latest trunk

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

  ViewVC Help
Powered by ViewVC 1.1.26