/[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 1371 - (hide annotations)
Thu Jan 3 06:11:21 2008 UTC (11 years, 7 months ago) by gross
File MIME type: text/plain
File size: 9534 byte(s)
some bugs in the updwing scheme but there are still problem with the cross wind direction.
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     register double d_ij;
52    
53     if (fc==NULL) return;
54 gross 1362 n=Paso_SystemMatrix_getTotalNumRows(fc->flux_matrix);
55 gross 1361
56     #pragma omp parallel private(color)
57     {
58     /* process main block */
59     for (color=0;color<fc->num_colors;++color) {
60     #pragma omp for private(i,iptr_ij,j,iptr_ji,d_ij) schedule(static)
61     for (i = 0; i < n; ++i) {
62     if (fc->colorOf[i]==color) {
63 gross 1370 fc->transport_matrix->mainBlock->val[fc->main_iptr[i]]+=alpha*fc->flux_matrix->mainBlock->val[fc->main_iptr[i]];
64    
65 gross 1362 for (iptr_ij=fc->flux_matrix->mainBlock->pattern->ptr[i];iptr_ij<fc->flux_matrix->mainBlock->pattern->ptr[i+1]; ++iptr_ij) {
66     j=fc->flux_matrix->mainBlock->pattern->index[iptr_ij];
67 gross 1370 if (j<i) {
68 gross 1361 /* find entry a[j,i] */
69 gross 1370 for (iptr_ji=fc->flux_matrix->mainBlock->pattern->ptr[j]; iptr_ji<fc->flux_matrix->mainBlock->pattern->ptr[j+1]; ++iptr_ji) {
70 gross 1362 if (fc->flux_matrix->mainBlock->pattern->index[iptr_ji]==i) {
71     d_ij=(-alpha)*MIN3(0.,fc->flux_matrix->mainBlock->val[iptr_ij],fc->flux_matrix->mainBlock->val[iptr_ji]);
72 gross 1370 printf("%d %d -> %e\n",i,j,d_ij);
73 gross 1363 fc->transport_matrix->mainBlock->val[iptr_ij]+=alpha*fc->flux_matrix->mainBlock->val[iptr_ij]+d_ij;
74     fc->transport_matrix->mainBlock->val[iptr_ji]+=alpha*fc->flux_matrix->mainBlock->val[iptr_ji]+d_ij;
75 gross 1370 printf("%d %d -> %e -> %e %e \n",i,j,d_ij,fc->transport_matrix->mainBlock->val[iptr_ij],fc->transport_matrix->mainBlock->val[iptr_ji]);
76 gross 1363 fc->transport_matrix->mainBlock->val[fc->main_iptr[i]]-=d_ij;
77     fc->transport_matrix->mainBlock->val[fc->main_iptr[j]]-=d_ij;
78 gross 1361 break;
79     }
80     }
81     }
82    
83     }
84     /* TODO process couple block */
85     }
86     }
87     #pragma omp barrier
88     }
89     }
90    
91     }
92 gross 1370
93     void Paso_FCTransportProblem_setFlux(Paso_FCTransportProblem * fc, double * u, double* fa) {
94    
95     double *remote_u=NULL;
96    
97     if (fc==NULL) return;
98     Paso_SystemMatrix_startCollect(fc->transport_matrix,u);
99     /* process main block */
100     Paso_SparseMatrix_MatrixVector_CSR_OFFSET0(1.,fc->transport_matrix->mainBlock,u,0.,fa);
101     /* finish exchange */
102     remote_u=Paso_SystemMatrix_finishCollect(fc->transport_matrix);
103     /* process couple block */
104     Paso_SparseMatrix_MatrixVector_CSR_OFFSET0(1.,fc->transport_matrix->coupleBlock,remote_u,1.,fa);
105    
106     Paso_FCTransportProblem_setAntiDiffusiveFlux(fc,u,remote_u,fa);
107     }
108 gross 1361 /**************************************************************/
109    
110     /* adds antidiffusion to fa
111    
112     P_p[i] + = sum_j min(0,a[i,j]) min(0,u[j]-u[i])
113     P_n[i] + = sum_j min(0,a[i,j]) max(0,u[j]-u[i])
114     Q_p[i] + = sum_j max(0,a[i,j]) max(0,u[j]-u[i])
115     Q_n[i] + = sum_j max(0,a[i,j]) min(0,u[j]-u[i])
116     d_ij=max(0,-a[i,j],-a[j,i])
117     l_ji=max(0,a[j,i],a[j,i]-a[i,j])
118     if a[j,i] >= a[i,j] and 0>a[i,j] : (i.e d_ij>0 and l_ji>=l_ij)
119     r_ij = u[i]>u[j] ? Q_p[i]/P_p[i] : Q_n[i]/Q_n[i]
120     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])
121     fa[i]+=f_ij
122     fa[j]-=f_ij
123    
124     */
125    
126 gross 1370 void Paso_FCTransportProblem_setAntiDiffusiveFlux(Paso_FCTransportProblem * fc, double * u, double *u_remote, double* fa) {
127 gross 1361
128 gross 1371 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;
129 gross 1361 index_t color, iptr_ij,j,iptr_ji, i;
130     dim_t n;
131    
132    
133     if (fc==NULL) return;
134 gross 1362 n=Paso_SystemMatrix_getTotalNumRows(fc->flux_matrix);
135 gross 1361 /* exchange */
136    
137 gross 1370
138 gross 1361 #pragma omp parallel private(color)
139     {
140     for (color=0;color<fc->num_colors;++color) {
141 gross 1371 #pragma omp for schedule(static) private(i, u_i,P_p,P_n,Q_p,Q_n,r_p,r_n,iptr_ij,a_ij,d,j,iptr_ji, u_j, r_ij, f_ij, a_ji, d_ij)
142 gross 1361 for (i = 0; i < n; ++i) {
143     if (fc->colorOf[i]==color) {
144     u_i=u[i];
145     /* gather the smoothness sensor */
146     P_p=0.;
147     P_n=0.;
148     Q_p=0.;
149     Q_n=0.;
150 gross 1371 r_p=0.;
151     r_n=0.;
152     /* #pragma ivdep */
153 gross 1362 for (iptr_ij=(fc->flux_matrix->mainBlock->pattern->ptr[i]);iptr_ij<(fc->flux_matrix->mainBlock->pattern->ptr[i+1]); ++iptr_ij) {
154     a_ij=fc->flux_matrix->mainBlock->val[iptr_ij];
155     j=fc->flux_matrix->mainBlock->pattern->index[iptr_ij];
156 gross 1361 d=u[j]-u_i;
157 gross 1371 printf("%d %d : %e %e :: %e \n",i,j,u_i,u[j],a_ij);
158 gross 1361 if (a_ij<0.) {
159     if (d<0.) {
160     P_p+=a_ij*d;
161     } else {
162     P_n+=a_ij*d;
163     }
164     } else {
165     if (d>0.) {
166     Q_p+=a_ij*d;
167     } else {
168     Q_n+=a_ij*d;
169     }
170     }
171     }
172     #pragma ivdep
173 gross 1362 for (iptr_ij=(fc->flux_matrix->coupleBlock->pattern->ptr[i]);iptr_ij<(fc->flux_matrix->coupleBlock->pattern->ptr[i+1]); ++iptr_ij) {
174     a_ij=fc->flux_matrix->coupleBlock->val[iptr_ij];
175     j=fc->flux_matrix->coupleBlock->pattern->index[iptr_ij];
176 gross 1361 d=u_remote[j]-u_i;
177 gross 1371
178 gross 1361 if (a_ij<0.) {
179     if (d<0.) {
180     P_p+=a_ij*d;
181     } else {
182     P_n+=a_ij*d;
183     }
184     } else {
185     if (d>0.) {
186     Q_p+=a_ij*d;
187     } else {
188     Q_n+=a_ij*d;
189     }
190     }
191     }
192     /* set the smoothness indicators */
193     r_p = (P_p > 0.) ? FLUX_LIMITER(Q_p/P_p) : 0.;
194 gross 1371 r_n = (P_n < 0.) ? FLUX_LIMITER(Q_n/P_n) : 0.;
195     printf("%d: %e %e %e : %e %e %e : %e\n",i,Q_p,P_p,r_p,Q_n,P_n,r_n,u_i);
196 gross 1361 /* anti diffusive flux from main block */
197 gross 1362 for (iptr_ij=fc->flux_matrix->mainBlock->pattern->ptr[i];iptr_ij<fc->flux_matrix->mainBlock->pattern->ptr[i+1]; ++iptr_ij) {
198     a_ij=fc->flux_matrix->mainBlock->val[iptr_ij];
199     j=fc->flux_matrix->mainBlock->pattern->index[iptr_ij];
200 gross 1371 if ( i!=j ) {
201 gross 1361 /* find entry a[j,i] */
202 gross 1371 for (iptr_ji=fc->flux_matrix->mainBlock->pattern->ptr[j];iptr_ji<fc->flux_matrix->mainBlock->pattern->ptr[j+1]; ++iptr_ji) {
203 gross 1362 if (fc->flux_matrix->mainBlock->pattern->index[iptr_ji]==i) {
204     a_ji=fc->flux_matrix->mainBlock->val[iptr_ji];
205 gross 1371 d_ij=-MIN3(0,a_ji,a_ij);
206     if ( (d_ij > 0.) && ( (a_ji > a_ij) || ( (a_ji == a_ij) && (j<i) ) ) ) {
207 gross 1361 u_j=u[j];
208     r_ij = u_i>u_j ? r_p : r_n;
209 gross 1371 f_ij =MIN(r_ij*d_ij,a_ji+d_ij)*(u_i-u_j);
210    
211 gross 1361 fa[i]+=f_ij;
212     fa[j]-=f_ij;
213 gross 1371 printf("%d %d => %e %e : %e %e : %e %e : fa[%d]=%e fa[%d]=%e\n",i,j,d_ij,(u_i-u_j), a_ij, a_ji, r_ij,f_ij,i,fa[i],j,fa[j]);
214    
215    
216 gross 1361 }
217 gross 1371 break;
218 gross 1361 }
219     }
220     }
221     }
222     /* anti diffusive flux from couple block */
223    
224     /* TODO */
225     }
226     }
227     #pragma omp barrier
228     }
229     }
230     }

  ViewVC Help
Powered by ViewVC 1.1.26