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

Annotation of /trunk/paso/src/SolverFCT_FluxControl.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1401 - (hide annotations)
Fri Jan 25 04:31:18 2008 UTC (12 years, 1 month ago) by gross
File MIME type: text/plain
File size: 8328 byte(s)
rewrite antidiffusion calculation to avoid coloring for OPENMP parallelization
1 gross 1361 /* $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: FluxControl */
17    
18     /**************************************************************/
19    
20     /* Author: l.gross@uq.edu.au */
21    
22     /**************************************************************/
23    
24    
25     #include "Paso.h"
26     #include "SolverFCT.h"
27     #include "PasoUtil.h"
28    
29     #define FLUX_S(a,b) ((SIGN(a)+SIGN(b))/2.)
30     #define MINMOD(a,b) (FLUX_S(a,b)*MIN(ABS(a),ABS(b)))
31     #define SUPERBEE(a,b) (FLUX_S(a,b)*MAX(MIN(2*ABS(a),ABS(b)),MIN(ABS(a),2*ABS(b))))
32 gross 1371 #define MC(a,b) (FLUX_S(a,b)*MIN3(ABS((a)+(b))/2,2*ABS(a),2*ABS(b)))
33 gross 1361
34 gross 1371 #define FLUX_L(a,b) MC(a,b) /* alter for other flux limiter */
35 gross 1361
36 gross 1371 #define FLUX_LIMITER(a) FLUX_L(1,a)
37 gross 1361
38     /**************************************************************/
39    
40     /* adds A plus stabelising diffusion into the matrix B */
41    
42     /* d_ij=alpha*max(0,-a[i,j],-a[j,i]) */
43     /* b[i,j]+=alpha*(a[i,j]+d_ij) */
44     /* b[j,i]+=alpha*(a[j,i]+d_ij) */
45     /* b[i,i]-=alpha*d_ij */
46     /* b[j,j]-=alpha*d_ij */
47    
48 gross 1363 void Paso_FCTransportProblem_addAdvectivePart(Paso_FCTransportProblem * fc, double alpha) {
49 gross 1361 dim_t n,i;
50     index_t color, iptr_ij,j,iptr_ji;
51 gross 1401 register double d_ij, sum;
52 gross 1361
53     if (fc==NULL) return;
54 gross 1362 n=Paso_SystemMatrix_getTotalNumRows(fc->flux_matrix);
55 gross 1361
56 gross 1401 #pragma omp parallel for private(i,iptr_ij,j,iptr_ji,d_ij,sum) schedule(static)
57     for (i = 0; i < n; ++i) {
58     sum=alpha*fc->flux_matrix->mainBlock->val[fc->main_iptr[i]];
59     for (iptr_ij=fc->flux_matrix->mainBlock->pattern->ptr[i];iptr_ij<fc->flux_matrix->mainBlock->pattern->ptr[i+1]; ++iptr_ij) {
60     j=fc->flux_matrix->mainBlock->pattern->index[iptr_ij];
61     if (j!=i) {
62     /* find entry a[j,i] */
63     for (iptr_ji=fc->flux_matrix->mainBlock->pattern->ptr[j]; iptr_ji<fc->flux_matrix->mainBlock->pattern->ptr[j+1]; ++iptr_ji) {
64     if (fc->flux_matrix->mainBlock->pattern->index[iptr_ji]==i) {
65     d_ij=(-alpha)*MIN3(0.,fc->flux_matrix->mainBlock->val[iptr_ij],
66     fc->flux_matrix->mainBlock->val[iptr_ji]);
67     fc->transport_matrix->mainBlock->val[iptr_ij]+=
68     alpha*fc->flux_matrix->mainBlock->val[iptr_ij]+d_ij;
69     sum-=d_ij;
70     break;
71     }
72     }
73     }
74     }
75     /* TODO process couple block */
76 gross 1370
77 gross 1401 /* update main diagonal */
78     fc->transport_matrix->mainBlock->val[fc->main_iptr[i]]+=sum;
79 gross 1361 }
80    
81     }
82 gross 1370
83     void Paso_FCTransportProblem_setFlux(Paso_FCTransportProblem * fc, double * u, double* fa) {
84 gross 1400 /*
85     * sets fa=transport_matrix*u+anti-diffuison_flux(u)
86     */
87 gross 1370
88     double *remote_u=NULL;
89    
90     if (fc==NULL) return;
91     Paso_SystemMatrix_startCollect(fc->transport_matrix,u);
92     /* process main block */
93     Paso_SparseMatrix_MatrixVector_CSR_OFFSET0(1.,fc->transport_matrix->mainBlock,u,0.,fa);
94     /* finish exchange */
95     remote_u=Paso_SystemMatrix_finishCollect(fc->transport_matrix);
96     /* process couple block */
97     Paso_SparseMatrix_MatrixVector_CSR_OFFSET0(1.,fc->transport_matrix->coupleBlock,remote_u,1.,fa);
98    
99     Paso_FCTransportProblem_setAntiDiffusiveFlux(fc,u,remote_u,fa);
100     }
101 gross 1361 /**************************************************************/
102    
103     /* adds antidiffusion to fa
104    
105     P_p[i] + = sum_j min(0,a[i,j]) min(0,u[j]-u[i])
106     P_n[i] + = sum_j min(0,a[i,j]) max(0,u[j]-u[i])
107     Q_p[i] + = sum_j max(0,a[i,j]) max(0,u[j]-u[i])
108     Q_n[i] + = sum_j max(0,a[i,j]) min(0,u[j]-u[i])
109     d_ij=max(0,-a[i,j],-a[j,i])
110     l_ji=max(0,a[j,i],a[j,i]-a[i,j])
111     if a[j,i] >= a[i,j] and 0>a[i,j] : (i.e d_ij>0 and l_ji>=l_ij)
112     r_ij = u[i]>u[j] ? Q_p[i]/P_p[i] : Q_n[i]/Q_n[i]
113     f_ij =min(FLUX_LIMITER(r_ij)*d_ij,l_ji) (u[i]-u[j])=min(FLUX_LIMITER(r_ij)*a[i,j],a[j,i]-a[i,j]) (u[i]-u[j])
114     fa[i]+=f_ij
115     fa[j]-=f_ij
116    
117     */
118    
119 gross 1370 void Paso_FCTransportProblem_setAntiDiffusiveFlux(Paso_FCTransportProblem * fc, double * u, double *u_remote, double* fa) {
120 gross 1361
121 gross 1401 register double u_i,P_p,P_n,Q_p,Q_n,r_p,r_n, a_ij, d, u_j, r_ij, f_ij, a_ji, d_ij, sum;
122 gross 1361 index_t color, iptr_ij,j,iptr_ji, i;
123     dim_t n;
124    
125    
126     if (fc==NULL) return;
127 gross 1362 n=Paso_SystemMatrix_getTotalNumRows(fc->flux_matrix);
128 gross 1361
129 gross 1370
130 gross 1401 #pragma omp parallel
131 gross 1361 {
132 gross 1401 /*
133     * calculate the smootness sensors
134     */
135     #pragma omp for schedule(static) private(i, u_i,P_p,P_n,Q_p,Q_n,iptr_ij,a_ij,j,d)
136     for (i = 0; i < n; ++i) {
137     u_i=u[i];
138     P_p=0.;
139     P_n=0.;
140     Q_p=0.;
141     Q_n=0.;
142     #pragma ivdep
143     for (iptr_ij=(fc->flux_matrix->mainBlock->pattern->ptr[i]);iptr_ij<(fc->flux_matrix->mainBlock->pattern->ptr[i+1]); ++iptr_ij) {
144     a_ij=fc->flux_matrix->mainBlock->val[iptr_ij];
145     j=fc->flux_matrix->mainBlock->pattern->index[iptr_ij];
146     d=u[j]-u_i;
147     if (a_ij<0.) {
148     if (d<0.) {
149     P_p+=a_ij*d;
150     } else {
151     P_n+=a_ij*d;
152     }
153     } else {
154     if (d>0.) {
155     Q_p+=a_ij*d;
156     } else {
157     Q_n+=a_ij*d;
158     }
159     }
160     }
161     #pragma ivdep
162     for (iptr_ij=(fc->flux_matrix->coupleBlock->pattern->ptr[i]);iptr_ij<(fc->flux_matrix->coupleBlock->pattern->ptr[i+1]); ++iptr_ij) {
163     a_ij=fc->flux_matrix->coupleBlock->val[iptr_ij];
164     j=fc->flux_matrix->coupleBlock->pattern->index[iptr_ij];
165     d=u_remote[j]-u_i;
166     if (a_ij<0.) {
167     if (d<0.) {
168     P_p+=a_ij*d;
169     } else {
170     P_n+=a_ij*d;
171     }
172     } else {
173     if (d>0.) {
174     Q_p+=a_ij*d;
175     } else {
176     Q_n+=a_ij*d;
177     }
178     }
179     }
180     /* set the smoothness indicators */
181     fc->r_p[i] = (P_p > 0.) ? FLUX_LIMITER(Q_p/P_p) : 0.;
182     fc->r_n[i] = (P_n < 0.) ? FLUX_LIMITER(Q_n/P_n) : 0.;
183 gross 1371
184 gross 1401 } /* end of row loop for smootheness indicator */
185    
186     /*
187     * calculate antidiffusion
188     */
189     #pragma omp for schedule(static) private(i, u_i, sum, iptr_ij, a_ij, j, iptr_ji, a_ji,d_ij, u_j, r_ij, f_ij)
190     for (i = 0; i < n; ++i) {
191     u_i=u[i];
192     sum=0;
193     /* anti diffusive flux from main block */
194     for (iptr_ij=fc->flux_matrix->mainBlock->pattern->ptr[i];iptr_ij<fc->flux_matrix->mainBlock->pattern->ptr[i+1]; ++iptr_ij) {
195     a_ij=fc->flux_matrix->mainBlock->val[iptr_ij];
196     j=fc->flux_matrix->mainBlock->pattern->index[iptr_ij];
197     if ( i!=j ) {
198     /* find entry a[j,i] */
199     for (iptr_ji=fc->flux_matrix->mainBlock->pattern->ptr[j];iptr_ji<fc->flux_matrix->mainBlock->pattern->ptr[j+1]; ++iptr_ji) {
200     if (fc->flux_matrix->mainBlock->pattern->index[iptr_ji]==i) {
201     a_ji=fc->flux_matrix->mainBlock->val[iptr_ji];
202     d_ij=-MIN3(0,a_ji,a_ij);
203     if (d_ij > 0.) {
204     u_j=u[j];
205     if (a_ji >= a_ij) {
206     r_ij = u_i>u_j ? fc->r_p[i] : fc->r_n[i];
207     f_ij =MIN(r_ij*d_ij,a_ji+d_ij);
208     } else {
209     r_ij = u_j>u_i ? fc->r_p[j] : fc->r_n[j];
210     f_ij =MIN(r_ij*d_ij,a_ij+d_ij);
211     }
212     sum+=f_ij*(u_i-u_j);
213     }
214     break;
215 gross 1361 }
216 gross 1401 }
217     }
218     }
219     /* anti diffusive flux from couple block */
220     /* TODO */
221 gross 1371
222    
223 gross 1375
224 gross 1401 fa[i]+=sum;
225     }
226     } /* end of parallel block */
227 gross 1361 }

  ViewVC Help
Powered by ViewVC 1.1.26