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