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

  ViewVC Help
Powered by ViewVC 1.1.26