/[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 1375 - (show annotations)
Wed Jan 9 00:15:05 2008 UTC (11 years, 10 months ago) by gross
File MIME type: text/plain
File size: 9535 byte(s)
bug in interpolation at reduced face elements fixed.
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 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 /**************************************************************/
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 void Paso_FCTransportProblem_setAntiDiffusiveFlux(Paso_FCTransportProblem * fc, double * u, double *u_remote, double* fa) {
127
128 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 index_t color, iptr_ij,j,iptr_ji, i;
130 dim_t n;
131
132
133 if (fc==NULL) return;
134 n=Paso_SystemMatrix_getTotalNumRows(fc->flux_matrix);
135 /* exchange */
136
137
138 #pragma omp parallel private(color)
139 {
140 for (color=0;color<fc->num_colors;++color) {
141 #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 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 r_p=0.;
151 r_n=0.;
152 /* #pragma ivdep */
153 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 d=u[j]-u_i;
157 printf("%d %d : %e %e :: %e \n",i,j,u_i,u[j],a_ij);
158 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 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 d=u_remote[j]-u_i;
177
178 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 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 /* anti diffusive flux from main block */
197 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 if ( i!=j ) {
201 /* find entry a[j,i] */
202 for (iptr_ji=fc->flux_matrix->mainBlock->pattern->ptr[j];iptr_ji<fc->flux_matrix->mainBlock->pattern->ptr[j+1]; ++iptr_ji) {
203 if (fc->flux_matrix->mainBlock->pattern->index[iptr_ji]==i) {
204 a_ji=fc->flux_matrix->mainBlock->val[iptr_ji];
205 d_ij=-MIN3(0,a_ji,a_ij);
206 if ( (d_ij > 0.) && ( (a_ji > a_ij) || ( (a_ji == a_ij) && (j<i) ) ) ) {
207 u_j=u[j];
208 r_ij = u_i>u_j ? r_p : r_n;
209 f_ij =MIN(r_ij*d_ij,a_ji+d_ij)*(u_i-u_j);
210
211 fa[i]+=f_ij;
212 fa[j]-=f_ij;
213 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
217 }
218 break;
219 }
220 }
221 }
222 }
223 /* anti diffusive flux from couple block */
224
225 /* TODO */
226 }
227 }
228 #pragma omp barrier
229 }
230 }
231 }

  ViewVC Help
Powered by ViewVC 1.1.26