1 static char help[] = "Variable-Viscosity Stokes Problem in 2d.\n\
2 Exact solutions provided by Mirko Velic.\n\n\n";
3
4 #include<petsc.h>
5
6 #include "ex75.h"
7
8 typedef struct {
9 PetscBool fem; /* Flag for FEM tests */
10 } AppCtx;
11
ProcessOptions(MPI_Comm comm,AppCtx * options)12 PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *options)
13 {
14 PetscErrorCode ierr;
15
16 PetscFunctionBeginUser;
17 options->fem = PETSC_FALSE;
18
19 ierr = PetscOptionsBegin(comm, "", "Stokes Problem Options", "DMPLEX");CHKERRQ(ierr);
20 ierr = PetscOptionsBool("-fem", "Run FEM tests", "ex75.c", options->fem, &options->fem, NULL);CHKERRQ(ierr);
21 ierr = PetscOptionsEnd();
22 PetscFunctionReturn(0);
23 }
24
25 /*
26 SolKxSolution - Exact Stokes solutions for exponentially varying viscosity
27
28 Input Parameters:
29 + x - The x coordinate at which to evaluate the solution
30 . z - The z coordinate at which to evaluate the solution
31 . kn - The constant defining the x-dependence of the forcing function
32 . km - The constant defining the z-dependence of the forcing function
33 - B - The viscosity coefficient
34
35 Output Parameters:
36 + vx - The x-velocity at (x,z)
37 . vz - The z-velocity at (x,z)
38 . p - The pressure at (x,z)
39 . sxx - The stress sigma_xx at (x,z)
40 . sxz - The stress sigma_xz at (x,z)
41 - szz - The stress sigma_zz at (x,z)
42
43 Note:
44 $ The domain is the square 0 <= x,z <= 1. We solve the Stokes equation for incompressible flow with free-slip boundary
45 $ conditions everywhere. The forcing term f is given by
46 $
47 $ fx = 0
48 $ fz = sigma*sin(km*z)*cos(kn*x)
49 $
50 $ where
51 $
52 $ km = m*Pi (m may be non-integral)
53 $ kn = n*Pi
54 $
55 $ meaning that the density rho is -sigma*sin(km*z)*cos(kn*x). The viscosity eta is exp(2*B*x).
56 */
SolKxSolution(PetscReal x,PetscReal z,PetscReal kn,PetscReal km,PetscReal B,PetscScalar * vx,PetscScalar * vz,PetscScalar * p,PetscScalar * sxx,PetscScalar * sxz,PetscScalar * szz)57 PetscErrorCode SolKxSolution(PetscReal x, PetscReal z, PetscReal kn, PetscReal km, PetscReal B, PetscScalar *vx, PetscScalar *vz, PetscScalar *p, PetscScalar *sxx, PetscScalar *sxz, PetscScalar *szz)
58 {
59 PetscScalar sigma;
60 PetscScalar _C1,_C2,_C3,_C4;
61 PetscScalar Rp, UU, VV;
62 PetscScalar a,b,r,_aa,_bb,AA,BB,Rm;
63 PetscScalar num1,num2,num3,num4,den1;
64
65 PetscScalar t1,t2,t3,t4,t5,t6,t7,t8,t9,t10;
66 PetscScalar t11,t12,t13,t14,t15,t16,t17,t18,t19,t20,t21;
67 PetscScalar t22,t23,t24,t25,t26,t28,t29,t30,t31,t32;
68 PetscScalar t33,t34,t35,t36,t37,t38,t39,t40,t41,t42;
69 PetscScalar t44,t45,t46,t47,t48,t49,t51,t52,t53,t54;
70 PetscScalar t56,t58,t61,t62,t63,t64,t65,t66,t67,t68;
71 PetscScalar t69,t70,t71,t72,t73,t74,t75,t76,t77,t78;
72 PetscScalar t79,t80,t81,t82,t83,t84,t85,t86,t87,t88;
73 PetscScalar t89,t90,t91,t92,t93,t94,t95,t96,t97,t98;
74 PetscScalar t99,t100,t101,t103,t104,t105,t106,t107,t108,t109;
75 PetscScalar t110,t111,t112,t113,t114,t115,t116,t117,t118,t119;
76 PetscScalar t120,t121,t123,t125,t127,t128,t130,t131,t132,t133;
77 PetscScalar t135,t136,t138,t140,t141,t142,t143,t152,t160,t162;
78
79 PetscFunctionBegin;
80 /*************************************************************************/
81 /*************************************************************************/
82 /* rho = -sin(km*z)*cos(kn*x) */
83 /* viscosity Z= exp(2*B*z) */
84 /* solution valid for km not zero -- should get trivial solution if km=0 */
85 sigma = 1.0;
86 /*************************************************************************/
87 /*************************************************************************/
88 a = B*B + km*km;
89 b = 2.0*km*B;
90 r = sqrt(a*a + b*b);
91 Rp = sqrt( (r+a)/2.0);
92 Rm = sqrt( (r-a)/2.0);
93 UU = Rp - B;
94 VV = Rp + B;
95
96 /*******************************************/
97 /* calculate the constants */
98 /*******************************************/
99 t1 = kn * kn;
100 t4 = km * km;
101 t6 = t4 * t4;
102 t7 = B * B;
103 t9 = 0.4e1 * t7 * t4;
104 t12 = 0.8e1 * t7 * kn * km;
105 t14 = 0.4e1 * t7 * t1;
106 t16 = 0.2e1 * t4 * t1;
107 t17 = t1 * t1;
108 _aa = -0.4e1 * B * t1 * sigma * (t4 + t1) / (t6 + t9 + t12 + t14 + t16 + t17) / (t6 + t9 - t12 + t14 + t16 + t17);
109
110 t2 = kn * kn;
111 t3 = t2 * t2;
112 t4 = B * B;
113 t6 = 0.4e1 * t4 * t2;
114 t7 = km * km;
115 t9 = 0.4e1 * t7 * t4;
116 t10 = t7 * t7;
117 t12 = 0.2e1 * t7 * t2;
118 t16 = 0.8e1 * t4 * kn * km;
119 _bb = sigma * kn * (t3 - t6 + t9 + t10 + t12) / (t10 + t9 + t16 + t6 + t12 + t3) / (t10 + t9 - t16 + t6 + t12 + t3);
120
121 AA = _aa;
122 BB = _bb;
123
124 t1 = Rm * Rm;
125 t2 = B - Rp;
126 t4 = Rp + B;
127 t6 = UU * x;
128 t9 = exp(t6 - 0.4e1 * Rp);
129 t13 = kn * kn;
130 t15 = B * B;
131 t18 = Rp * Rp;
132 t19 = t18 * B;
133 t20 = t15 * Rp;
134 t22 = t1 * Rp;
135 t24 = B * t1;
136 t32 = 0.8e1 * t15 * BB * kn * Rp;
137 t34 = 0.2e1 * Rm;
138 t35 = cos(t34);
139 t37 = Rm * Rp;
140 t49 = sin(t34);
141 t63 = exp(t6 - 0.2e1 * Rp);
142 t65 = Rm * t2;
143 t67 = 0.2e1 * B * kn;
144 t68 = B * Rm;
145 t69 = t67 + t68 + t37;
146 t73 = 0.3e1 * t15;
147 t75 = 0.2e1 * B * Rp;
148 t76 = t73 - t75 + t1 - t13 - t18;
149 t78 = t65 * t76 * BB;
150 t80 = Rm - kn;
151 t81 = cos(t80);
152 t83 = t68 - t67 + t37;
153 t88 = Rm + kn;
154 t89 = cos(t88);
155 t92 = t65 * t76 * AA;
156 t97 = sin(t80);
157 t103 = sin(t88);
158 t108 = exp(t6 - 0.3e1 * Rp - B);
159 t110 = Rm * t4;
160 t111 = t67 + t68 - t37;
161 t115 = t73 + t75 + t1 - t13 - t18;
162 t117 = t110 * t115 * BB;
163 t120 = -t67 + t68 - t37;
164 t127 = t110 * t115 * AA;
165 t140 = exp(t6 - Rp - B);
166 num1 = -0.4e1 * t1 * t2 * t4 * AA * t9 + ((0.2e1 * Rp * (-B * t13 + 0.3e1 * t15 * B - t19 - 0.2e1 * t20 - 0.2e1 * t22 - t24) * AA - t32) * t35 + (0.2e1 * t37 * (t1 - t13 + 0.5e1 * t15 - t18) * AA - 0.8e1 * B * BB * kn * Rm * Rp) * t49 - 0.2e1 * B * (0.3e1 * t20 - t18 * Rp - 0.2e1 * t19 - Rp * t13 - t22 - 0.2e1 * t24) * AA + t32) * t63 + ((0.2e1 * t65 * t69 * AA + t78) * t81 + (0.2e1 * t65 * t83 * AA - t78) * t89 + (t92 - 0.2e1 * t65 * t69 * BB) * t97 + (t92 + 0.2e1 * t65 * t83 * BB) * t103) * t108 + ((-0.2e1 * t110 * t111 * AA - t117) * t81 + (-0.2e1 * t110 * t120 * AA + t117) * t89 + (-t127 + 0.2e1 * t110 * t111 * BB) * t97 + (-t127 - 0.2e1 * t110 * t120 * BB) * t103) * t140;
167
168 t1 = Rp + B;
169 t2 = Rm * t1;
170 t3 = B * B;
171 t4 = 0.3e1 * t3;
172 t5 = B * Rp;
173 t7 = Rm * Rm;
174 t8 = kn * kn;
175 t9 = Rp * Rp;
176 t10 = t4 + 0.2e1 * t5 + t7 - t8 - t9;
177 t12 = t2 * t10 * AA;
178 t14 = B * Rm;
179 t20 = UU * x;
180 t23 = exp(t20 - 0.4e1 * Rp);
181 t25 = Rm * Rp;
182 t32 = Rm * kn;
183 t37 = 0.2e1 * Rm;
184 t38 = cos(t37);
185 t41 = t3 * B;
186 t44 = t3 * Rp;
187 t48 = B * t7;
188 t53 = t3 * BB;
189 t54 = kn * Rp;
190 t58 = sin(t37);
191 t69 = exp(t20 - 0.2e1 * Rp);
192 t71 = t9 * Rp;
193 t72 = Rm * t71;
194 t73 = t3 * Rm;
195 t75 = 0.5e1 * t73 * Rp;
196 t77 = 0.8e1 * t44 * kn;
197 t78 = t25 * t8;
198 t79 = t7 * Rm;
199 t80 = B * t79;
200 t81 = t14 * t8;
201 t82 = t79 * Rp;
202 t84 = 0.3e1 * t41 * Rm;
203 t85 = t14 * t9;
204 t86 = -t72 + t75 + t77 - t78 + t80 - t81 + t82 + t84 + t85;
205 t88 = t7 * t9;
206 t89 = t5 * t8;
207 t90 = t7 * t3;
208 t91 = B * t71;
209 t92 = t48 * Rp;
210 t94 = 0.2e1 * t14 * t54;
211 t96 = 0.3e1 * Rp * t41;
212 t98 = 0.2e1 * t73 * kn;
213 t100 = 0.2e1 * t9 * t3;
214 t101 = -t88 - t89 - t90 - t91 - t92 - t94 + t96 - t98 - t100;
215 t105 = Rm - kn;
216 t106 = cos(t105);
217 t108 = t75 - t77 - t78 + t85 - t72 - t81 + t80 + t84 + t82;
218 t110 = -t100 + t96 - t91 + t94 + t98 - t92 - t89 - t88 - t90;
219 t114 = Rm + kn;
220 t115 = cos(t114);
221 t121 = sin(t105);
222 t127 = sin(t114);
223 t132 = exp(t20 - 0.3e1 * Rp - B);
224 t135 = 0.2e1 * B * kn;
225 t136 = t135 + t14 - t25;
226 t142 = -t135 + t14 - t25;
227 t152 = t2 * t10 * BB;
228 t162 = exp(t20 - Rp - B);
229 num2 = (0.2e1 * t12 - 0.8e1 * t14 * kn * t1 * BB) * t23 + ((-0.2e1 * t25 * (t7 - t8 + 0.5e1 * t3 - t9) * AA + 0.8e1 * B * BB * t32 * Rp) * t38 + (0.2e1 * Rp * (-B * t8 + 0.3e1 * t41 - t9 * B - 0.2e1 * t44 - 0.2e1 * t7 * Rp - t48) * AA - 0.8e1 * t53 * t54) * t58 - 0.2e1 * t14 * (t4 + t9 - t8 + t7) * AA + 0.8e1 * t53 * t32) * t69 + ((-t86 * AA - 0.2e1 * t101 * BB) * t106 + (-t108 * AA + 0.2e1 * t110 * BB) * t115 + (-0.2e1 * t101 * AA + t86 * BB) * t121 + (-0.2e1 * t110 * AA - t108 * BB) * t127) * t132 + ((t12 - 0.2e1 * t2 * t136 * BB) * t106 + (t12 + 0.2e1 * t2 * t142 * BB) * t115 + (-0.2e1 * t2 * t136 * AA - t152) * t121 + (-0.2e1 * t2 * t142 * AA + t152) * t127) * t162;
230
231 t1 = Rm * Rm;
232 t2 = B - Rp;
233 t4 = Rp + B;
234 t6 = VV * x;
235 t7 = exp(-t6);
236 t11 = B * t1;
237 t12 = Rp * Rp;
238 t13 = t12 * B;
239 t14 = B * B;
240 t15 = t14 * Rp;
241 t19 = kn * kn;
242 t21 = t1 * Rp;
243 t30 = 0.8e1 * t14 * BB * kn * Rp;
244 t32 = 0.2e1 * Rm;
245 t33 = cos(t32);
246 t35 = Rm * Rp;
247 t47 = sin(t32);
248 t61 = exp(-t6 - 0.2e1 * Rp);
249 t63 = Rm * t2;
250 t65 = 0.2e1 * B * kn;
251 t66 = B * Rm;
252 t67 = t65 + t66 + t35;
253 t71 = 0.3e1 * t14;
254 t73 = 0.2e1 * B * Rp;
255 t74 = t71 - t73 + t1 - t19 - t12;
256 t76 = t63 * t74 * BB;
257 t78 = Rm - kn;
258 t79 = cos(t78);
259 t81 = t66 - t65 + t35;
260 t86 = Rm + kn;
261 t87 = cos(t86);
262 t90 = t63 * t74 * AA;
263 t95 = sin(t78);
264 t101 = sin(t86);
265 t106 = exp(-t6 - 0.3e1 * Rp - B);
266 t108 = Rm * t4;
267 t109 = t65 + t66 - t35;
268 t113 = t71 + t73 + t1 - t19 - t12;
269 t115 = t108 * t113 * BB;
270 t118 = -t65 + t66 - t35;
271 t125 = t108 * t113 * AA;
272 t138 = exp(-t6 - Rp - B);
273 num3 = -0.4e1 * t1 * t2 * t4 * AA * t7 + ((-0.2e1 * Rp * (-t11 - t13 + 0.2e1 * t15 + 0.3e1 * t14 * B - B * t19 + 0.2e1 * t21) * AA + t30) * t33 + (-0.2e1 * t35 * (t1 - t19 + 0.5e1 * t14 - t12) * AA + 0.8e1 * B * BB * kn * Rm * Rp) * t47 + 0.2e1 * B * (-t12 * Rp + 0.2e1 * t11 + 0.3e1 * t15 + 0.2e1 * t13 - t21 - Rp * t19) * AA - t30) * t61 + ((-0.2e1 * t63 * t67 * AA - t76) * t79 + (-0.2e1 * t63 * t81 * AA + t76) * t87 + (-t90 + 0.2e1 * t63 * t67 * BB) * t95 + (-t90 - 0.2e1 * t63 * t81 * BB) * t101) * t106 + ((0.2e1 * t108 * t109 * AA + t115) * t79 + (0.2e1 * t108 * t118 * AA - t115) * t87 + (t125 - 0.2e1 * t108 * t109 * BB) * t95 + (t125 + 0.2e1 * t108 * t118 * BB) * t101) * t138;
274
275 t1 = B - Rp;
276 t2 = Rm * t1;
277 t3 = B * B;
278 t4 = 0.3e1 * t3;
279 t5 = B * Rp;
280 t7 = Rm * Rm;
281 t8 = kn * kn;
282 t9 = Rp * Rp;
283 t10 = t4 - 0.2e1 * t5 + t7 - t8 - t9;
284 t12 = t2 * t10 * AA;
285 t14 = B * Rm;
286 t20 = VV * x;
287 t21 = exp(-t20);
288 t23 = Rm * Rp;
289 t30 = Rm * kn;
290 t35 = 0.2e1 * Rm;
291 t36 = cos(t35);
292 t38 = B * t7;
293 t40 = t3 * Rp;
294 t42 = t3 * B;
295 t51 = t3 * BB;
296 t52 = kn * Rp;
297 t56 = sin(t35);
298 t67 = exp(-t20 - 0.2e1 * Rp);
299 t70 = 0.2e1 * B * kn;
300 t71 = t70 + t14 + t23;
301 t76 = Rm - kn;
302 t77 = cos(t76);
303 t79 = t14 - t70 + t23;
304 t84 = Rm + kn;
305 t85 = cos(t84);
306 t91 = t2 * t10 * BB;
307 t93 = sin(t76);
308 t99 = sin(t84);
309 t104 = exp(-t20 - 0.3e1 * Rp - B);
310 t106 = t9 * Rp;
311 t107 = Rm * t106;
312 t108 = t3 * Rm;
313 t110 = 0.5e1 * t108 * Rp;
314 t112 = 0.8e1 * t40 * kn;
315 t113 = t23 * t8;
316 t114 = t7 * Rm;
317 t115 = B * t114;
318 t116 = t14 * t8;
319 t117 = t114 * Rp;
320 t119 = 0.3e1 * t42 * Rm;
321 t120 = t14 * t9;
322 t121 = t107 - t110 - t112 + t113 + t115 - t116 - t117 + t119 + t120;
323 t123 = t38 * Rp;
324 t125 = 0.2e1 * t14 * t52;
325 t127 = 0.3e1 * Rp * t42;
326 t128 = t7 * t3;
327 t130 = 0.2e1 * t9 * t3;
328 t131 = t7 * t9;
329 t132 = B * t106;
330 t133 = t5 * t8;
331 t135 = 0.2e1 * t108 * kn;
332 t136 = -t123 - t125 + t127 + t128 + t130 + t131 - t132 - t133 + t135;
333 t141 = -t110 + t112 + t113 + t120 + t107 - t116 + t115 + t119 - t117;
334 t143 = t125 - t132 + t130 - t135 + t127 + t131 - t123 + t128 - t133;
335 t160 = exp(-t20 - Rp - B);
336 num4 = (0.2e1 * t12 - 0.8e1 * t14 * kn * t1 * BB) * t21 + ((0.2e1 * t23 * (t7 - t8 + 0.5e1 * t3 - t9) * AA - 0.8e1 * B * BB * t30 * Rp) * t36 + (-0.2e1 * Rp * (-t38 - t9 * B + 0.2e1 * t40 + 0.3e1 * t42 - B * t8 + 0.2e1 * t7 * Rp) * AA + 0.8e1 * t51 * t52) * t56 - 0.2e1 * t14 * (t4 + t9 - t8 + t7) * AA + 0.8e1 * t51 * t30) * t67 + ((t12 - 0.2e1 * t2 * t71 * BB) * t77 + (t12 + 0.2e1 * t2 * t79 * BB) * t85 + (-0.2e1 * t2 * t71 * AA - t91) * t93 + (-0.2e1 * t2 * t79 * AA + t91) * t99) * t104 + ((-t121 * AA + 0.2e1 * t136 * BB) * t77 + (-t141 * AA - 0.2e1 * t143 * BB) * t85 + (0.2e1 * t136 * AA + t121 * BB) * t93 + (0.2e1 * t143 * AA - t141 * BB) * t99) * t160;
337
338
339 t1 = Rm * Rm;
340 t2 = Rp * Rp;
341 t3 = t1 * t2;
342 t4 = B * B;
343 t5 = t1 * t4;
344 t9 = exp(-0.4e1 * Rp);
345 t15 = cos(0.2e1 * Rm);
346 t22 = exp(-0.2e1 * Rp);
347 den1 = (-0.4e1 * t3 + 0.4e1 * t5) * t9 + ((0.8e1 * t1 + 0.8e1 * t4) * t2 * t15 - 0.8e1 * t5 - 0.8e1 * t2 * t4) * t22 - 0.4e1 * t3 + 0.4e1 * t5;
348
349 _C1=num1/den1; _C2=num2/den1; _C3=num3/den1; _C4=num4/den1;
350
351 /*******************************************/
352 /* calculate solution */
353 /*******************************************/
354 t1 = Rm * x;
355 t2 = cos(t1);
356 t4 = sin(t1);
357 t10 = exp(-0.2e1 * x * B);
358 t12 = kn * x;
359 t13 = cos(t12);
360 t16 = sin(t12);
361 *vx = -km * (_C1 * t2 + _C2 * t4 + _C3 * t2 + _C4 * t4 + t10 * AA * t13 + t10 * BB * t16);
362
363 t2 = Rm * x;
364 t3 = cos(t2);
365 t6 = sin(t2);
366 t22 = exp(-0.2e1 * x * B);
367 t23 = B * t22;
368 t24 = kn * x;
369 t25 = cos(t24);
370 t29 = sin(t24);
371 *vz = UU * _C1 * t3 + UU * _C2 * t6 - _C1 * t6 * Rm + _C2 * t3 * Rm - VV * _C3 * t3 - VV * _C4 * t6 - _C3 * t6 * Rm + _C4 * t3 * Rm - 0.2e1 * t23 * AA * t25 - 0.2e1 * t23 * BB * t29 - t22 * AA * t29 * kn + t22 * BB * t25 * kn;
372
373 t3 = exp(0.2e1 * x * B);
374 t4 = t3 * B;
375 t8 = km * km;
376 t9 = t3 * t8;
377 t11 = 0.3e1 * t9 * Rm;
378 t12 = Rm * Rm;
379 t14 = t3 * t12 * Rm;
380 t15 = UU * UU;
381 t19 = 0.4e1 * t4 * UU * Rm - t11 - t14 + 0.3e1 * t3 * t15 * Rm;
382 t20 = Rm * x;
383 t21 = sin(t20);
384 t26 = 0.2e1 * t9 * B;
385 t33 = 0.2e1 * t4 * t12;
386 t36 = -t3 * t15 * UU - t26 + 0.3e1 * t9 * UU + 0.3e1 * t3 * UU * t12 + t33 - 0.2e1 * t4 * t15;
387 t37 = cos(t20);
388 t46 = VV * VV;
389 t53 = -t11 - t14 + 0.3e1 * t3 * t46 * Rm - 0.4e1 * t4 * VV * Rm;
390 t64 = -t26 + t33 + t3 * t46 * VV - 0.3e1 * t9 * VV - 0.2e1 * t4 * t46 - 0.3e1 * t3 * VV * t12;
391 t73 = kn * kn;
392 t74 = t73 * kn;
393 t79 = B * B;
394 t86 = B * t8;
395 t90 = kn * x;
396 t91 = sin(t90);
397 t106 = cos(t90);
398 *sxx = -((t19 * t21 + t36 * t37) * _C1 + (t36 * t21 - t19 * t37) * _C2 + (t53 * t21 + t64 * t37) * _C3 + (t64 * t21 - t53 * t37) * _C4 + (-AA * t74 - 0.4e1 * BB * t73 * B + 0.4e1 * t79 * AA * kn - 0.3e1 * t8 * AA * kn - 0.8e1 * t86 * BB) * t91 + (-0.8e1 * t86 * AA - 0.4e1 * AA * t73 * B - 0.4e1 * t79 * BB * kn + 0.3e1 * t8 * BB * kn + BB * t74) * t106) / km;
399
400 t3 = exp(0.2e1 * x * B);
401 t4 = km * km;
402 t5 = t3 * t4;
403 t6 = Rm * x;
404 t7 = cos(t6);
405 t8 = _C1 * t7;
406 t10 = sin(t6);
407 t11 = _C2 * t10;
408 t13 = _C3 * t7;
409 t15 = _C4 * t10;
410 t18 = kn * x;
411 t19 = cos(t18);
412 t22 = sin(t18);
413 t24 = UU * UU;
414 t25 = t3 * t24;
415 t28 = t3 * UU;
416 t38 = Rm * Rm;
417 t39 = t7 * t38;
418 t42 = t10 * t38;
419 t44 = t5 * t8 + t5 * t11 + t5 * t13 + t5 * t15 + t4 * AA * t19 + t4 * BB * t22 + t25 * t8 + t25 * t11 - 0.2e1 * t28 * _C1 * t10 * Rm + 0.2e1 * t28 * _C2 * t7 * Rm - t3 * _C1 * t39 - t3 * _C2 * t42;
420 t45 = VV * VV;
421 t46 = t3 * t45;
422 t49 = t3 * VV;
423 t62 = B * B;
424 t78 = kn * kn;
425 t82 = t46 * t13 + t46 * t15 + 0.2e1 * t49 * _C3 * t10 * Rm - 0.2e1 * t49 * _C4 * t7 * Rm - t3 * _C3 * t39 - t3 * _C4 * t42 + 0.4e1 * t62 * AA * t19 + 0.4e1 * t62 * BB * t22 + 0.4e1 * B * AA * t22 * kn - 0.4e1 * B * BB * t19 * kn - AA * t19 * t78 - BB * t22 * t78;
426 *sxz = t44 + t82;
427
428 t3 = exp(0.2e1 * x * B);
429 t4 = t3 * B;
430 t8 = km * km;
431 t9 = t3 * t8;
432 t10 = t9 * Rm;
433 t11 = Rm * Rm;
434 t13 = t3 * t11 * Rm;
435 t14 = UU * UU;
436 t18 = 0.4e1 * t4 * UU * Rm - t10 - t13 + 0.3e1 * t3 * t14 * Rm;
437 t19 = Rm * x;
438 t20 = sin(t19);
439 t25 = 0.2e1 * t9 * B;
440 t31 = 0.2e1 * t4 * t11;
441 t34 = -t3 * t14 * UU - t25 + t9 * UU + 0.3e1 * t3 * UU * t11 + t31 - 0.2e1 * t4 * t14;
442 t35 = cos(t19);
443 t44 = VV * VV;
444 t51 = -t10 - t13 + 0.3e1 * t3 * t44 * Rm - 0.4e1 * t4 * VV * Rm;
445 t61 = -t25 + t31 + t3 * t44 * VV - t9 * VV - 0.2e1 * t4 * t44 - 0.3e1 * t3 * VV * t11;
446 t70 = kn * kn;
447 t71 = t70 * kn;
448 t76 = B * B;
449 t82 = B * t8;
450 t86 = kn * x;
451 t87 = sin(t86);
452 t101 = cos(t86);
453 *p = ((t18 * t20 + t34 * t35) * _C1 + (t34 * t20 - t18 * t35) * _C2 + (t51 * t20 + t61 * t35) * _C3 + (t61 * t20 - t51 * t35) * _C4 + (-AA * t71 - 0.4e1 * BB * t70 * B + 0.4e1 * t76 * AA * kn - t8 * AA * kn - 0.4e1 * t82 * BB) * t87 + (-0.4e1 * t82 * AA - 0.4e1 * AA * t70 * B - 0.4e1 * t76 * BB * kn + t8 * BB * kn + BB * t71) * t101) / km;
454
455 t3 = exp(0.2e1 * x * B);
456 t4 = UU * UU;
457 t8 = km * km;
458 t9 = t3 * t8;
459 t10 = t9 * Rm;
460 t11 = Rm * Rm;
461 t13 = t3 * t11 * Rm;
462 t14 = t3 * B;
463 t18 = 0.3e1 * t3 * t4 * Rm + t10 - t13 + 0.4e1 * t14 * UU * Rm;
464 t19 = Rm * x;
465 t20 = sin(t19);
466 t23 = 0.2e1 * t9 * B;
467 t33 = 0.2e1 * t14 * t11;
468 t34 = -t23 + 0.3e1 * t3 * UU * t11 - t9 * UU - t3 * t4 * UU - 0.2e1 * t4 * t14 + t33;
469 t35 = cos(t19);
470 t47 = VV * VV;
471 t51 = t10 - 0.4e1 * t14 * VV * Rm + 0.3e1 * t3 * t47 * Rm - t13;
472 t61 = t9 * VV - t23 + t3 * t47 * VV - 0.2e1 * t14 * t47 + t33 - 0.3e1 * t3 * VV * t11;
473 t70 = B * B;
474 t74 = kn * kn;
475 t75 = t74 * kn;
476 t83 = kn * x;
477 t84 = sin(t83);
478 t96 = cos(t83);
479 *szz = -((t18 * t20 + t34 * t35) * _C1 + (t34 * t20 - t18 * t35) * _C2 + (t51 * t20 + t61 * t35) * _C3 + (t61 * t20 - t51 * t35) * _C4 + (0.4e1 * t70 * AA * kn - AA * t75 - 0.4e1 * BB * t74 * B + t8 * AA * kn) * t84 + (-t8 * BB * kn - 0.4e1 * AA * t74 * B - 0.4e1 * t70 * BB * kn + BB * t75) * t96) / km;
480
481 /* vx = Vx, vz = Vz, sxx = xx-component of stress tensor, sxz = xz-component of stress tensor, p = pressure, szz = zz-component of stress tensor */
482 *vx *= cos(km*z); /* Vx */
483 *vz *= sin(km*z); /* Vz */
484 *p *= cos(km*z); /* p */
485 *sxx *= cos(km*z); /* sxx total stress */
486 *sxz *= sin(km*z); /* tzx stress */
487 *szz *= cos(km*z); /* szz total stress */
488
489 /* rho = -sigma*sin(km*z)*cos(kn*x); */ /* density */
490 PetscFunctionReturn(0);
491 }
492
SolKxWrapperV(PetscInt dim,const PetscReal x[],PetscInt Nf,PetscScalar v[],void * ctx)493 PetscErrorCode SolKxWrapperV(PetscInt dim, const PetscReal x[], PetscInt Nf, PetscScalar v[], void *ctx)
494 {
495 PetscReal B = 100.0;
496 PetscReal kn = 100*M_PI;
497 PetscReal km = 100*M_PI;
498 PetscScalar p, sxx, sxz, szz;
499
500 PetscFunctionBeginUser;
501 SolKxSolution(x[0], x[1], kn, km, B, &v[0], &v[1], &p, &sxx, &sxz, &szz);
502 PetscFunctionReturn(0);
503 }
504
SolKxWrapperP(PetscInt dim,const PetscReal x[],PetscInt Nf,PetscScalar v[],void * ctx)505 PetscErrorCode SolKxWrapperP(PetscInt dim, const PetscReal x[], PetscInt Nf, PetscScalar v[], void *ctx)
506 {
507 PetscReal B = 100.0;
508 PetscReal kn = 100*M_PI;
509 PetscReal km = 100*M_PI;
510 PetscScalar vx, vz, sxx, sxz, szz;
511
512 PetscFunctionBeginUser;
513 SolKxSolution(x[0], x[1], kn, km, B, &vx, &vz, &v[0], &sxx, &sxz, &szz);
514 PetscFunctionReturn(0);
515 }
516
517 /*
518 Compare the C implementation with generated data from Maple
519 */
MapleTest(MPI_Comm comm,AppCtx * ctx)520 PetscErrorCode MapleTest(MPI_Comm comm, AppCtx *ctx)
521 {
522 const PetscInt n = 41;
523 PetscScalar vxMaple[41][41], vzMaple[41][41], pMaple[41][41], sxxMaple[41][41], sxzMaple[41][41], szzMaple[41][41];
524 PetscReal x[41], z[41];
525 PetscReal kn, km, B;
526 PetscInt i, j;
527 PetscErrorCode ierr;
528
529 PetscFunctionBegin;
530 ierr = SolKxData5(x, z, &kn, &km, &B, vxMaple, vzMaple, pMaple, sxxMaple, sxzMaple, szzMaple);CHKERRQ(ierr);
531 for (i = 0; i < n; ++i) {
532 for (j = 0; j < n; ++j) {
533 PetscScalar vx, vz, p, sxx, sxz, szz;
534 PetscReal norm;
535
536 ierr = SolKxSolution(x[i], z[j], kn, km, B, &vx, &vz, &p, &sxx, &sxz, &szz);CHKERRQ(ierr);
537 norm = sqrt(PetscSqr(PetscAbsScalar(vx - vxMaple[i][j])) + PetscSqr(PetscAbsScalar(vz - vzMaple[i][j])));
538 if (norm > 1.0e-10) {
539 ierr = PetscPrintf(PETSC_COMM_SELF, "%0.17e %0.17e %0.17e %0.17e %0.17e %0.17e %0.17e %0.17e %0.17e\n",
540 (double)x[i], (double)z[j], (double)PetscAbsScalar(vx - vxMaple[i][j]), (double)PetscAbsScalar(vz - vzMaple[i][j]), (double)PetscAbsScalar(p - pMaple[i][j]),
541 (double)PetscAbsScalar(sxx - sxxMaple[i][j]), (double)PetscAbsScalar(sxz - sxzMaple[i][j]), (double)PetscAbsScalar(szz - szzMaple[i][j]), (double)norm);
542 SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid solution, error %g", (double)norm);
543 }
544 }
545 }
546 ierr = PetscPrintf(comm, "Verified Maple test 5\n");CHKERRQ(ierr);
547 PetscFunctionReturn(0);
548 }
549
FEMTest(MPI_Comm comm,AppCtx * ctx)550 PetscErrorCode FEMTest(MPI_Comm comm, AppCtx *ctx)
551 {
552 DM dm;
553 Vec u;
554 PetscErrorCode (*funcs[2])(PetscInt, const PetscReal [], PetscInt, PetscScalar *, void *) = {SolKxWrapperV, SolKxWrapperP};
555 PetscReal discError;
556 PetscErrorCode ierr;
557
558 PetscFunctionBegin;
559 if (!ctx->fem) PetscFunctionReturn(0);
560 /* Create DM */
561 ierr = DMPlexCreateBoxMesh(comm, 2, PETSC_TRUE, NULL, NULL, NULL, NULL, PETSC_FALSE, &dm);CHKERRQ(ierr);
562 ierr = DMSetFromOptions(dm);CHKERRQ(ierr);
563 /* Project solution into FE space */
564 ierr = DMGetGlobalVector(dm, &u);CHKERRQ(ierr);
565 ierr = DMProjectFunction(dm, 0.0, funcs, NULL, INSERT_VALUES, u);CHKERRQ(ierr);
566 ierr = DMComputeL2Diff(dm, 0.0, funcs, NULL, u, &discError);CHKERRQ(ierr);
567 ierr = VecViewFromOptions(u, NULL, "-vec_view");CHKERRQ(ierr);
568 /* Cleanup */
569 ierr = DMRestoreGlobalVector(dm, &u);CHKERRQ(ierr);
570 ierr = DMDestroy(&dm);CHKERRQ(ierr);
571 PetscFunctionReturn(0);
572 }
573
main(int argc,char ** argv)574 int main(int argc, char **argv)
575 {
576 AppCtx user; /* user-defined work context */
577 PetscErrorCode ierr;
578
579 ierr = PetscInitialize(&argc, &argv, NULL,help);if (ierr) return ierr;
580 ierr = ProcessOptions(PETSC_COMM_WORLD, &user);CHKERRQ(ierr);
581 ierr = MapleTest(PETSC_COMM_WORLD, &user);CHKERRQ(ierr);
582 ierr = FEMTest(PETSC_COMM_WORLD, &user);CHKERRQ(ierr);
583 ierr = PetscFinalize();
584 return ierr;
585 }
586