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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1400 - (show annotations)
Thu Jan 24 06:04:31 2008 UTC (12 years, 1 month ago) by gross
File MIME type: text/plain
File size: 9647 byte(s)
better test example for upwinding added
1 /* $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 #define MC(a,b) (FLUX_S(a,b)*MIN3(ABS((a)+(b))/2,2*ABS(a),2*ABS(b)))
33
34 #define FLUX_L(a,b) MC(a,b) /* alter for other flux limiter */
35
36 #define FLUX_LIMITER(a) FLUX_L(1,a)
37
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 void Paso_FCTransportProblem_addAdvectivePart(Paso_FCTransportProblem * fc, double alpha) {
49 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 n=Paso_SystemMatrix_getTotalNumRows(fc->flux_matrix);
55
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 fc->transport_matrix->mainBlock->val[fc->main_iptr[i]]+=alpha*fc->flux_matrix->mainBlock->val[fc->main_iptr[i]];
64
65 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 if (j<i) {
68 /* find entry a[j,i] */
69 for (iptr_ji=fc->flux_matrix->mainBlock->pattern->ptr[j]; iptr_ji<fc->flux_matrix->mainBlock->pattern->ptr[j+1]; ++iptr_ji) {
70 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 /* printf("%d %d -> %e\n",i,j,d_ij); */
73 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 /* 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 fc->transport_matrix->mainBlock->val[fc->main_iptr[i]]-=d_ij;
77 fc->transport_matrix->mainBlock->val[fc->main_iptr[j]]-=d_ij;
78 break;
79 }
80 }
81 }
82
83 }
84 /* TODO process couple block */
85 }
86 }
87 #pragma omp barrier
88 }
89 }
90
91 }
92
93 void Paso_FCTransportProblem_setFlux(Paso_FCTransportProblem * fc, double * u, double* fa) {
94 /*
95 * sets fa=transport_matrix*u+anti-diffuison_flux(u)
96 */
97
98 double *remote_u=NULL;
99
100 if (fc==NULL) return;
101 Paso_SystemMatrix_startCollect(fc->transport_matrix,u);
102 /* process main block */
103 Paso_SparseMatrix_MatrixVector_CSR_OFFSET0(1.,fc->transport_matrix->mainBlock,u,0.,fa);
104 /* finish exchange */
105 remote_u=Paso_SystemMatrix_finishCollect(fc->transport_matrix);
106 /* process couple block */
107 Paso_SparseMatrix_MatrixVector_CSR_OFFSET0(1.,fc->transport_matrix->coupleBlock,remote_u,1.,fa);
108
109 Paso_FCTransportProblem_setAntiDiffusiveFlux(fc,u,remote_u,fa);
110 }
111 /**************************************************************/
112
113 /* adds antidiffusion to fa
114
115 P_p[i] + = sum_j min(0,a[i,j]) min(0,u[j]-u[i])
116 P_n[i] + = sum_j min(0,a[i,j]) max(0,u[j]-u[i])
117 Q_p[i] + = sum_j max(0,a[i,j]) max(0,u[j]-u[i])
118 Q_n[i] + = sum_j max(0,a[i,j]) min(0,u[j]-u[i])
119 d_ij=max(0,-a[i,j],-a[j,i])
120 l_ji=max(0,a[j,i],a[j,i]-a[i,j])
121 if a[j,i] >= a[i,j] and 0>a[i,j] : (i.e d_ij>0 and l_ji>=l_ij)
122 r_ij = u[i]>u[j] ? Q_p[i]/P_p[i] : Q_n[i]/Q_n[i]
123 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])
124 fa[i]+=f_ij
125 fa[j]-=f_ij
126
127 */
128
129 void Paso_FCTransportProblem_setAntiDiffusiveFlux(Paso_FCTransportProblem * fc, double * u, double *u_remote, double* fa) {
130
131 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;
132 index_t color, iptr_ij,j,iptr_ji, i;
133 dim_t n;
134
135
136 if (fc==NULL) return;
137 n=Paso_SystemMatrix_getTotalNumRows(fc->flux_matrix);
138 /* exchange */
139
140
141 #pragma omp parallel private(color)
142 {
143 for (color=0;color<fc->num_colors;++color) {
144 #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)
145 for (i = 0; i < n; ++i) {
146 if (fc->colorOf[i]==color) {
147 u_i=u[i];
148 /* gather the smoothness sensor */
149 P_p=0.;
150 P_n=0.;
151 Q_p=0.;
152 Q_n=0.;
153 r_p=0.;
154 r_n=0.;
155 /* #pragma ivdep */
156 for (iptr_ij=(fc->flux_matrix->mainBlock->pattern->ptr[i]);iptr_ij<(fc->flux_matrix->mainBlock->pattern->ptr[i+1]); ++iptr_ij) {
157 a_ij=fc->flux_matrix->mainBlock->val[iptr_ij];
158 j=fc->flux_matrix->mainBlock->pattern->index[iptr_ij];
159 d=u[j]-u_i;
160 /* printf("%d %d : %e %e :: %e \n",i,j,u_i,u[j],a_ij); */
161 if (a_ij<0.) {
162 if (d<0.) {
163 P_p+=a_ij*d;
164 } else {
165 P_n+=a_ij*d;
166 }
167 } else {
168 if (d>0.) {
169 Q_p+=a_ij*d;
170 } else {
171 Q_n+=a_ij*d;
172 }
173 }
174 }
175 #pragma ivdep
176 for (iptr_ij=(fc->flux_matrix->coupleBlock->pattern->ptr[i]);iptr_ij<(fc->flux_matrix->coupleBlock->pattern->ptr[i+1]); ++iptr_ij) {
177 a_ij=fc->flux_matrix->coupleBlock->val[iptr_ij];
178 j=fc->flux_matrix->coupleBlock->pattern->index[iptr_ij];
179 d=u_remote[j]-u_i;
180
181 if (a_ij<0.) {
182 if (d<0.) {
183 P_p+=a_ij*d;
184 } else {
185 P_n+=a_ij*d;
186 }
187 } else {
188 if (d>0.) {
189 Q_p+=a_ij*d;
190 } else {
191 Q_n+=a_ij*d;
192 }
193 }
194 }
195 /* set the smoothness indicators */
196 r_p = (P_p > 0.) ? FLUX_LIMITER(Q_p/P_p) : 0.;
197 r_n = (P_n < 0.) ? FLUX_LIMITER(Q_n/P_n) : 0.;
198 /* printf("Flux control %d: %e %e %e : %e %e %e : %e\n",i,Q_p,P_p,r_p,Q_n,P_n,r_n,u_i); */
199 /* anti diffusive flux from main block */
200 for (iptr_ij=fc->flux_matrix->mainBlock->pattern->ptr[i];iptr_ij<fc->flux_matrix->mainBlock->pattern->ptr[i+1]; ++iptr_ij) {
201 a_ij=fc->flux_matrix->mainBlock->val[iptr_ij];
202 j=fc->flux_matrix->mainBlock->pattern->index[iptr_ij];
203 if ( i!=j ) {
204 /* find entry a[j,i] */
205 for (iptr_ji=fc->flux_matrix->mainBlock->pattern->ptr[j];iptr_ji<fc->flux_matrix->mainBlock->pattern->ptr[j+1]; ++iptr_ji) {
206 if (fc->flux_matrix->mainBlock->pattern->index[iptr_ji]==i) {
207 a_ji=fc->flux_matrix->mainBlock->val[iptr_ji];
208 d_ij=-MIN3(0,a_ji,a_ij);
209 if ( (d_ij > 0.) && ( (a_ji > a_ij) || ( (a_ji == a_ij) && (j<i) ) ) ) {
210 u_j=u[j];
211 r_ij = u_i>u_j ? r_p : r_n;
212 f_ij =MIN(r_ij*d_ij,a_ji+d_ij)*(u_i-u_j);
213
214 fa[i]+=f_ij;
215 fa[j]-=f_ij;
216 /* 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]); */
217
218
219
220 }
221 break;
222 }
223 }
224 }
225 }
226 /* anti diffusive flux from couple block */
227
228 /* TODO */
229 }
230 }
231 #pragma omp barrier
232 }
233 }
234 }

  ViewVC Help
Powered by ViewVC 1.1.26