/[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 1361 - (show annotations)
Fri Dec 14 09:26:51 2007 UTC (12 years, 10 months ago) by gross
Original Path: trunk/paso/src/Solver_FluxControl.c
File MIME type: text/plain
File size: 10444 byte(s)
first steps towards a flux controlled transport solver.
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
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 /* free all memory used by FluxControl */
40
41 void Paso_Solver_FluxControl_free(Paso_Solver_FluxControl* in) {
42 if (in!=NULL) {
43 Paso_SystemMatrix_freeBuffer(in->matrix);
44 Paso_SystemMatrix_free(in->matrix);
45 MEMFREE(in->colorOf);
46 MEMFREE(in->main_iptr);
47 MEMFREE(in);
48 }
49 }
50
51 /**************************************************************/
52
53 /* constructs a flux control mechanism */
54
55 Paso_Solver_FluxControl* Paso_SolverFCT_getFluxControl(Paso_SystemMatrix * A) {
56
57 Paso_Solver_FluxControl* out=NULL;
58 dim_t n,i;
59 index_t iptr,iptr_main,k;
60
61 if (A==NULL) return out;
62 n=Paso_SystemMatrix_getTotalNumRows(A);
63 if (A->block_size!=1) {
64 Paso_setError(TYPE_ERROR,"Paso_SolverFCT_getFluxControl: block size > 1 is not supported.");
65 return NULL;
66 }
67 out=MEMALLOC(1,Paso_Solver_FluxControl);
68 if (Paso_checkPtr(out)) return NULL;
69
70 out->matrix=Paso_SystemMatrix_reference(A);
71 out->colorOf=NULL;
72 out->main_iptr=NULL;
73
74
75 /* allocations: */
76 out->colorOf=MEMALLOC(n,index_t);
77 out->main_iptr=MEMALLOC(n,index_t);
78 if ( ! (Paso_checkPtr(out->colorOf) || Paso_checkPtr(out->main_iptr) ) ) {
79 printf("Paso_SolverFCT_getFluxControl: Revise coloring!!\n");
80 Paso_Pattern_color(A->mainBlock->pattern,&(out->num_colors),out->colorOf);
81 Paso_SystemMatrix_allocBuffer(A);
82
83 #pragma omp parallel for schedule(static) private(i,iptr,iptr_main,k)
84 for (i = 0; i < n; ++i) {
85 for (iptr=A->mainBlock->pattern->ptr[i];iptr<A->mainBlock->pattern->ptr[i+1]; ++iptr) {
86 iptr_main=A->mainBlock->pattern->ptr[0]-1;
87 for (iptr=A->mainBlock->pattern->ptr[i];iptr<A->mainBlock->pattern->ptr[i+1]; iptr++) {
88 if (A->mainBlock->pattern->index[iptr]==i) {
89 iptr_main=iptr;
90 break;
91 }
92 }
93 out->main_iptr[i]=iptr_main;
94 if (iptr_main==A->mainBlock->pattern->ptr[0]-1)
95 Paso_setError(VALUE_ERROR, "Paso_SolverFCT_getFluxControl: no main diagonal");
96 }
97 }
98
99 }
100 if (Paso_noError()) {
101 return out;
102 } else {
103 Paso_Solver_FluxControl_free(out);
104 return NULL;
105 }
106 }
107
108 /**************************************************************/
109
110 /* adds A plus stabelising diffusion into the matrix B */
111
112 /* d_ij=alpha*max(0,-a[i,j],-a[j,i]) */
113 /* b[i,j]+=alpha*(a[i,j]+d_ij) */
114 /* b[j,i]+=alpha*(a[j,i]+d_ij) */
115 /* b[i,i]-=alpha*d_ij */
116 /* b[j,j]-=alpha*d_ij */
117
118 void Paso_Solver_FluxControl_addDiffusion(Paso_Solver_FluxControl * fc, double alpha, Paso_SystemMatrix * B) {
119 dim_t n,i;
120 index_t color, iptr_ij,j,iptr_ji;
121 register double d_ij;
122
123 if (fc==NULL) return;
124 n=Paso_SystemMatrix_getTotalNumRows(fc->matrix);
125 /* TODO test - same pattern + block size */
126
127 #pragma omp parallel private(color)
128 {
129 /* process main block */
130 for (color=0;color<fc->num_colors;++color) {
131 #pragma omp for private(i,iptr_ij,j,iptr_ji,d_ij) schedule(static)
132 for (i = 0; i < n; ++i) {
133 if (fc->colorOf[i]==color) {
134 for (iptr_ij=fc->matrix->mainBlock->pattern->ptr[i];iptr_ij<fc->matrix->mainBlock->pattern->ptr[i+1]; ++iptr_ij) {
135 j=fc->matrix->mainBlock->pattern->index[iptr_ij];
136 if (i<j) {
137 /* find entry a[j,i] */
138 for (iptr_ji=fc->matrix->mainBlock->pattern->ptr[i];iptr_ji<fc->matrix->mainBlock->pattern->ptr[j+1]-1; ++iptr_ji) {
139 if (fc->matrix->mainBlock->pattern->index[iptr_ji]==i) {
140 d_ij=(-alpha)*MIN3(0.,fc->matrix->mainBlock->val[iptr_ij],fc->matrix->mainBlock->val[iptr_ji]);
141 B->mainBlock->val[iptr_ij]+=alpha*fc->matrix->mainBlock->val[iptr_ij]+d_ij;
142 B->mainBlock->val[iptr_ji]+=alpha*fc->matrix->mainBlock->val[iptr_ji]+d_ij;
143 B->mainBlock->val[fc->main_iptr[i]]-=d_ij;
144 B->mainBlock->val[fc->main_iptr[j]]-=d_ij;
145 break;
146 }
147 }
148 }
149
150 }
151 /* TODO process couple block */
152 }
153 }
154 #pragma omp barrier
155 }
156 }
157
158 }
159 /**************************************************************/
160
161 /* adds antidiffusion to fa
162
163 P_p[i] + = sum_j min(0,a[i,j]) min(0,u[j]-u[i])
164 P_n[i] + = sum_j min(0,a[i,j]) max(0,u[j]-u[i])
165 Q_p[i] + = sum_j max(0,a[i,j]) max(0,u[j]-u[i])
166 Q_n[i] + = sum_j max(0,a[i,j]) min(0,u[j]-u[i])
167 d_ij=max(0,-a[i,j],-a[j,i])
168 l_ji=max(0,a[j,i],a[j,i]-a[i,j])
169 if a[j,i] >= a[i,j] and 0>a[i,j] : (i.e d_ij>0 and l_ji>=l_ij)
170 r_ij = u[i]>u[j] ? Q_p[i]/P_p[i] : Q_n[i]/Q_n[i]
171 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])
172 fa[i]+=f_ij
173 fa[j]-=f_ij
174
175 */
176
177 void Paso_Solver_FluxControl_setAntiDiffusiveFlux(Paso_Solver_FluxControl * fc, double * u, double* fa) {
178
179 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;
180 double *u_remote=NULL;
181 index_t color, iptr_ij,j,iptr_ji, i;
182 dim_t n;
183
184
185 if (fc==NULL) return;
186 n=Paso_SystemMatrix_getTotalNumRows(fc->matrix);
187 /* exchange */
188 Paso_SystemMatrix_startCollect(fc->matrix,u);
189 u_remote=Paso_SystemMatrix_finishCollect(fc->matrix);
190
191 #pragma omp parallel private(color)
192 {
193 for (color=0;color<fc->num_colors;++color) {
194 #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)
195 for (i = 0; i < n; ++i) {
196 if (fc->colorOf[i]==color) {
197 u_i=u[i];
198 /* gather the smoothness sensor */
199 P_p=0.;
200 P_n=0.;
201 Q_p=0.;
202 Q_n=0.;
203 #pragma ivdep
204 for (iptr_ij=(fc->matrix->mainBlock->pattern->ptr[i]);iptr_ij<(fc->matrix->mainBlock->pattern->ptr[i+1]); ++iptr_ij) {
205 a_ij=fc->matrix->mainBlock->val[iptr_ij];
206 j=fc->matrix->mainBlock->pattern->index[iptr_ij];
207 d=u[j]-u_i;
208 if (a_ij<0.) {
209 if (d<0.) {
210 P_p+=a_ij*d;
211 } else {
212 P_n+=a_ij*d;
213 }
214 } else {
215 if (d>0.) {
216 Q_p+=a_ij*d;
217 } else {
218 Q_n+=a_ij*d;
219 }
220 }
221 }
222 #pragma ivdep
223 for (iptr_ij=(fc->matrix->coupleBlock->pattern->ptr[i]);iptr_ij<(fc->matrix->coupleBlock->pattern->ptr[i+1]); ++iptr_ij) {
224 a_ij=fc->matrix->coupleBlock->val[iptr_ij];
225 j=fc->matrix->coupleBlock->pattern->index[iptr_ij];
226 d=u_remote[j]-u_i;
227 if (a_ij<0.) {
228 if (d<0.) {
229 P_p+=a_ij*d;
230 } else {
231 P_n+=a_ij*d;
232 }
233 } else {
234 if (d>0.) {
235 Q_p+=a_ij*d;
236 } else {
237 Q_n+=a_ij*d;
238 }
239 }
240 }
241 /* set the smoothness indicators */
242 r_p = (P_p > 0.) ? FLUX_LIMITER(Q_p/P_p) : 0.;
243 r_n = (P_n < 0 ) ? FLUX_LIMITER(Q_n/P_n) : 0.;
244 /* anti diffusive flux from main block */
245 for (iptr_ij=fc->matrix->mainBlock->pattern->ptr[i];iptr_ij<fc->matrix->mainBlock->pattern->ptr[i+1]; ++iptr_ij) {
246 a_ij=fc->matrix->mainBlock->val[iptr_ij];
247 j=fc->matrix->mainBlock->pattern->index[iptr_ij];
248 if (a_ij < 0 && i!=j) {
249 /* find entry a[j,i] */
250 for (iptr_ji=fc->matrix->mainBlock->pattern->ptr[i];iptr_ji<fc->matrix->mainBlock->pattern->ptr[j+1]-1; ++iptr_ji) {
251 if (fc->matrix->mainBlock->pattern->index[iptr_ji]==i) {
252 a_ji=fc->matrix->mainBlock->val[iptr_ji];
253 if (a_ji > a_ij || (a_ji == a_ij && j<i) ) {
254 u_j=u[j];
255 r_ij = u_i>u_j ? r_p : r_n;
256 f_ij =MIN(r_ij*a_ij,a_ji-a_ij)*(u_i-u_j);
257 fa[i]+=f_ij;
258 fa[j]-=f_ij;
259 break;
260 }
261 }
262 }
263 }
264 }
265 /* anti diffusive flux from couple block */
266
267 /* TODO */
268 }
269 }
270 #pragma omp barrier
271 }
272 }
273 }

  ViewVC Help
Powered by ViewVC 1.1.26