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