Actual source code: ex16.c
1: static char help[] = "Large-deformation Elasticity Buckling Example";
3: /*F-----------------------------------------------------------------------
5: This example solves the 3D large deformation elasticity problem
7: \begin{equation}
8: \int_{\Omega}F \cdot S : \nabla v d\Omega + \int_{\Omega} (loading)\mathbf{e}_y \cdot v d\Omega = 0
9: \end{equation}
11: F is the deformation gradient, and S is the second Piola-Kirchhoff tensor from the Saint Venant-Kirchhoff model of
12: hyperelasticity. \Omega is a (arc) angle subsection of a cylindrical shell of thickness (height), inner radius
13: (rad) and width (width). The problem is discretized using Q1 finite elements on a logically structured grid.
14: Homogeneous Dirichlet boundary conditions are applied at the centers of the ends of the sphere.
16: This example is tunable with the following options:
17: -rad : the radius of the circle
18: -arc : set the angle of the arch represented
19: -loading : set the bulk loading (the mass)
20: -ploading : set the point loading at the top
21: -height : set the height of the arch
22: -width : set the width of the arch
23: -view_line : print initial and final offsets of the centerline of the
24: beam along the x direction
26: The material properties may be modified using either:
27: -mu : the first lame' parameter
28: -lambda : the second lame' parameter
30: Or:
31: -young : the Young's modulus
32: -poisson : the poisson ratio
34: This example is meant to show the strain placed upon the nonlinear solvers when trying to "snap through" the arch
35: using the loading. Under certain parameter regimes, the arch will invert under the load, and the number of Newton
36: steps will jump considerably. Composed nonlinear solvers may be used to mitigate this difficulty. This example
37: also demonstrates the use of the arc length continuation method NEWTONAL, which avoids the numerical difficulties
38: of the snap-through via tracing the equilibrium path through load increments.
40: The initial setup follows the example in pg. 268 of "Nonlinear Finite Element Methods" by Peter Wriggers, but is a
41: 3D extension.
43: ------------------------------------------------------------------------F*/
45: #include <petscsnes.h>
46: #include <petscdm.h>
47: #include <petscdmda.h>
49: #define QP0 0.2113248654051871
50: #define QP1 0.7886751345948129
51: #define NQ 2
52: #define NB 2
53: #define NEB 8
54: #define NEQ 8
55: #define NPB 24
57: #define NVALS NEB *NEQ
58: const PetscReal pts[NQ] = {QP0, QP1};
59: const PetscReal wts[NQ] = {0.5, 0.5};
61: PetscScalar vals[NVALS];
62: PetscScalar grad[3 * NVALS];
64: typedef PetscScalar Field[3];
65: typedef PetscScalar CoordField[3];
67: typedef PetscScalar JacField[9];
69: PetscErrorCode FormFunctionLocal(DMDALocalInfo *, Field ***, Field ***, void *);
70: PetscErrorCode FormJacobianLocal(DMDALocalInfo *, Field ***, Mat, Mat, void *);
71: PetscErrorCode DisplayLine(SNES, Vec);
72: PetscErrorCode FormElements(void);
74: typedef struct {
75: PetscReal loading;
76: PetscReal mu;
77: PetscReal lambda;
78: PetscReal rad;
79: PetscReal height;
80: PetscReal width;
81: PetscReal arc;
82: PetscReal ploading;
83: PetscReal load_factor;
84: } AppCtx;
86: PetscErrorCode InitialGuess(DM, AppCtx *, Vec);
87: PetscErrorCode ArcLengthScaling(DM, AppCtx *, Vec);
88: PetscErrorCode FormRHS(DM, AppCtx *, Vec);
89: PetscErrorCode FormCoordinates(DM, AppCtx *);
90: PetscErrorCode TangentLoad(SNES, Vec, Vec, void *);
92: int main(int argc, char **argv)
93: {
94: AppCtx user; /* user-defined work context */
95: PetscInt mx, my, its;
96: MPI_Comm comm;
97: SNES snes;
98: DM da;
99: Vec x, X, b;
100: PetscBool youngflg, poissonflg, muflg, lambdaflg, alflg, view = PETSC_FALSE, viewline = PETSC_FALSE, al_rescale = PETSC_FALSE;
101: PetscReal poisson = 0.2, young = 4e4;
102: char filename[PETSC_MAX_PATH_LEN] = "ex16.vts";
103: char filename_def[PETSC_MAX_PATH_LEN] = "ex16_def.vts";
105: PetscFunctionBeginUser;
106: PetscCall(PetscInitialize(&argc, &argv, NULL, help));
107: PetscCall(FormElements());
108: comm = PETSC_COMM_WORLD;
109: PetscCall(SNESCreate(comm, &snes));
110: PetscCall(DMDACreate3d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DMDA_STENCIL_BOX, 11, 2, 2, PETSC_DECIDE, PETSC_DECIDE, PETSC_DECIDE, 3, 1, NULL, NULL, NULL, &da));
111: PetscCall(DMSetFromOptions(da));
112: PetscCall(DMSetUp(da));
113: PetscCall(SNESSetDM(snes, (DM)da));
115: PetscCall(DMDAGetInfo(da, 0, &mx, &my, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE));
116: user.loading = 0.0;
117: user.arc = PETSC_PI / 3.;
118: user.mu = 4.0;
119: user.lambda = 1.0;
120: user.rad = 100.0;
121: user.height = 3.;
122: user.width = 1.;
123: user.ploading = -5e3;
124: user.load_factor = 1.0;
126: PetscCall(PetscOptionsGetReal(NULL, NULL, "-arc", &user.arc, NULL));
127: PetscCall(PetscOptionsGetReal(NULL, NULL, "-mu", &user.mu, &muflg));
128: PetscCall(PetscOptionsGetReal(NULL, NULL, "-lambda", &user.lambda, &lambdaflg));
129: PetscCall(PetscOptionsGetReal(NULL, NULL, "-rad", &user.rad, NULL));
130: PetscCall(PetscOptionsGetReal(NULL, NULL, "-height", &user.height, NULL));
131: PetscCall(PetscOptionsGetReal(NULL, NULL, "-width", &user.width, NULL));
132: PetscCall(PetscOptionsGetReal(NULL, NULL, "-loading", &user.loading, NULL));
133: PetscCall(PetscOptionsGetReal(NULL, NULL, "-ploading", &user.ploading, NULL));
134: PetscCall(PetscOptionsGetReal(NULL, NULL, "-poisson", &poisson, &poissonflg));
135: PetscCall(PetscOptionsGetReal(NULL, NULL, "-young", &young, &youngflg));
136: if (youngflg || poissonflg || !(muflg || lambdaflg)) {
137: /* set the lame' parameters based upon the poisson ratio and young's modulus */
138: user.lambda = poisson * young / ((1. + poisson) * (1. - 2. * poisson));
139: user.mu = young / (2. * (1. + poisson));
140: }
141: PetscCall(PetscOptionsGetBool(NULL, NULL, "-view", &view, NULL));
142: PetscCall(PetscOptionsGetBool(NULL, NULL, "-view_line", &viewline, NULL));
143: PetscCall(PetscOptionsGetBool(NULL, NULL, "-rescale", &al_rescale, NULL));
145: PetscCall(DMDASetFieldName(da, 0, "x_disp"));
146: PetscCall(DMDASetFieldName(da, 1, "y_disp"));
147: PetscCall(DMDASetFieldName(da, 2, "z_disp"));
149: PetscCall(DMSetApplicationContext(da, &user));
150: PetscCall(DMDASNESSetFunctionLocal(da, INSERT_VALUES, (DMDASNESFunctionFn *)FormFunctionLocal, &user));
151: PetscCall(DMDASNESSetJacobianLocal(da, (DMDASNESJacobianFn *)FormJacobianLocal, &user));
152: PetscCall(SNESSetFromOptions(snes));
153: PetscCall(SNESNewtonALSetFunction(snes, TangentLoad, &user));
154: PetscCall(PetscObjectTypeCompare((PetscObject)snes, SNESNEWTONAL, &alflg));
155: if (alflg) {
156: user.load_factor = 0.0;
158: if (al_rescale) {
159: Vec scaling;
161: PetscCall(DMCreateGlobalVector(da, &scaling));
162: PetscCall(ArcLengthScaling(da, &user, scaling));
163: PetscCall(VecPointwiseMult(scaling, scaling, scaling));
164: PetscCall(SNESNewtonALSetDiagonalScaling(snes, scaling));
165: PetscCall(VecDestroy(&scaling));
166: }
167: }
169: PetscCall(FormCoordinates(da, &user));
171: PetscCall(DMCreateGlobalVector(da, &x));
172: PetscCall(DMCreateGlobalVector(da, &b));
173: PetscCall(InitialGuess(da, &user, x));
174: PetscCall(FormRHS(da, &user, b));
176: PetscCall(PetscPrintf(comm, "lambda: %f mu: %f\n", (double)user.lambda, (double)user.mu));
178: /* show a cross-section of the initial state */
179: if (viewline) PetscCall(DisplayLine(snes, x));
181: /* get the loaded configuration */
182: PetscCall(SNESSolve(snes, b, x));
184: PetscCall(SNESGetIterationNumber(snes, &its));
185: PetscCall(PetscPrintf(comm, "Number of SNES iterations = %" PetscInt_FMT "\n", its));
186: PetscCall(SNESGetSolution(snes, &X));
187: /* show a cross-section of the final state */
188: if (viewline) PetscCall(DisplayLine(snes, X));
190: if (view) {
191: PetscViewer viewer;
192: Vec coords;
193: PetscCall(PetscViewerVTKOpen(PETSC_COMM_WORLD, filename, FILE_MODE_WRITE, &viewer));
194: PetscCall(VecView(x, viewer));
195: PetscCall(PetscViewerDestroy(&viewer));
196: PetscCall(DMGetCoordinates(da, &coords));
197: PetscCall(VecAXPY(coords, 1.0, x));
198: PetscCall(PetscViewerVTKOpen(PETSC_COMM_WORLD, filename_def, FILE_MODE_WRITE, &viewer));
199: PetscCall(VecView(x, viewer));
200: PetscCall(PetscViewerDestroy(&viewer));
201: }
203: PetscCall(VecDestroy(&x));
204: PetscCall(VecDestroy(&b));
205: PetscCall(DMDestroy(&da));
206: PetscCall(SNESDestroy(&snes));
207: PetscCall(PetscFinalize());
208: return 0;
209: }
211: PetscInt OnBoundary(PetscInt i, PetscInt j, PetscInt k, PetscInt mx, PetscInt my, PetscInt mz)
212: {
213: if ((i == 0 || i == mx - 1) && j == my - 1) return 1;
214: return 0;
215: }
217: void BoundaryValue(PetscInt i, PetscInt j, PetscInt k, PetscInt mx, PetscInt my, PetscInt mz, PetscScalar *val, AppCtx *user)
218: {
219: /* reference coordinates */
220: PetscReal p_x = ((PetscReal)i) / ((PetscReal)(mx - 1));
221: PetscReal p_y = ((PetscReal)j) / ((PetscReal)(my - 1));
222: PetscReal p_z = ((PetscReal)k) / ((PetscReal)(mz - 1));
223: PetscReal o_x = p_x;
224: PetscReal o_y = p_y;
225: PetscReal o_z = p_z;
226: val[0] = o_x - p_x;
227: val[1] = o_y - p_y;
228: val[2] = o_z - p_z;
229: }
231: void InvertTensor(PetscScalar *t, PetscScalar *ti, PetscReal *dett)
232: {
233: const PetscScalar a = t[0];
234: const PetscScalar b = t[1];
235: const PetscScalar c = t[2];
236: const PetscScalar d = t[3];
237: const PetscScalar e = t[4];
238: const PetscScalar f = t[5];
239: const PetscScalar g = t[6];
240: const PetscScalar h = t[7];
241: const PetscScalar i = t[8];
242: const PetscReal det = PetscRealPart(a * (e * i - f * h) - b * (i * d - f * g) + c * (d * h - e * g));
243: const PetscReal di = 1. / det;
244: if (ti) {
245: const PetscScalar A = (e * i - f * h);
246: const PetscScalar B = -(d * i - f * g);
247: const PetscScalar C = (d * h - e * g);
248: const PetscScalar D = -(b * i - c * h);
249: const PetscScalar E = (a * i - c * g);
250: const PetscScalar F = -(a * h - b * g);
251: const PetscScalar G = (b * f - c * e);
252: const PetscScalar H = -(a * f - c * d);
253: const PetscScalar II = (a * e - b * d);
254: ti[0] = di * A;
255: ti[1] = di * D;
256: ti[2] = di * G;
257: ti[3] = di * B;
258: ti[4] = di * E;
259: ti[5] = di * H;
260: ti[6] = di * C;
261: ti[7] = di * F;
262: ti[8] = di * II;
263: }
264: if (dett) *dett = det;
265: }
267: void TensorTensor(PetscScalar *a, PetscScalar *b, PetscScalar *c)
268: {
269: PetscInt i, j, m;
270: for (i = 0; i < 3; i++) {
271: for (j = 0; j < 3; j++) {
272: c[i + 3 * j] = 0;
273: for (m = 0; m < 3; m++) c[i + 3 * j] += a[m + 3 * j] * b[i + 3 * m];
274: }
275: }
276: }
278: void TensorTransposeTensor(PetscScalar *a, PetscScalar *b, PetscScalar *c)
279: {
280: PetscInt i, j, m;
281: for (i = 0; i < 3; i++) {
282: for (j = 0; j < 3; j++) {
283: c[i + 3 * j] = 0;
284: for (m = 0; m < 3; m++) c[i + 3 * j] += a[3 * m + j] * b[i + 3 * m];
285: }
286: }
287: }
289: void TensorVector(PetscScalar *rot, PetscScalar *vec, PetscScalar *tvec)
290: {
291: tvec[0] = rot[0] * vec[0] + rot[1] * vec[1] + rot[2] * vec[2];
292: tvec[1] = rot[3] * vec[0] + rot[4] * vec[1] + rot[5] * vec[2];
293: tvec[2] = rot[6] * vec[0] + rot[7] * vec[1] + rot[8] * vec[2];
294: }
296: void DeformationGradient(Field *ex, PetscInt qi, PetscInt qj, PetscInt qk, PetscScalar *invJ, PetscScalar *F)
297: {
298: PetscInt ii, jj, kk, l;
299: for (l = 0; l < 9; l++) F[l] = 0.;
300: F[0] = 1.;
301: F[4] = 1.;
302: F[8] = 1.;
303: /* form the deformation gradient at this basis function -- loop over element unknowns */
304: for (kk = 0; kk < NB; kk++) {
305: for (jj = 0; jj < NB; jj++) {
306: for (ii = 0; ii < NB; ii++) {
307: PetscInt idx = ii + jj * NB + kk * NB * NB;
308: PetscInt bidx = NEB * idx + qi + NQ * qj + NQ * NQ * qk;
309: PetscScalar lgrad[3];
310: TensorVector(invJ, &grad[3 * bidx], lgrad);
311: F[0] += lgrad[0] * ex[idx][0];
312: F[1] += lgrad[1] * ex[idx][0];
313: F[2] += lgrad[2] * ex[idx][0];
314: F[3] += lgrad[0] * ex[idx][1];
315: F[4] += lgrad[1] * ex[idx][1];
316: F[5] += lgrad[2] * ex[idx][1];
317: F[6] += lgrad[0] * ex[idx][2];
318: F[7] += lgrad[1] * ex[idx][2];
319: F[8] += lgrad[2] * ex[idx][2];
320: }
321: }
322: }
323: }
325: void DeformationGradientJacobian(PetscInt qi, PetscInt qj, PetscInt qk, PetscInt ii, PetscInt jj, PetscInt kk, PetscInt fld, PetscScalar *invJ, PetscScalar *dF)
326: {
327: PetscInt l;
328: PetscScalar lgrad[3];
329: PetscInt idx = ii + jj * NB + kk * NB * NB;
330: PetscInt bidx = NEB * idx + qi + NQ * qj + NQ * NQ * qk;
331: for (l = 0; l < 9; l++) dF[l] = 0.;
332: /* form the deformation gradient at this basis function -- loop over element unknowns */
333: TensorVector(invJ, &grad[3 * bidx], lgrad);
334: dF[3 * fld] = lgrad[0];
335: dF[3 * fld + 1] = lgrad[1];
336: dF[3 * fld + 2] = lgrad[2];
337: }
339: void LagrangeGreenStrain(PetscScalar *F, PetscScalar *E)
340: {
341: PetscInt i, j, m;
342: for (i = 0; i < 3; i++) {
343: for (j = 0; j < 3; j++) {
344: E[i + 3 * j] = 0;
345: for (m = 0; m < 3; m++) E[i + 3 * j] += 0.5 * F[3 * m + j] * F[i + 3 * m];
346: }
347: }
348: for (i = 0; i < 3; i++) E[i + 3 * i] -= 0.5;
349: }
351: void SaintVenantKirchoff(PetscReal lambda, PetscReal mu, PetscScalar *F, PetscScalar *S)
352: {
353: PetscInt i, j;
354: PetscScalar E[9];
355: PetscScalar trE = 0;
356: LagrangeGreenStrain(F, E);
357: for (i = 0; i < 3; i++) trE += E[i + 3 * i];
358: for (i = 0; i < 3; i++) {
359: for (j = 0; j < 3; j++) {
360: S[i + 3 * j] = 2. * mu * E[i + 3 * j];
361: if (i == j) S[i + 3 * j] += trE * lambda;
362: }
363: }
364: }
366: void SaintVenantKirchoffJacobian(PetscReal lambda, PetscReal mu, PetscScalar *F, PetscScalar *dF, PetscScalar *dS)
367: {
368: PetscScalar FtdF[9], dE[9];
369: PetscInt i, j;
370: PetscScalar dtrE = 0.;
371: TensorTransposeTensor(dF, F, dE);
372: TensorTransposeTensor(F, dF, FtdF);
373: for (i = 0; i < 9; i++) dE[i] += FtdF[i];
374: for (i = 0; i < 9; i++) dE[i] *= 0.5;
375: for (i = 0; i < 3; i++) dtrE += dE[i + 3 * i];
376: for (i = 0; i < 3; i++) {
377: for (j = 0; j < 3; j++) {
378: dS[i + 3 * j] = 2. * mu * dE[i + 3 * j];
379: if (i == j) dS[i + 3 * j] += dtrE * lambda;
380: }
381: }
382: }
384: PetscErrorCode FormElements(void)
385: {
386: PetscInt i, j, k, ii, jj, kk;
387: PetscReal bx, by, bz, dbx, dby, dbz;
389: PetscFunctionBeginUser;
390: /* construct the basis function values and derivatives */
391: for (k = 0; k < NB; k++) {
392: for (j = 0; j < NB; j++) {
393: for (i = 0; i < NB; i++) {
394: /* loop over the quadrature points */
395: for (kk = 0; kk < NQ; kk++) {
396: for (jj = 0; jj < NQ; jj++) {
397: for (ii = 0; ii < NQ; ii++) {
398: PetscInt idx = ii + NQ * jj + NQ * NQ * kk + NEQ * i + NEQ * NB * j + NEQ * NB * NB * k;
399: bx = pts[ii];
400: by = pts[jj];
401: bz = pts[kk];
402: dbx = 1.;
403: dby = 1.;
404: dbz = 1.;
405: if (i == 0) {
406: bx = 1. - bx;
407: dbx = -1;
408: }
409: if (j == 0) {
410: by = 1. - by;
411: dby = -1;
412: }
413: if (k == 0) {
414: bz = 1. - bz;
415: dbz = -1;
416: }
417: vals[idx] = bx * by * bz;
418: grad[3 * idx + 0] = dbx * by * bz;
419: grad[3 * idx + 1] = dby * bx * bz;
420: grad[3 * idx + 2] = dbz * bx * by;
421: }
422: }
423: }
424: }
425: }
426: }
427: PetscFunctionReturn(PETSC_SUCCESS);
428: }
430: void GatherElementData(PetscInt mx, PetscInt my, PetscInt mz, Field ***x, CoordField ***c, PetscInt i, PetscInt j, PetscInt k, Field *ex, CoordField *ec, AppCtx *user)
431: {
432: PetscInt m;
433: PetscInt ii, jj, kk;
434: /* gather the data -- loop over element unknowns */
435: for (kk = 0; kk < NB; kk++) {
436: for (jj = 0; jj < NB; jj++) {
437: for (ii = 0; ii < NB; ii++) {
438: PetscInt idx = ii + jj * NB + kk * NB * NB;
439: /* decouple the boundary nodes for the displacement variables */
440: if (OnBoundary(i + ii, j + jj, k + kk, mx, my, mz)) {
441: BoundaryValue(i + ii, j + jj, k + kk, mx, my, mz, ex[idx], user);
442: } else {
443: for (m = 0; m < 3; m++) ex[idx][m] = x[k + kk][j + jj][i + ii][m];
444: }
445: for (m = 0; m < 3; m++) ec[idx][m] = c[k + kk][j + jj][i + ii][m];
446: }
447: }
448: }
449: }
451: void QuadraturePointGeometricJacobian(CoordField *ec, PetscInt qi, PetscInt qj, PetscInt qk, PetscScalar *J)
452: {
453: PetscInt ii, jj, kk;
454: /* construct the gradient at the given quadrature point named by i,j,k */
455: for (ii = 0; ii < 9; ii++) J[ii] = 0;
456: for (kk = 0; kk < NB; kk++) {
457: for (jj = 0; jj < NB; jj++) {
458: for (ii = 0; ii < NB; ii++) {
459: PetscInt idx = ii + jj * NB + kk * NB * NB;
460: PetscInt bidx = NEB * idx + qi + NQ * qj + NQ * NQ * qk;
461: J[0] += grad[3 * bidx + 0] * ec[idx][0];
462: J[1] += grad[3 * bidx + 1] * ec[idx][0];
463: J[2] += grad[3 * bidx + 2] * ec[idx][0];
464: J[3] += grad[3 * bidx + 0] * ec[idx][1];
465: J[4] += grad[3 * bidx + 1] * ec[idx][1];
466: J[5] += grad[3 * bidx + 2] * ec[idx][1];
467: J[6] += grad[3 * bidx + 0] * ec[idx][2];
468: J[7] += grad[3 * bidx + 1] * ec[idx][2];
469: J[8] += grad[3 * bidx + 2] * ec[idx][2];
470: }
471: }
472: }
473: }
475: void FormElementJacobian(Field *ex, CoordField *ec, Field *ef, Field *eq, PetscScalar *ej, AppCtx *user)
476: {
477: PetscReal vol;
478: PetscScalar J[9];
479: PetscScalar invJ[9];
480: PetscScalar F[9], S[9], dF[9], dS[9], dFS[9], FdS[9], FS[9];
481: PetscReal scl;
482: PetscInt i, j, k, l, ii, jj, kk, ll, qi, qj, qk, m;
484: if (ej)
485: for (i = 0; i < NPB * NPB; i++) ej[i] = 0.;
486: if (ef)
487: for (i = 0; i < NEB; i++) {
488: ef[i][0] = 0.;
489: ef[i][1] = 0.;
490: ef[i][2] = 0.;
491: }
492: if (eq)
493: for (i = 0; i < NEB; i++) {
494: eq[i][0] = 0.;
495: eq[i][1] = 0.;
496: eq[i][2] = 0.;
497: }
498: /* loop over quadrature */
499: for (qk = 0; qk < NQ; qk++) {
500: for (qj = 0; qj < NQ; qj++) {
501: for (qi = 0; qi < NQ; qi++) {
502: QuadraturePointGeometricJacobian(ec, qi, qj, qk, J);
503: InvertTensor(J, invJ, &vol);
504: scl = vol * wts[qi] * wts[qj] * wts[qk];
505: DeformationGradient(ex, qi, qj, qk, invJ, F);
506: SaintVenantKirchoff(user->lambda, user->mu, F, S);
507: /* form the function */
508: if (ef) {
509: TensorTensor(F, S, FS);
510: for (kk = 0; kk < NB; kk++) {
511: for (jj = 0; jj < NB; jj++) {
512: for (ii = 0; ii < NB; ii++) {
513: PetscInt idx = ii + jj * NB + kk * NB * NB;
514: PetscInt bidx = NEB * idx + qi + NQ * qj + NQ * NQ * qk;
515: PetscScalar lgrad[3];
516: TensorVector(invJ, &grad[3 * bidx], lgrad);
517: /* mu*F : grad phi_{u,v,w} */
518: for (m = 0; m < 3; m++) ef[idx][m] += scl * (lgrad[0] * FS[3 * m + 0] + lgrad[1] * FS[3 * m + 1] + lgrad[2] * FS[3 * m + 2]);
519: ef[idx][1] -= user->load_factor * scl * user->loading * vals[bidx];
520: }
521: }
522: }
523: }
524: if (eq) {
525: for (kk = 0; kk < NB; kk++) {
526: for (jj = 0; jj < NB; jj++) {
527: for (ii = 0; ii < NB; ii++) {
528: PetscInt idx = ii + jj * NB + kk * NB * NB;
529: PetscInt bidx = NEB * idx + qi + NQ * qj + NQ * NQ * qk;
530: /* external force vector */
531: eq[idx][1] += scl * user->loading * vals[bidx];
532: }
533: }
534: }
535: }
536: /* form the jacobian */
537: if (ej) {
538: /* loop over trialfunctions */
539: for (k = 0; k < NB; k++) {
540: for (j = 0; j < NB; j++) {
541: for (i = 0; i < NB; i++) {
542: for (l = 0; l < 3; l++) {
543: PetscInt tridx = l + 3 * (i + j * NB + k * NB * NB);
544: DeformationGradientJacobian(qi, qj, qk, i, j, k, l, invJ, dF);
545: SaintVenantKirchoffJacobian(user->lambda, user->mu, F, dF, dS);
546: TensorTensor(dF, S, dFS);
547: TensorTensor(F, dS, FdS);
548: for (m = 0; m < 9; m++) dFS[m] += FdS[m];
549: /* loop over testfunctions */
550: for (kk = 0; kk < NB; kk++) {
551: for (jj = 0; jj < NB; jj++) {
552: for (ii = 0; ii < NB; ii++) {
553: PetscInt idx = ii + jj * NB + kk * NB * NB;
554: PetscInt bidx = 8 * idx + qi + NQ * qj + NQ * NQ * qk;
555: PetscScalar lgrad[3];
556: TensorVector(invJ, &grad[3 * bidx], lgrad);
557: for (ll = 0; ll < 3; ll++) {
558: PetscInt teidx = ll + 3 * (ii + jj * NB + kk * NB * NB);
559: ej[teidx + NPB * tridx] += scl * (lgrad[0] * dFS[3 * ll + 0] + lgrad[1] * dFS[3 * ll + 1] + lgrad[2] * dFS[3 * ll + 2]);
560: }
561: }
562: }
563: } /* end of testfunctions */
564: }
565: }
566: }
567: } /* end of trialfunctions */
568: }
569: }
570: }
571: } /* end of quadrature points */
572: }
574: void ApplyBCsElement(PetscInt mx, PetscInt my, PetscInt mz, PetscInt i, PetscInt j, PetscInt k, PetscScalar *jacobian)
575: {
576: PetscInt ii, jj, kk, ll, ei, ej, ek, el;
577: for (kk = 0; kk < NB; kk++) {
578: for (jj = 0; jj < NB; jj++) {
579: for (ii = 0; ii < NB; ii++) {
580: for (ll = 0; ll < 3; ll++) {
581: PetscInt tridx = ll + 3 * (ii + jj * NB + kk * NB * NB);
582: for (ek = 0; ek < NB; ek++) {
583: for (ej = 0; ej < NB; ej++) {
584: for (ei = 0; ei < NB; ei++) {
585: for (el = 0; el < 3; el++) {
586: if (OnBoundary(i + ii, j + jj, k + kk, mx, my, mz) || OnBoundary(i + ei, j + ej, k + ek, mx, my, mz)) {
587: PetscInt teidx = el + 3 * (ei + ej * NB + ek * NB * NB);
588: if (teidx == tridx) {
589: jacobian[tridx + NPB * teidx] = 1.;
590: } else {
591: jacobian[tridx + NPB * teidx] = 0.;
592: }
593: }
594: }
595: }
596: }
597: }
598: }
599: }
600: }
601: }
602: }
604: PetscErrorCode FormJacobianLocal(DMDALocalInfo *info, Field ***x, Mat jacpre, Mat jac, void *ptr)
605: {
606: /* values for each basis function at each quadrature point */
607: AppCtx *user = (AppCtx *)ptr;
608: PetscInt i, j, k, m, l;
609: PetscInt ii, jj, kk;
610: PetscScalar ej[NPB * NPB];
611: PetscScalar vals[NPB * NPB];
612: Field ex[NEB];
613: CoordField ec[NEB];
615: PetscInt xs = info->xs, ys = info->ys, zs = info->zs;
616: PetscInt xm = info->xm, ym = info->ym, zm = info->zm;
617: PetscInt xes, yes, zes, xee, yee, zee;
618: PetscInt mx = info->mx, my = info->my, mz = info->mz;
619: DM cda;
620: CoordField ***c;
621: Vec C;
622: PetscInt nrows;
623: MatStencil col[NPB], row[NPB];
624: PetscScalar v[9];
626: PetscFunctionBeginUser;
627: PetscCall(DMGetCoordinateDM(info->da, &cda));
628: PetscCall(DMGetCoordinatesLocal(info->da, &C));
629: PetscCall(DMDAVecGetArray(cda, C, &c));
630: PetscCall(MatScale(jac, 0.0));
632: xes = xs;
633: yes = ys;
634: zes = zs;
635: xee = xs + xm;
636: yee = ys + ym;
637: zee = zs + zm;
638: if (xs > 0) xes = xs - 1;
639: if (ys > 0) yes = ys - 1;
640: if (zs > 0) zes = zs - 1;
641: if (xs + xm == mx) xee = xs + xm - 1;
642: if (ys + ym == my) yee = ys + ym - 1;
643: if (zs + zm == mz) zee = zs + zm - 1;
644: for (k = zes; k < zee; k++) {
645: for (j = yes; j < yee; j++) {
646: for (i = xes; i < xee; i++) {
647: GatherElementData(mx, my, mz, x, c, i, j, k, ex, ec, user);
648: FormElementJacobian(ex, ec, NULL, NULL, ej, user);
649: ApplyBCsElement(mx, my, mz, i, j, k, ej);
650: nrows = 0.;
651: for (kk = 0; kk < NB; kk++) {
652: for (jj = 0; jj < NB; jj++) {
653: for (ii = 0; ii < NB; ii++) {
654: PetscInt idx = ii + jj * 2 + kk * 4;
655: for (m = 0; m < 3; m++) {
656: col[3 * idx + m].i = i + ii;
657: col[3 * idx + m].j = j + jj;
658: col[3 * idx + m].k = k + kk;
659: col[3 * idx + m].c = m;
660: if (i + ii >= xs && i + ii < xm + xs && j + jj >= ys && j + jj < ys + ym && k + kk >= zs && k + kk < zs + zm) {
661: row[nrows].i = i + ii;
662: row[nrows].j = j + jj;
663: row[nrows].k = k + kk;
664: row[nrows].c = m;
665: for (l = 0; l < NPB; l++) vals[NPB * nrows + l] = ej[NPB * (3 * idx + m) + l];
666: nrows++;
667: }
668: }
669: }
670: }
671: }
672: PetscCall(MatSetValuesStencil(jac, nrows, row, NPB, col, vals, ADD_VALUES));
673: }
674: }
675: }
677: PetscCall(MatAssemblyBegin(jac, MAT_FLUSH_ASSEMBLY));
678: PetscCall(MatAssemblyEnd(jac, MAT_FLUSH_ASSEMBLY));
680: /* set the diagonal */
681: v[0] = 1.;
682: v[1] = 0.;
683: v[2] = 0.;
684: v[3] = 0.;
685: v[4] = 1.;
686: v[5] = 0.;
687: v[6] = 0.;
688: v[7] = 0.;
689: v[8] = 1.;
690: for (k = zs; k < zs + zm; k++) {
691: for (j = ys; j < ys + ym; j++) {
692: for (i = xs; i < xs + xm; i++) {
693: if (OnBoundary(i, j, k, mx, my, mz)) {
694: for (m = 0; m < 3; m++) {
695: col[m].i = i;
696: col[m].j = j;
697: col[m].k = k;
698: col[m].c = m;
699: }
700: PetscCall(MatSetValuesStencil(jac, 3, col, 3, col, v, INSERT_VALUES));
701: }
702: }
703: }
704: }
706: PetscCall(MatAssemblyBegin(jac, MAT_FINAL_ASSEMBLY));
707: PetscCall(MatAssemblyEnd(jac, MAT_FINAL_ASSEMBLY));
709: PetscCall(DMDAVecRestoreArray(cda, C, &c));
710: PetscFunctionReturn(PETSC_SUCCESS);
711: }
713: PetscErrorCode FormFunctionLocal(DMDALocalInfo *info, Field ***x, Field ***f, void *ptr)
714: {
715: /* values for each basis function at each quadrature point */
716: AppCtx *user = (AppCtx *)ptr;
717: PetscInt i, j, k, l;
718: PetscInt ii, jj, kk;
720: Field ef[NEB];
721: Field ex[NEB];
722: CoordField ec[NEB];
724: PetscInt xs = info->xs, ys = info->ys, zs = info->zs;
725: PetscInt xm = info->xm, ym = info->ym, zm = info->zm;
726: PetscInt xes, yes, zes, xee, yee, zee;
727: PetscInt mx = info->mx, my = info->my, mz = info->mz;
728: DM cda;
729: CoordField ***c;
730: Vec C;
732: PetscFunctionBeginUser;
733: PetscCall(DMGetCoordinateDM(info->da, &cda));
734: PetscCall(DMGetCoordinatesLocal(info->da, &C));
735: PetscCall(DMDAVecGetArray(cda, C, &c));
736: PetscCall(DMDAGetInfo(info->da, 0, &mx, &my, &mz, 0, 0, 0, 0, 0, 0, 0, 0, 0));
737: PetscCall(DMDAGetCorners(info->da, &xs, &ys, &zs, &xm, &ym, &zm));
739: /* loop over elements */
740: for (k = zs; k < zs + zm; k++) {
741: for (j = ys; j < ys + ym; j++) {
742: for (i = xs; i < xs + xm; i++) {
743: for (l = 0; l < 3; l++) f[k][j][i][l] = 0.;
744: }
745: }
746: }
747: /* element starts and ends */
748: xes = xs;
749: yes = ys;
750: zes = zs;
751: xee = xs + xm;
752: yee = ys + ym;
753: zee = zs + zm;
754: if (xs > 0) xes = xs - 1;
755: if (ys > 0) yes = ys - 1;
756: if (zs > 0) zes = zs - 1;
757: if (xs + xm == mx) xee = xs + xm - 1;
758: if (ys + ym == my) yee = ys + ym - 1;
759: if (zs + zm == mz) zee = zs + zm - 1;
760: for (k = zes; k < zee; k++) {
761: for (j = yes; j < yee; j++) {
762: for (i = xes; i < xee; i++) {
763: GatherElementData(mx, my, mz, x, c, i, j, k, ex, ec, user);
764: FormElementJacobian(ex, ec, ef, NULL, NULL, user);
765: /* put this element's additions into the residuals */
766: for (kk = 0; kk < NB; kk++) {
767: for (jj = 0; jj < NB; jj++) {
768: for (ii = 0; ii < NB; ii++) {
769: PetscInt idx = ii + jj * NB + kk * NB * NB;
770: if (k + kk >= zs && j + jj >= ys && i + ii >= xs && k + kk < zs + zm && j + jj < ys + ym && i + ii < xs + xm) {
771: if (OnBoundary(i + ii, j + jj, k + kk, mx, my, mz)) {
772: for (l = 0; l < 3; l++) f[k + kk][j + jj][i + ii][l] = x[k + kk][j + jj][i + ii][l] - ex[idx][l];
773: } else {
774: for (l = 0; l < 3; l++) f[k + kk][j + jj][i + ii][l] += ef[idx][l];
775: }
776: }
777: }
778: }
779: }
780: }
781: }
782: }
783: PetscCall(DMDAVecRestoreArray(cda, C, &c));
784: PetscFunctionReturn(PETSC_SUCCESS);
785: }
787: PetscErrorCode TangentLoad(SNES snes, Vec X, Vec Q, void *ptr)
788: {
789: /* values for each basis function at each quadrature point */
790: AppCtx *user = (AppCtx *)ptr;
791: PetscInt xs, ys, zs;
792: PetscInt xm, ym, zm;
793: PetscInt mx, my, mz;
794: DM da;
795: Vec Xl, Ql;
796: Field ***x, ***q;
797: PetscInt i, j, k, l;
798: PetscInt ii, jj, kk;
800: Field eq[NEB];
801: Field ex[NEB];
802: CoordField ec[NEB];
804: PetscInt xes, yes, zes, xee, yee, zee;
805: DM cda;
806: CoordField ***c;
807: Vec C;
809: PetscFunctionBeginUser;
810: /* update user context with current load parameter */
811: PetscCall(SNESNewtonALGetLoadParameter(snes, &user->load_factor));
813: PetscCall(SNESGetDM(snes, &da));
814: PetscCall(DMGetLocalVector(da, &Xl));
815: PetscCall(DMGetLocalVector(da, &Ql));
816: PetscCall(DMGlobalToLocal(da, X, INSERT_VALUES, Xl));
818: PetscCall(DMDAVecGetArray(da, Xl, &x));
819: PetscCall(DMDAVecGetArray(da, Ql, &q));
821: PetscCall(DMGetCoordinateDM(da, &cda));
822: PetscCall(DMGetCoordinatesLocal(da, &C));
823: PetscCall(DMDAVecGetArray(cda, C, &c));
824: PetscCall(DMDAGetInfo(da, 0, &mx, &my, &mz, 0, 0, 0, 0, 0, 0, 0, 0, 0));
825: PetscCall(DMDAGetCorners(da, &xs, &ys, &zs, &xm, &ym, &zm));
827: /* loop over elements */
828: for (k = zs; k < zs + zm; k++) {
829: for (j = ys; j < ys + ym; j++) {
830: for (i = xs; i < xs + xm; i++) {
831: for (l = 0; l < 3; l++) q[k][j][i][l] = 0.;
832: }
833: }
834: }
835: /* element starts and ends */
836: xes = xs;
837: yes = ys;
838: zes = zs;
839: xee = xs + xm;
840: yee = ys + ym;
841: zee = zs + zm;
842: if (xs > 0) xes = xs - 1;
843: if (ys > 0) yes = ys - 1;
844: if (zs > 0) zes = zs - 1;
845: if (xs + xm == mx) xee = xs + xm - 1;
846: if (ys + ym == my) yee = ys + ym - 1;
847: if (zs + zm == mz) zee = zs + zm - 1;
848: for (k = zes; k < zee; k++) {
849: for (j = yes; j < yee; j++) {
850: for (i = xes; i < xee; i++) {
851: GatherElementData(mx, my, mz, x, c, i, j, k, ex, ec, user);
852: FormElementJacobian(ex, ec, NULL, eq, NULL, user);
853: /* put this element's additions into the residuals */
854: for (kk = 0; kk < NB; kk++) {
855: for (jj = 0; jj < NB; jj++) {
856: for (ii = 0; ii < NB; ii++) {
857: PetscInt idx = ii + jj * NB + kk * NB * NB;
858: if (k + kk >= zs && j + jj >= ys && i + ii >= xs && k + kk < zs + zm && j + jj < ys + ym && i + ii < xs + xm) {
859: if (!OnBoundary(i + ii, j + jj, k + kk, mx, my, mz)) {
860: for (l = 0; l < 3; l++) q[k + kk][j + jj][i + ii][l] += eq[idx][l];
861: }
862: }
863: }
864: }
865: }
866: }
867: }
868: }
869: PetscCall(DMDAVecRestoreArray(da, Xl, &x));
870: PetscCall(DMDAVecRestoreArray(da, Ql, &q));
871: PetscCall(VecZeroEntries(Q));
872: PetscCall(DMLocalToGlobal(da, Ql, INSERT_VALUES, Q));
873: PetscCall(DMRestoreLocalVector(da, &Ql));
874: PetscCall(DMRestoreLocalVector(da, &Xl));
875: PetscCall(DMDAVecRestoreArray(cda, C, &c));
876: PetscFunctionReturn(PETSC_SUCCESS);
877: }
879: PetscErrorCode FormCoordinates(DM da, AppCtx *user)
880: {
881: Vec coords;
882: DM cda;
883: PetscInt mx, my, mz;
884: PetscInt i, j, k, xs, ys, zs, xm, ym, zm;
885: CoordField ***x;
887: PetscFunctionBeginUser;
888: PetscCall(DMGetCoordinateDM(da, &cda));
889: PetscCall(DMCreateGlobalVector(cda, &coords));
890: PetscCall(DMDAGetInfo(da, 0, &mx, &my, &mz, 0, 0, 0, 0, 0, 0, 0, 0, 0));
891: PetscCall(DMDAGetCorners(da, &xs, &ys, &zs, &xm, &ym, &zm));
892: PetscCall(DMDAVecGetArray(da, coords, &x));
893: for (k = zs; k < zs + zm; k++) {
894: for (j = ys; j < ys + ym; j++) {
895: for (i = xs; i < xs + xm; i++) {
896: PetscReal cx = ((PetscReal)i) / ((PetscReal)(mx - 1));
897: PetscReal cy = ((PetscReal)j) / ((PetscReal)(my - 1));
898: PetscReal cz = ((PetscReal)k) / ((PetscReal)(mz - 1));
899: PetscReal rad = user->rad + cy * user->height;
900: PetscReal ang = (cx - 0.5) * user->arc;
901: x[k][j][i][0] = rad * PetscSinReal(ang);
902: x[k][j][i][1] = rad * PetscCosReal(ang) - (user->rad + 0.5 * user->height) * PetscCosReal(-0.5 * user->arc);
903: x[k][j][i][2] = user->width * (cz - 0.5);
904: }
905: }
906: }
907: PetscCall(DMDAVecRestoreArray(da, coords, &x));
908: PetscCall(DMSetCoordinates(da, coords));
909: PetscCall(VecDestroy(&coords));
910: PetscFunctionReturn(PETSC_SUCCESS);
911: }
913: PetscErrorCode InitialGuess(DM da, AppCtx *user, Vec X)
914: {
915: PetscInt i, j, k, xs, ys, zs, xm, ym, zm;
916: PetscInt mx, my, mz;
917: Field ***x;
919: PetscFunctionBeginUser;
920: PetscCall(DMDAGetCorners(da, &xs, &ys, &zs, &xm, &ym, &zm));
921: PetscCall(DMDAGetInfo(da, 0, &mx, &my, &mz, 0, 0, 0, 0, 0, 0, 0, 0, 0));
922: PetscCall(DMDAVecGetArray(da, X, &x));
924: for (k = zs; k < zs + zm; k++) {
925: for (j = ys; j < ys + ym; j++) {
926: for (i = xs; i < xs + xm; i++) {
927: /* reference coordinates */
928: PetscReal p_x = ((PetscReal)i) / ((PetscReal)(mx - 1));
929: PetscReal p_y = ((PetscReal)j) / ((PetscReal)(my - 1));
930: PetscReal p_z = ((PetscReal)k) / ((PetscReal)(mz - 1));
931: PetscReal o_x = p_x;
932: PetscReal o_y = p_y;
933: PetscReal o_z = p_z;
934: x[k][j][i][0] = o_x - p_x;
935: x[k][j][i][1] = o_y - p_y;
936: x[k][j][i][2] = o_z - p_z;
937: }
938: }
939: }
940: PetscCall(DMDAVecRestoreArray(da, X, &x));
941: PetscFunctionReturn(PETSC_SUCCESS);
942: }
944: PetscErrorCode ArcLengthScaling(DM da, AppCtx *user, Vec V)
945: {
946: PetscInt i, j, k, xs, ys, zs, xm, ym, zm;
947: PetscInt mx, my, mz;
948: Field ***v;
949: PetscReal rad = user->rad + user->height;
951: PetscFunctionBeginUser;
952: PetscCall(DMDAGetCorners(da, &xs, &ys, &zs, &xm, &ym, &zm));
953: PetscCall(DMDAGetInfo(da, 0, &mx, &my, &mz, 0, 0, 0, 0, 0, 0, 0, 0, 0));
954: PetscCall(DMDAVecGetArray(da, V, &v));
956: for (k = zs; k < zs + zm; k++) {
957: for (j = ys; j < ys + ym; j++) {
958: for (i = xs; i < xs + xm; i++) {
959: v[k][j][i][0] = 2. / rad;
960: v[k][j][i][1] = 2. / rad;
961: v[k][j][i][2] = 2. / user->width;
962: }
963: }
964: }
965: PetscCall(DMDAVecRestoreArray(da, V, &v));
966: PetscFunctionReturn(PETSC_SUCCESS);
967: }
969: PetscErrorCode FormRHS(DM da, AppCtx *user, Vec X)
970: {
971: PetscInt i, j, k, xs, ys, zs, xm, ym, zm;
972: PetscInt mx, my, mz;
973: Field ***x;
975: PetscFunctionBeginUser;
976: PetscCall(DMDAGetCorners(da, &xs, &ys, &zs, &xm, &ym, &zm));
977: PetscCall(DMDAGetInfo(da, 0, &mx, &my, &mz, 0, 0, 0, 0, 0, 0, 0, 0, 0));
978: PetscCall(DMDAVecGetArray(da, X, &x));
980: for (k = zs; k < zs + zm; k++) {
981: for (j = ys; j < ys + ym; j++) {
982: for (i = xs; i < xs + xm; i++) {
983: x[k][j][i][0] = 0.;
984: x[k][j][i][1] = 0.;
985: x[k][j][i][2] = 0.;
986: if (i == (mx - 1) / 2 && j == (my - 1)) x[k][j][i][1] = user->ploading / (mz - 1);
987: }
988: }
989: }
990: PetscCall(DMDAVecRestoreArray(da, X, &x));
991: PetscFunctionReturn(PETSC_SUCCESS);
992: }
994: PetscErrorCode DisplayLine(SNES snes, Vec X)
995: {
996: PetscInt r, i, j = 0, k = 0, xs, xm, ys, ym, zs, zm, mx, my, mz;
997: Field ***x;
998: CoordField ***c;
999: DM da, cda;
1000: Vec C;
1001: PetscMPIInt size, rank;
1003: PetscFunctionBegin;
1004: PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
1005: PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
1006: PetscCall(SNESGetDM(snes, &da));
1007: PetscCall(DMDAGetInfo(da, 0, &mx, &my, &mz, 0, 0, 0, 0, 0, 0, 0, 0, 0));
1008: PetscCall(DMGetCoordinateDM(da, &cda));
1009: PetscCall(DMGetCoordinates(da, &C));
1010: j = my / 2;
1011: k = mz / 2;
1012: PetscCall(DMDAGetCorners(da, &xs, &ys, &zs, &xm, &ym, &zm));
1013: PetscCall(DMDAVecGetArray(da, X, &x));
1014: PetscCall(DMDAVecGetArray(cda, C, &c));
1015: for (r = 0; r < size; r++) {
1016: if (rank == r) {
1017: if (j >= ys && j < ys + ym && k >= zs && k < zs + zm) {
1018: for (i = xs; i < xs + xm; i++) {
1019: PetscCall(PetscPrintf(PETSC_COMM_SELF, "%" PetscInt_FMT " %d %d: %f %f %f\n", i, 0, 0, (double)PetscRealPart(c[k][j][i][0] + x[k][j][i][0]), (double)PetscRealPart(c[k][j][i][1] + x[k][j][i][1]), (double)PetscRealPart(c[k][j][i][2] + x[k][j][i][2])));
1020: }
1021: }
1022: }
1023: PetscCallMPI(MPI_Barrier(PETSC_COMM_WORLD));
1024: }
1025: PetscCall(DMDAVecRestoreArray(da, X, &x));
1026: PetscCall(DMDAVecRestoreArray(cda, C, &c));
1027: PetscFunctionReturn(PETSC_SUCCESS);
1028: }
1030: /*TEST
1032: test:
1033: nsize: 2
1034: args: -da_refine 2 -pc_type mg -rad 10.0 -young 10. -ploading 0.0 -loading -1. -mg_levels_ksp_max_it 2 -snes_monitor_short -ksp_monitor_short
1035: requires: !single
1036: timeoutfactor: 3
1038: test:
1039: suffix: 2
1040: args: -da_refine 2 -pc_type mg -rad 10.0 -young 10. -ploading 0.0 -loading -1. -mg_levels_ksp_max_it 2 -snes_monitor_short -ksp_monitor_short -npc_snes_type fas -npc_fas_levels_snes_type ncg -npc_fas_levels_snes_max_it 3 -npc_snes_monitor_short -npc_fas_levels_snes_linesearch_maxlambda 2.0
1041: requires: !single
1043: test:
1044: suffix: 3
1045: args: -da_refine 1 -da_overlap 3 -da_local_subdomains 4 -snes_type aspin -rad 10.0 -young 10. -ploading 0.0 -loading -0.5 -snes_monitor_short -ksp_monitor_short -npc_sub_snes_rtol 1e-2 -ksp_rtol 1e-2 -ksp_max_it 14 -snes_converged_reason -snes_max_linear_solve_fail 100 -npc_sub_ksp_type preonly -npc_sub_pc_type lu
1046: requires: !single
1048: test:
1049: suffix: 4
1050: args: -da_refine 2 -pc_type mg -rad 10.0 -young 10. -ploading -1. -loading -1. -mg_levels_ksp_max_it 2 -snes_monitor_short -ksp_monitor_short -snes_type newtonal -snes_newtonal_step_size 30 -ksp_rtol 1e-4 -info
1051: requires: !single defined(PETSC_USE_INFO)
1052: filter: grep -h -e "KSP Residual norm" -e "SNES Function norm" -e "Number of SNES iterations" -e "mu:"
1054: test:
1055: suffix: 5
1056: args: -da_refine 2 -pc_type mg -rad 10.0 -young 10. -ploading -1. -loading -1. -mg_levels_ksp_max_it 2 -snes_monitor_short -ksp_monitor_short -snes_type newtonal -snes_newtonal_step_size 30 -snes_newtonal_correction_type normal -ksp_rtol 1e-4
1057: requires: !single
1059: test:
1060: suffix: 6
1061: args: -da_refine 2 -pc_type mg -rad 10.0 -young 10. -ploading -0.5 -loading -1 -mg_levels_ksp_max_it 2 -snes_monitor_short -ksp_monitor_short -snes_type newtonal -snes_newtonal_step_size 0.5 -snes_newtonal_max_continuation_steps 3 -snes_newtonal_scale_rhs false -snes_newtonal_lambda_min -0.1 -ksp_rtol 1e-4
1062: requires: !single
1064: test:
1065: nsize: {{1 2}}
1066: suffix: 7
1067: args: -dm_vec_type kokkos -dm_mat_type aijkokkos -da_refine 1 -rad 10.0 -young 10. -ploading -1. -loading -1. -rescale -pc_type mg -mg_levels_ksp_max_it 2 -snes_monitor_short -snes_type newtonal -snes_newtonal_step_size 10 -snes_newtonal_correction_type exact -ksp_rtol 1e-4
1068: requires: kokkos_kernels !complex !single
1069: timeoutfactor: 3
1071: test:
1072: nsize: {{1 2}}
1073: suffix: 8
1074: args: -dm_vec_type cuda -dm_mat_type aijcusparse -da_refine 1 -rad 10.0 -young 10. -ploading -1. -loading -1. -rescale -pc_type mg -mg_levels_ksp_max_it 2 -snes_monitor_short -snes_type newtonal -snes_newtonal_step_size 10 -snes_newtonal_correction_type exact -ksp_rtol 1e-4
1075: requires: cuda !complex !single
1076: timeoutfactor: 3
1078: test:
1079: nsize: {{1 2}}
1080: suffix: 9
1081: args: -da_refine 2 -rad 10.0 -young 10. -ploading -1. -loading -1. -rescale -pc_type mg -mg_levels_ksp_max_it 2 -snes_monitor_short -snes_type newtonal -snes_newtonal_step_size 10 -snes_newtonal_correction_type normal -ksp_rtol 1e-4
1082: requires: !single
1084: TEST*/