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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 4286 - (show annotations)
Thu Mar 7 04:28:11 2013 UTC (6 years, 5 months ago) by caltinay
File MIME type: text/plain
File size: 9931 byte(s)
Assorted spelling fixes.

1
2 /*****************************************************************************
3 *
4 * Copyright (c) 2003-2013 by University of Queensland
5 * http://www.uq.edu.au
6 *
7 * Primary Business: Queensland, Australia
8 * Licensed under the Open Software License version 3.0
9 * http://www.opensource.org/licenses/osl-3.0.php
10 *
11 * Development until 2012 by Earth Systems Science Computational Center (ESSCC)
12 * Development since 2012 by School of Earth Sciences
13 *
14 *****************************************************************************/
15
16
17 /************************************************************************************/
18
19 /* Paso: CFT: fluxlimiter for a transport problem
20 *
21 */
22 /************************************************************************************/
23
24 /* Author: l.gross@uq.edu.au */
25
26 /************************************************************************************/
27
28 #include "FluxLimiter.h"
29
30
31 Paso_FCT_FluxLimiter* Paso_FCT_FluxLimiter_alloc(Paso_TransportProblem *fctp)
32 {
33 Paso_FCT_FluxLimiter* out=NULL;
34 const dim_t n=Paso_SystemMatrix_getTotalNumRows(fctp->transport_matrix);
35 const dim_t blockSize=Paso_TransportProblem_getBlockSize(fctp);
36
37 out=MEMALLOC(1, Paso_FCT_FluxLimiter);
38 if (! Esys_checkPtr(out)) {
39
40 out->mpi_info = Esys_MPIInfo_getReference(fctp->mpi_info);
41 out->u_tilde=MEMALLOC(n,double);
42 Esys_checkPtr(out->u_tilde);
43
44 out->MQ=MEMALLOC(2*n,double);
45 Esys_checkPtr(out->MQ);
46
47 out->R=MEMALLOC(2*n,double);
48 Esys_checkPtr(out->R);
49
50 out->R_coupler=Paso_Coupler_alloc(Paso_TransportProblem_borrowConnector(fctp),2*blockSize);
51 out->u_tilde_coupler=Paso_Coupler_alloc(Paso_TransportProblem_borrowConnector(fctp),blockSize);
52 out->antidiffusive_fluxes=Paso_SystemMatrix_alloc(fctp->transport_matrix->type,
53 fctp->transport_matrix->pattern,
54 fctp->transport_matrix->row_block_size,
55 fctp->transport_matrix->col_block_size, TRUE);
56 out->borrowed_lumped_mass_matrix=fctp->lumped_mass_matrix;
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->R_coupler);
75 Paso_Coupler_free(in->u_tilde_coupler);
76 MEMFREE(in);
77 }
78 }
79
80 /* sets the predictor u_tilda from Mu_tilda by solving M_C * u_tilda = Mu_tilda
81 * and calculates the limiters QP and QN: : */
82 void Paso_FCT_FluxLimiter_setU_tilda(Paso_FCT_FluxLimiter* flux_limiter, const double *Mu_tilda) {
83
84 const dim_t n=Paso_FCT_FluxLimiter_getTotalNumRows(flux_limiter);
85 double * remote_u_tilde =NULL;
86 const Paso_SystemMatrixPattern *pattern = Paso_FCT_FluxLimiter_getFluxPattern(flux_limiter);
87 const double *lumped_mass_matrix = flux_limiter->borrowed_lumped_mass_matrix;
88 index_t i, iptr_ij;
89
90 #pragma omp parallel private(i)
91 {
92 #pragma omp for schedule(static)
93 for (i = 0; i < n; ++i) {
94 const double m=lumped_mass_matrix[i];
95 if (m > 0 ) {
96 flux_limiter->u_tilde[i]=Mu_tilda[i]/m;
97 } else {
98 flux_limiter->u_tilde[i]=Mu_tilda[i];
99 }
100
101 }
102 }
103 /* distribute u_tilde: */
104 Paso_Coupler_startCollect(flux_limiter->u_tilde_coupler,flux_limiter->u_tilde);
105 /*
106 * calculate MQ_P[i] = lumped_mass_matrix[i] * max_{j} (\tilde{u}[j]- \tilde{u}[i] )
107 * MQ_N[i] = lumped_mass_matrix[i] * min_{j} (\tilde{u}[j]- \tilde{u}[i] )
108 *
109 */
110 /* first we calculate the min and max of u_tilda in the main block
111 QP, QN are used to hold the result */
112 #pragma omp parallel private(i)
113 {
114 #pragma omp for schedule(static)
115 for (i = 0; i < n; ++i) {
116 if ( flux_limiter->borrowed_lumped_mass_matrix[i]> 0) { /* no constraint */
117 double u_min_i=LARGE_POSITIVE_FLOAT;
118 double u_max_i=-LARGE_POSITIVE_FLOAT;
119 #pragma ivdep
120 for (iptr_ij=(pattern->mainPattern->ptr[i]);iptr_ij<pattern->mainPattern->ptr[i+1]; ++iptr_ij) {
121 const index_t j=pattern->mainPattern->index[iptr_ij];
122 const double u_j=flux_limiter->u_tilde[j];
123 u_min_i=MIN(u_min_i,u_j);
124 u_max_i=MAX(u_max_i,u_j);
125 }
126 flux_limiter->MQ[2*i] = u_min_i;
127 flux_limiter->MQ[2*i+1] = u_max_i;
128
129 } else {
130 flux_limiter->MQ[2*i ] =LARGE_POSITIVE_FLOAT;
131 flux_limiter->MQ[2*i+1] =LARGE_POSITIVE_FLOAT;
132 }
133 }
134 } /* parallel region */
135 /* complete distribute u_tilde: */
136 Paso_Coupler_finishCollect(flux_limiter->u_tilde_coupler);
137 remote_u_tilde=Paso_Coupler_borrowRemoteData(flux_limiter->u_tilde_coupler);
138
139 /* now we look at the couple matrix and set the final value for QP, QN */
140 #pragma omp parallel private(i, iptr_ij)
141 {
142 #pragma omp for schedule(static)
143 for (i = 0; i < n; ++i) {
144 if ( flux_limiter->borrowed_lumped_mass_matrix[i]> 0) { /* no constraint */
145 const double u_i = flux_limiter->u_tilde[i];
146 double u_min_i = flux_limiter->MQ[2*i];
147 double u_max_i = flux_limiter->MQ[2*i+1];
148 #pragma ivdep
149 for (iptr_ij=(pattern->col_couplePattern->ptr[i]);iptr_ij<pattern->col_couplePattern->ptr[i+1]; ++iptr_ij) {
150 const index_t j = pattern->col_couplePattern->index[iptr_ij];
151 const double u_j = remote_u_tilde[j];
152 u_min_i=MIN(u_min_i,u_j);
153 u_max_i=MAX(u_max_i,u_j);
154 }
155 flux_limiter->MQ[2*i ] = ( u_min_i-u_i ) * lumped_mass_matrix[i]; /* M_C*Q_min*/
156 flux_limiter->MQ[2*i+1] = ( u_max_i-u_i ) * lumped_mass_matrix[i] ; /* M_C*Q_max */
157 }
158
159 }
160 } /* parallel region */
161 }
162
163 /* starts to update a vector (not given yet) from the antidiffusion fluxes in flux_limiter->antidiffusive_fluxes
164 (needs u_tilde and Q */
165 void Paso_FCT_FluxLimiter_addLimitedFluxes_Start(Paso_FCT_FluxLimiter* flux_limiter) {
166
167 const dim_t n=Paso_FCT_FluxLimiter_getTotalNumRows(flux_limiter);
168 const Paso_SystemMatrixPattern *pattern = Paso_FCT_FluxLimiter_getFluxPattern(flux_limiter);
169 const double* u_tilde = flux_limiter->u_tilde;
170 const double* remote_u_tilde=Paso_Coupler_borrowRemoteData(flux_limiter->u_tilde_coupler);
171 Paso_SystemMatrix * adf=flux_limiter->antidiffusive_fluxes;
172 index_t iptr_ij;
173 dim_t i;
174
175 #pragma omp parallel private(i, iptr_ij)
176 {
177 #pragma omp for schedule(static)
178 for (i = 0; i < n; ++i) {
179 double R_N_i =1;
180 double R_P_i =1;
181 if ( flux_limiter->borrowed_lumped_mass_matrix[i]> 0) { /* no constraint */
182 const double u_tilde_i=u_tilde[i];
183 double P_P_i=0.;
184 double P_N_i=0.;
185 const double MQ_min=flux_limiter->MQ[2*i];
186 const double MQ_max=flux_limiter->MQ[2*i+1];
187 #pragma ivdep
188 for (iptr_ij=(pattern->mainPattern->ptr[i]);iptr_ij<pattern->mainPattern->ptr[i+1]; ++iptr_ij) {
189 const index_t j=pattern->mainPattern->index[iptr_ij];
190 if (i != j ) {
191 const double f_ij=adf->mainBlock->val[iptr_ij];
192 const double u_tilde_j=u_tilde[j];
193 /* pre-limiter */
194 if (f_ij * (u_tilde_j-u_tilde_i) >= 0) {
195 adf->mainBlock->val[iptr_ij]=0;
196 } else {
197 if (f_ij <=0) {
198 P_N_i+=f_ij;
199 } else {
200 P_P_i+=f_ij;
201 }
202 }
203 }
204 }
205
206 /* now the couple matrix: */
207 #pragma ivdep
208 for (iptr_ij=(pattern->col_couplePattern->ptr[i]);iptr_ij<pattern->col_couplePattern->ptr[i+1]; ++iptr_ij) {
209 const index_t j=pattern->col_couplePattern->index[iptr_ij];
210 const double f_ij=adf->col_coupleBlock->val[iptr_ij];
211 const double u_tilde_j=remote_u_tilde[j];
212 /* pre-limiter */
213 if (f_ij * (u_tilde_j-u_tilde_i) >= 0) {
214 adf->col_coupleBlock->val[iptr_ij]=0;
215 } else {
216 if (f_ij <=0) {
217 P_N_i+=f_ij;
218 } else {
219 P_P_i+=f_ij;
220 }
221 }
222 }
223 /* finally the R+ and R- are calculated */
224
225 if (P_N_i<0) R_N_i=MIN(1,MQ_min/P_N_i);
226 if (P_P_i>0) R_P_i=MIN(1,MQ_max/P_P_i);
227 }
228 flux_limiter->R[2*i] = R_N_i;
229 flux_limiter->R[2*i+1] = R_P_i;
230 }
231 } /* end parallel region */
232
233 /* now we kick off the distribution of the R's */
234 Paso_Coupler_startCollect(flux_limiter->R_coupler,flux_limiter->R);
235 }
236
237
238 /* completes the exchange of the R factors and adds the weighted anti-diffusion fluxes to the residual b */
239 void Paso_FCT_FluxLimiter_addLimitedFluxes_Complete(Paso_FCT_FluxLimiter* flux_limiter, double* b)
240 {
241 const dim_t n=Paso_FCT_FluxLimiter_getTotalNumRows(flux_limiter);
242 const Paso_SystemMatrixPattern *pattern = Paso_FCT_FluxLimiter_getFluxPattern(flux_limiter);
243 const Paso_SystemMatrix * adf=flux_limiter->antidiffusive_fluxes;
244 double *R = flux_limiter->R;
245 double *remote_R=NULL;
246 dim_t i;
247 index_t iptr_ij;
248
249 remote_R=Paso_Coupler_finishCollect(flux_limiter->R_coupler);
250
251 #pragma omp parallel for schedule(static) private(i, iptr_ij)
252 for (i = 0; i < n; ++i) {
253
254 const double R_N_i=R[2*i];
255 const double R_P_i=R[2*i+1];
256 double f_i=b[i];
257
258 #pragma ivdep
259 for (iptr_ij=(pattern->mainPattern->ptr[i]);iptr_ij<pattern->mainPattern->ptr[i+1]; ++iptr_ij) {
260 const index_t j = pattern->mainPattern->index[iptr_ij];
261 const double f_ij = adf->mainBlock->val[iptr_ij];
262 const double R_P_j = R[2*j+1];
263 const double R_N_j = R[2*j];
264 const double rtmp=((f_ij >=0 ) ? MIN(R_P_i,R_N_j) : MIN(R_N_i,R_P_j) ) ;
265 f_i+=f_ij*rtmp;
266 }
267 #pragma ivdep
268 for (iptr_ij=(pattern->col_couplePattern->ptr[i]);iptr_ij<pattern->col_couplePattern->ptr[i+1]; ++iptr_ij) {
269 const index_t j = pattern->col_couplePattern->index[iptr_ij];
270 const double f_ij = adf->col_coupleBlock->val[iptr_ij];
271 const double R_P_j = remote_R[2*j+1];
272 const double R_N_j = remote_R[2*j];
273 const double rtmp=(f_ij >=0 ) ? MIN(R_P_i, R_N_j) : MIN(R_N_i, R_P_j) ;
274 f_i+=f_ij*rtmp;
275 }
276 b[i]=f_i;
277 } /* and of i loop */
278 }

  ViewVC Help
Powered by ViewVC 1.1.26