| 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 |
break; |
| 218 |
} |
| 219 |
} |
| 220 |
} |
| 221 |
} |
| 222 |
/* anti diffusive flux from couple block */ |
| 223 |
|
| 224 |
/* TODO */ |
| 225 |
} |
| 226 |
} |
| 227 |
#pragma omp barrier |
| 228 |
} |
| 229 |
} |
| 230 |
} |
Parent Directory
| 
Revision Log