CG  Version 25
selfAdjointArtificialDiffusion.h
Go to the documentation of this file.
1  ! Face centered derivatives for the self-adjoint artificial diffusion
2  ! p=plus, m=minus, z=zero
3  ! Rectangular grid
4  uxmzzR(i1,i2,i3,c)=(u(i1,i2,i3,c)-u(i1-1,i2,i3,c))*dxi
5  uymzzR(i1,i2,i3,c)=(u(i1,i2+1,i3,c)-u(i1,i2-1,i3,c)+u(i1-1,i2+1,i3,c)-u(i1-1,i2-1,i3,c))*dyi*.25
6  uzmzzR(i1,i2,i3,c)=(u(i1,i2,i3+1,c)-u(i1,i2,i3-1,c)+u(i1-1,i2,i3+1,c)-u(i1-1,i2,i3-1,c))*dzi*.25
7 
8  uxzmzR(i1,i2,i3,c)=(u(i1+1,i2,i3,c)-u(i1-1,i2,i3,c)+u(i1+1,i2-1,i3,c)-u(i1-1,i2-1,i3,c))*dxi*.25
9  uyzmzR(i1,i2,i3,c)=(u(i1,i2,i3,c)-u(i1,i2-1,i3,c))*dyi
10  uzzmzR(i1,i2,i3,c)=(u(i1,i2,i3+1,c)-u(i1,i2,i3-1,c)+u(i1,i2-1,i3+1,c)-u(i1,i2-1,i3-1,c))*dzi*.25
11 
12  uxzzmR(i1,i2,i3,c)=(u(i1+1,i2,i3,c)-u(i1-1,i2,i3,c)+u(i1+1,i2,i3-1,c)-u(i1-1,i2,i3-1,c))*dxi*.25
13  uyzzmR(i1,i2,i3,c)=(u(i1,i2+1,i3,c)-u(i1,i2-1,i3,c)+u(i1,i2+1,i3-1,c)-u(i1,i2-1,i3-1,c))*dyi*.25
14  uzzzmR(i1,i2,i3,c)=(u(i1,i2,i3,c)-u(i1,i2,i3-1,c))*dzi
15 
16  ! curvilinear grid
17  udmzC(i1,i2,i3,m,c)=(rsxy(i1,i2,i3,0,m)+rsxy(i1-1,i2,i3,0,m))*(u(i1,i2,i3,c)-u(i1-1,i2 ,i3,c))*dr2i +
18  & (rsxy(i1,i2,i3,1,m)+rsxy(i1-1,i2,i3,1,m))*(
19  & u(i1,i2+1,i3,c)-u(i1,i2-1,i3,c)+ u(i1-1,i2+1,i3,c)-u(i1-1,i2-1,i3,c))*dsi*.125
20  udzmC(i1,i2,i3,m,c)=(rsxy(i1,i2,i3,1,m)+rsxy(i1,i2-1,i3,1,m))*(u(i1,i2,i3,c)-u(i1,i2-1,i3,c))*ds2i +
21  & (rsxy(i1,i2,i3,0,m)+rsxy(i1,i2-1,i3,0,m))*(
22  & u(i1+1,i2,i3,c)-u(i1-1,i2,i3,c)+ u(i1+1,i2-1,i3,c)-u(i1-1,i2-1,i3,c))*dri*.125
23 
24  udmzzC(i1,i2,i3,m,c)=(rsxy(i1,i2,i3,0,m)+rsxy(i1-1,i2,i3,0,m))*(u(i1,i2,i3,c)-u(i1-1,i2 ,i3,c))*dr2i +
25  & (rsxy(i1,i2,i3,1,m)+rsxy(i1-1,i2,i3,1,m))*(
26  & u(i1,i2+1,i3,c)-u(i1,i2-1,i3,c)+ u(i1-1,i2+1,i3,c)-u(i1-1,i2-1,i3,c))*dsi*.125+
27  & (rsxy(i1,i2,i3,2,m)+rsxy(i1-1,i2,i3,2,m))*(
28  & u(i1,i2,i3+1,c)-u(i1,i2,i3-1,c)+ u(i1-1,i2,i3+1,c)-u(i1-1,i2,i3-1,c))*dti*.125
29  udzmzC(i1,i2,i3,m,c)=(rsxy(i1,i2,i3,1,m)+rsxy(i1,i2-1,i3,1,m))*(u(i1,i2,i3,c)-u(i1,i2-1,i3,c))*ds2i +
30  & (rsxy(i1,i2,i3,0,m)+rsxy(i1,i2-1,i3,0,m))*(
31  & u(i1+1,i2,i3,c)-u(i1-1,i2,i3,c)+ u(i1+1,i2-1,i3,c)-u(i1-1,i2-1,i3,c))*dri*.125+
32  & (rsxy(i1,i2,i3,2,m)+rsxy(i1,i2-1,i3,2,m))*(
33  & u(i1,i2,i3+1,c)-u(i1,i2,i3-1,c)+ u(i1,i2-1,i3+1,c)-u(i1,i2-1,i3-1,c))*dti*.125
34 
35  udzzmC(i1,i2,i3,m,c)=(rsxy(i1,i2,i3,2,m)+rsxy(i1,i2,i3-1,2,m))*(u(i1,i2,i3,c)-u(i1,i2,i3-1,c))*dt2i +
36  & (rsxy(i1,i2,i3,0,m)+rsxy(i1,i2,i3-1,0,m))*(
37  & u(i1+1,i2,i3,c)-u(i1-1,i2,i3,c)+ u(i1+1,i2,i3-1,c)-u(i1-1,i2,i3-1,c))*dri*.125+
38  & (rsxy(i1,i2,i3,1,m)+rsxy(i1,i2,i3-1,1,m))*(
39  & u(i1,i2+1,i3,c)-u(i1,i2-1,i3,c)+ u(i1,i2+1,i3-1,c)-u(i1,i2-1,i3-1,c))*dsi*.125
40 
41  ! Coefficients of the artificial diffusion for the momentum equations
42  ! 2D - rectangular
43  admzR(i1,i2,i3)=ad21+cd22*( abs(uxmzzR(i1,i2,i3,uc))+abs(uxmzzR(i1,i2,i3,vc))+
44  & abs(uymzzR(i1,i2,i3,uc))+abs(uymzzR(i1,i2,i3,vc)) )
45 
46  adzmR(i1,i2,i3)=ad21+cd22*( abs(uxzmzR(i1,i2,i3,uc))+abs(uxzmzR(i1,i2,i3,vc))+
47  & abs(uyzmzR(i1,i2,i3,uc))+abs(uyzmzR(i1,i2,i3,vc)) )
48 
49  ! 3D
50  admzzR(i1,i2,i3)=ad21+cd22*( abs(uxmzzR(i1,i2,i3,uc))+abs(uxmzzR(i1,i2,i3,vc))+abs(uxmzzR(i1,i2,i3,wc))+
51  & abs(uymzzR(i1,i2,i3,uc))+abs(uymzzR(i1,i2,i3,vc))+abs(uymzzR(i1,i2,i3,wc))+
52  & abs(uzmzzR(i1,i2,i3,uc))+abs(uzmzzR(i1,i2,i3,vc))+abs(uzmzzR(i1,i2,i3,wc)) )
53 
54  adzmzR(i1,i2,i3)=ad21+cd22*( abs(uxzmzR(i1,i2,i3,uc))+abs(uxzmzR(i1,i2,i3,vc))+abs(uxzmzR(i1,i2,i3,wc))+
55  & abs(uyzmzR(i1,i2,i3,uc))+abs(uyzmzR(i1,i2,i3,vc))+abs(uyzmzR(i1,i2,i3,wc))+
56  & abs(uzzmzR(i1,i2,i3,uc))+abs(uzzmzR(i1,i2,i3,vc))+abs(uzzmzR(i1,i2,i3,wc)) )
57 
58  adzzmR(i1,i2,i3)=ad21+cd22*( abs(uxzzmR(i1,i2,i3,uc))+abs(uxzzmR(i1,i2,i3,vc))+abs(uxzzmR(i1,i2,i3,wc))+
59  & abs(uyzzmR(i1,i2,i3,uc))+abs(uyzzmR(i1,i2,i3,vc))+abs(uyzzmR(i1,i2,i3,wc))+
60  & abs(uzzzmR(i1,i2,i3,uc))+abs(uzzzmR(i1,i2,i3,vc))+abs(uzzzmR(i1,i2,i3,wc)) )
61  ! 2D - curvilinear
62  admzC(i1,i2,i3)=ad21+cd22*( abs(udmzC(i1,i2,i3,0,uc))+abs(udmzC(i1,i2,i3,0,vc))+
63  & abs(udmzC(i1,i2,i3,1,uc))+abs(udmzC(i1,i2,i3,1,vc)) )
64 
65  adzmC(i1,i2,i3)=ad21+cd22*( abs(udzmC(i1,i2,i3,0,uc))+abs(udzmC(i1,i2,i3,0,vc))+
66  & abs(udzmC(i1,i2,i3,1,uc))+abs(udzmC(i1,i2,i3,1,vc)) )
67 
68  ! 3D
69  admzzC(i1,i2,i3)=ad21+cd22*( abs(udmzzC(i1,i2,i3,0,uc))+abs(udmzzC(i1,i2,i3,0,vc))+abs(udmzzC(i1,i2,i3,0,wc))+
70  & abs(udmzzC(i1,i2,i3,1,uc))+abs(udmzzC(i1,i2,i3,1,vc))+abs(udmzzC(i1,i2,i3,1,wc))+
71  & abs(udmzzC(i1,i2,i3,2,uc))+abs(udmzzC(i1,i2,i3,2,vc))+abs(udmzzC(i1,i2,i3,2,wc)) )
72 
73  adzmzC(i1,i2,i3)=ad21+cd22*( abs(udzmzC(i1,i2,i3,0,uc))+abs(udzmzC(i1,i2,i3,0,vc))+abs(udzmzC(i1,i2,i3,0,wc))+
74  & abs(udzmzC(i1,i2,i3,1,uc))+abs(udzmzC(i1,i2,i3,1,vc))+abs(udzmzC(i1,i2,i3,1,wc))+
75  & abs(udzmzC(i1,i2,i3,2,uc))+abs(udzmzC(i1,i2,i3,2,vc))+abs(udzmzC(i1,i2,i3,2,wc)) )
76 
77  adzzmC(i1,i2,i3)=ad21+cd22*( abs(udzzmC(i1,i2,i3,0,uc))+abs(udzzmC(i1,i2,i3,0,vc))+abs(udzzmC(i1,i2,i3,0,wc))+
78  & abs(udzzmC(i1,i2,i3,1,uc))+abs(udzzmC(i1,i2,i3,1,vc))+abs(udzzmC(i1,i2,i3,1,wc))+
79  & abs(udzzmC(i1,i2,i3,2,uc))+abs(udzzmC(i1,i2,i3,2,vc))+abs(udzzmC(i1,i2,i3,2,wc)) )
80 
81  ! Coefficients of the artificial diffusion for the SA turbulence model
82  ! 2D - rectangular
83  admzRSA(i1,i2,i3)=ad21n+cd22n*( abs(uxmzzR(i1,i2,i3,nc))+abs(uymzzR(i1,i2,i3,nc)) )
84  adzmRSA(i1,i2,i3)=ad21n+cd22n*( abs(uxzmzR(i1,i2,i3,nc))+abs(uyzmzR(i1,i2,i3,nc)) )
85  ! 3D
86  admzzRSA(i1,i2,i3)=ad21n+cd22n*( abs(uxmzzR(i1,i2,i3,nc))+abs(uymzzR(i1,i2,i3,nc))+abs(uzmzzR(i1,i2,i3,nc)) )
87  adzmzRSA(i1,i2,i3)=ad21n+cd22n*( abs(uxzmzR(i1,i2,i3,nc))+abs(uyzmzR(i1,i2,i3,nc))+abs(uzzmzR(i1,i2,i3,nc)) )
88  adzzmRSA(i1,i2,i3)=ad21n+cd22n*( abs(uxzzmR(i1,i2,i3,nc))+abs(uyzzmR(i1,i2,i3,nc))+abs(uzzzmR(i1,i2,i3,nc)) )
89  ! 2D - curvilinear
90  admzCSA(i1,i2,i3)=ad21n+cd22n*( abs(udmzC(i1,i2,i3,0,nc))+abs(udmzC(i1,i2,i3,1,nc)) )
91  adzmCSA(i1,i2,i3)=ad21n+cd22n*( abs(udzmC(i1,i2,i3,0,nc))+abs(udzmC(i1,i2,i3,1,nc)) )
92  ! 3D
93  admzzCSA(i1,i2,i3)=ad21n+cd22n*( abs(udmzzC(i1,i2,i3,0,nc))+abs(udmzzC(i1,i2,i3,1,nc))+abs(udmzzC(i1,i2,i3,2,nc)))
94  adzmzCSA(i1,i2,i3)=ad21n+cd22n*( abs(udzmzC(i1,i2,i3,0,nc))+abs(udzmzC(i1,i2,i3,1,nc))+abs(udzmzC(i1,i2,i3,2,nc)))
95  adzzmCSA(i1,i2,i3)=ad21n+cd22n*( abs(udzzmC(i1,i2,i3,0,nc))+abs(udzzmC(i1,i2,i3,1,nc))+abs(udzzmC(i1,i2,i3,2,nc)))
96 
97 
98  ! Here are the parts of the artificial diffusion that are explicit (appear on the RHS)
99  adE0(i1,i2,i3,c) = cdzm*u(i1,i2-1,i3,c)+cdzp*u(i1,i2+1,i3,c)
100  adE1(i1,i2,i3,c) = cdmz*u(i1-1,i2,i3,c)+cdpz*u(i1+1,i2,i3,c)
101  adE2(i1,i2,i3,c) = 0.
102 
103  adE3d0(i1,i2,i3,c) = cdzmz*u(i1,i2-1,i3,c)+cdzpz*u(i1,i2+1,i3,c)+cdzzm*u(i1,i2,i3-1,c)+cdzzp*u(i1,i2,i3+1,c)
104  adE3d1(i1,i2,i3,c) = cdmzz*u(i1-1,i2,i3,c)+cdpzz*u(i1+1,i2,i3,c)+cdzzm*u(i1,i2,i3-1,c)+cdzzp*u(i1,i2,i3+1,c)
105  adE3d2(i1,i2,i3,c) = cdmzz*u(i1-1,i2,i3,c)+cdpzz*u(i1+1,i2,i3,c)+cdzmz*u(i1,i2-1,i3,c)+cdzpz*u(i1,i2+1,i3,c)
106 
107  ad2f(i1,i2,i3,m)= -cdDiag*u(i1,i2,i3,m)+cdmz*u(i1-1,i2,i3,m)+cdpz*u(i1+1,i2,i3,m)+
108  & cdzm*u(i1,i2-1,i3,m)+cdzp*u(i1,i2+1,i3,m)
109 
110  ad3f(i1,i2,i3,m)= -cdDiag*u(i1,i2,i3,m)+cdmzz*u(i1-1,i2,i3,m)+cdpzz*u(i1+1,i2,i3,m)+
111  & cdzmz*u(i1,i2-1,i3,m)+cdzpz*u(i1,i2+1,i3,m)+
112  & cdzzm*u(i1,i2,i3-1,m)+cdzzp*u(i1,i2,i3+1,m)
113 
114  ! Here are the full artificial diffusion terms
115  adSelfAdjoint2dR(i1,i2,i3,c)=admzR(i1 ,i2 ,i3 )*(u(i1-1,i2,i3,c)-u(i1,i2,i3,c))+
116  & admzR(i1+1,i2 ,i3 )*(u(i1+1,i2,i3,c)-u(i1,i2,i3,c))+
117  & adzmR(i1 ,i2 ,i3 )*(u(i1,i2-1,i3,c)-u(i1,i2,i3,c))+
118  & adzmR(i1 ,i2+1,i3 )*(u(i1,i2+1,i3,c)-u(i1,i2,i3,c))
119 
120  adSelfAdjoint3dR(i1,i2,i3,c)=admzzR(i1 ,i2 ,i3 )*(u(i1-1,i2,i3,c)-u(i1,i2,i3,c))+
121  & admzzR(i1+1,i2 ,i3 )*(u(i1+1,i2,i3,c)-u(i1,i2,i3,c))+
122  & adzmzR(i1 ,i2 ,i3 )*(u(i1,i2-1,i3,c)-u(i1,i2,i3,c))+
123  & adzmzR(i1 ,i2+1,i3 )*(u(i1,i2+1,i3,c)-u(i1,i2,i3,c))+
124  & adzzmR(i1 ,i2 ,i3 )*(u(i1,i2,i3-1,c)-u(i1,i2,i3,c))+
125  & adzzmR(i1 ,i2 ,i3+1)*(u(i1,i2,i3+1,c)-u(i1,i2,i3,c))
126 
127 
128  adSelfAdjoint2dC(i1,i2,i3,c)=admzC(i1 ,i2 ,i3 )*(u(i1-1,i2,i3,c)-u(i1,i2,i3,c))+
129  & admzC(i1+1,i2 ,i3 )*(u(i1+1,i2,i3,c)-u(i1,i2,i3,c))+
130  & adzmC(i1 ,i2 ,i3 )*(u(i1,i2-1,i3,c)-u(i1,i2,i3,c))+
131  & adzmC(i1 ,i2+1,i3 )*(u(i1,i2+1,i3,c)-u(i1,i2,i3,c))
132 
133  adSelfAdjoint3dC(i1,i2,i3,c)=admzzC(i1 ,i2 ,i3 )*(u(i1-1,i2,i3,c)-u(i1,i2,i3,c))+
134  & admzzC(i1+1,i2 ,i3 )*(u(i1+1,i2,i3,c)-u(i1,i2,i3,c))+
135  & adzmzC(i1 ,i2 ,i3 )*(u(i1,i2-1,i3,c)-u(i1,i2,i3,c))+
136  & adzmzC(i1 ,i2+1,i3 )*(u(i1,i2+1,i3,c)-u(i1,i2,i3,c))+
137  & adzzmC(i1 ,i2 ,i3 )*(u(i1,i2,i3-1,c)-u(i1,i2,i3,c))+
138  & adzzmC(i1 ,i2 ,i3+1)*(u(i1,i2,i3+1,c)-u(i1,i2,i3,c))
139 
140  ! Here are versions for the turbulence model
141  adSelfAdjoint2dRSA(i1,i2,i3,c)=admzRSA(i1 ,i2 ,i3 )*(u(i1-1,i2,i3,c)-u(i1,i2,i3,c))+
142  & admzRSA(i1+1,i2 ,i3 )*(u(i1+1,i2,i3,c)-u(i1,i2,i3,c))+
143  & adzmRSA(i1 ,i2 ,i3 )*(u(i1,i2-1,i3,c)-u(i1,i2,i3,c))+
144  & adzmRSA(i1 ,i2+1,i3 )*(u(i1,i2+1,i3,c)-u(i1,i2,i3,c))
145 
146  adSelfAdjoint3dRSA(i1,i2,i3,c)=admzzRSA(i1 ,i2 ,i3 )*(u(i1-1,i2,i3,c)-u(i1,i2,i3,c))+
147  & admzzRSA(i1+1,i2 ,i3 )*(u(i1+1,i2,i3,c)-u(i1,i2,i3,c))+
148  & adzmzRSA(i1 ,i2 ,i3 )*(u(i1,i2-1,i3,c)-u(i1,i2,i3,c))+
149  & adzmzRSA(i1 ,i2+1,i3 )*(u(i1,i2+1,i3,c)-u(i1,i2,i3,c))+
150  & adzzmRSA(i1 ,i2 ,i3 )*(u(i1,i2,i3-1,c)-u(i1,i2,i3,c))+
151  & adzzmRSA(i1 ,i2 ,i3+1)*(u(i1,i2,i3+1,c)-u(i1,i2,i3,c))
152 
153 
154  adSelfAdjoint2dCSA(i1,i2,i3,c)=admzCSA(i1 ,i2 ,i3 )*(u(i1-1,i2,i3,c)-u(i1,i2,i3,c))+
155  & admzCSA(i1+1,i2 ,i3 )*(u(i1+1,i2,i3,c)-u(i1,i2,i3,c))+
156  & adzmCSA(i1 ,i2 ,i3 )*(u(i1,i2-1,i3,c)-u(i1,i2,i3,c))+
157  & adzmCSA(i1 ,i2+1,i3 )*(u(i1,i2+1,i3,c)-u(i1,i2,i3,c))
158 
159  adSelfAdjoint3dCSA(i1,i2,i3,c)=admzzCSA(i1 ,i2 ,i3 )*(u(i1-1,i2,i3,c)-u(i1,i2,i3,c))+
160  & admzzCSA(i1+1,i2 ,i3 )*(u(i1+1,i2,i3,c)-u(i1,i2,i3,c))+
161  & adzmzCSA(i1 ,i2 ,i3 )*(u(i1,i2-1,i3,c)-u(i1,i2,i3,c))+
162  & adzmzCSA(i1 ,i2+1,i3 )*(u(i1,i2+1,i3,c)-u(i1,i2,i3,c))+
163  & adzzmCSA(i1 ,i2 ,i3 )*(u(i1,i2,i3-1,c)-u(i1,i2,i3,c))+
164  & adzzmCSA(i1 ,i2 ,i3+1)*(u(i1,i2,i3+1,c)-u(i1,i2,i3,c))
165 
166 
Generated on Fri Jan 4 2013 10:19:27 for CG by   doxygen 1.8.3