Line data Source code
1 : /*
2 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3 : SLEPc - Scalable Library for Eigenvalue Problem Computations
4 : Copyright (c) 2002-, Universitat Politecnica de Valencia, Spain
5 :
6 : This file is part of SLEPc.
7 : SLEPc is distributed under a 2-clause BSD license (see LICENSE).
8 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
9 : */
10 : /*
11 : This is a nonlinear eigenvalue problem. When p=2, it is reduced to a linear Laplace eigenvalue
12 : problem.
13 :
14 : -\nabla\cdot(|\nabla u|^{p-2} \nabla u) = k |u|^{p-2} u in (0,1)x(0,1),
15 :
16 : u = 0 on the entire boundary.
17 :
18 : The code is implemented based on DMPlex using Q1 FEM on a quadrilateral mesh. In this code, we consider p=3.
19 :
20 : Contributed by Fande Kong fdkong.jd@gmail.com
21 : */
22 :
23 : static char help[] = "Nonlinear inverse iteration for A(x)*x=lambda*B(x)*x.\n\n";
24 :
25 : #include <slepceps.h>
26 : #include <petscdmplex.h>
27 : #include <petscds.h>
28 :
29 : PetscErrorCode CreateSquareMesh(MPI_Comm,DM*);
30 : PetscErrorCode SetupDiscretization(DM);
31 : PetscErrorCode FormJacobianA(SNES,Vec,Mat,Mat,void*);
32 : PetscErrorCode FormFunctionA(SNES,Vec,Vec,void*);
33 : PetscErrorCode MatMult_A(Mat A,Vec x,Vec y);
34 : PetscErrorCode FormJacobianB(SNES,Vec,Mat,Mat,void*);
35 : PetscErrorCode FormFunctionB(SNES,Vec,Vec,void*);
36 : PetscErrorCode MatMult_B(Mat A,Vec x,Vec y);
37 : PetscErrorCode FormFunctionAB(SNES,Vec,Vec,Vec,void*);
38 : PetscErrorCode BoundaryGlobalIndex(DM,const char*,IS*);
39 : PetscErrorCode FormNorm(SNES,Vec,PetscReal*,void*);
40 :
41 : typedef struct {
42 : IS bdis; /* global indices for boundary DoFs */
43 : SNES snes;
44 : EPS eps;
45 : } AppCtx;
46 :
47 19 : int main(int argc,char **argv)
48 : {
49 19 : DM dm;
50 19 : MPI_Comm comm;
51 19 : AppCtx user;
52 19 : EPS eps; /* eigenproblem solver context */
53 19 : ST st;
54 19 : EPSType type;
55 19 : Mat A,B,P;
56 19 : Vec v0;
57 19 : PetscContainer container;
58 19 : PetscInt nev,nconv,m,n,M,N;
59 19 : PetscBool nonlin,flg=PETSC_FALSE,update;
60 19 : SNES snes;
61 19 : PetscReal tol,relerr;
62 19 : PetscBool use_shell_matrix=PETSC_FALSE,test_init_sol=PETSC_FALSE,use_custom_norm=PETSC_FALSE,sign_normalization=PETSC_TRUE;
63 :
64 19 : PetscFunctionBeginUser;
65 19 : PetscCall(SlepcInitialize(&argc,&argv,(char*)0,help));
66 19 : comm = PETSC_COMM_WORLD;
67 : /* Create a quadrilateral mesh on domain (0,1)x(0,1) */
68 19 : PetscCall(CreateSquareMesh(comm,&dm));
69 : /* Setup basis function */
70 19 : PetscCall(SetupDiscretization(dm));
71 19 : PetscCall(BoundaryGlobalIndex(dm,"marker",&user.bdis));
72 : /* Check if we are going to use shell matrices */
73 19 : PetscCall(PetscOptionsGetBool(NULL,NULL,"-use_shell_matrix",&use_shell_matrix,NULL));
74 19 : if (use_shell_matrix) {
75 7 : PetscCall(DMCreateMatrix(dm,&P));
76 7 : PetscCall(MatGetLocalSize(P,&m,&n));
77 7 : PetscCall(MatGetSize(P,&M,&N));
78 7 : PetscCall(MatCreateShell(comm,m,n,M,N,&user,&A));
79 7 : PetscCall(MatShellSetOperation(A,MATOP_MULT,(void(*)(void))MatMult_A));
80 7 : PetscCall(MatCreateShell(comm,m,n,M,N,&user,&B));
81 7 : PetscCall(MatShellSetOperation(B,MATOP_MULT,(void(*)(void))MatMult_B));
82 : } else {
83 12 : PetscCall(DMCreateMatrix(dm,&A));
84 12 : PetscCall(MatDuplicate(A,MAT_COPY_VALUES,&B));
85 : }
86 : /* Check whether we should use a custom normalization */
87 19 : PetscCall(PetscOptionsGetBool(NULL,NULL,"-use_custom_norm",&use_custom_norm,NULL));
88 : /* Check whether we should normalize Bx by the sign of its first nonzero element */
89 19 : PetscCall(PetscOptionsGetBool(NULL,NULL,"-sign_normalization",&sign_normalization,NULL));
90 :
91 : /*
92 : Compose callback functions and context that will be needed by the solver
93 : */
94 19 : PetscCall(PetscObjectComposeFunction((PetscObject)A,"formFunction",FormFunctionA));
95 19 : PetscCall(PetscOptionsGetBool(NULL,NULL,"-form_function_ab",&flg,NULL));
96 19 : if (flg) PetscCall(PetscObjectComposeFunction((PetscObject)A,"formFunctionAB",FormFunctionAB));
97 19 : PetscCall(PetscObjectComposeFunction((PetscObject)A,"formJacobian",FormJacobianA));
98 19 : PetscCall(PetscObjectComposeFunction((PetscObject)B,"formFunction",FormFunctionB));
99 19 : if (use_custom_norm) PetscCall(PetscObjectComposeFunction((PetscObject)B,"formNorm",FormNorm));
100 19 : PetscCall(PetscContainerCreate(comm,&container));
101 19 : PetscCall(PetscContainerSetPointer(container,&user));
102 19 : PetscCall(PetscObjectCompose((PetscObject)A,"formFunctionCtx",(PetscObject)container));
103 19 : PetscCall(PetscObjectCompose((PetscObject)A,"formJacobianCtx",(PetscObject)container));
104 19 : PetscCall(PetscObjectCompose((PetscObject)B,"formFunctionCtx",(PetscObject)container));
105 19 : if (use_custom_norm) PetscCall(PetscObjectCompose((PetscObject)B,"formNormCtx",(PetscObject)container));
106 19 : PetscCall(PetscContainerDestroy(&container));
107 :
108 : /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
109 : Create the eigensolver and set various options
110 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
111 :
112 19 : PetscCall(EPSCreate(comm,&eps));
113 19 : PetscCall(EPSSetOperators(eps,A,B));
114 19 : PetscCall(EPSSetProblemType(eps,EPS_GNHEP));
115 19 : user.eps = eps;
116 : /*
117 : Use nonlinear inverse iteration
118 : */
119 19 : PetscCall(EPSSetType(eps,EPSPOWER));
120 19 : PetscCall(EPSPowerSetNonlinear(eps,PETSC_TRUE));
121 : /* Set the Bx sign normalization (or not) */
122 19 : PetscCall(EPSPowerSetSignNormalization(eps,sign_normalization));
123 : /*
124 : Attach DM to SNES
125 : */
126 19 : PetscCall(EPSPowerGetSNES(eps,&snes));
127 19 : user.snes = snes;
128 19 : PetscCall(SNESSetDM(snes,dm));
129 19 : PetscCall(EPSSetFromOptions(eps));
130 :
131 : /* Set a preconditioning matrix to ST */
132 19 : if (use_shell_matrix) {
133 7 : PetscCall(EPSGetST(eps,&st));
134 7 : PetscCall(STSetPreconditionerMat(st,P));
135 : }
136 :
137 : /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
138 : Solve the eigensystem
139 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
140 :
141 19 : PetscCall(EPSSolve(eps));
142 :
143 19 : PetscCall(EPSGetConverged(eps,&nconv));
144 19 : PetscCall(PetscOptionsGetBool(NULL,NULL,"-test_init_sol",&test_init_sol,NULL));
145 19 : if (nconv && test_init_sol) {
146 2 : PetscScalar k;
147 2 : PetscReal norm0;
148 2 : PetscInt nits;
149 :
150 2 : PetscCall(MatCreateVecs(A,&v0,NULL));
151 2 : PetscCall(EPSGetEigenpair(eps,0,&k,NULL,v0,NULL));
152 2 : PetscCall(EPSSetInitialSpace(eps,1,&v0));
153 2 : PetscCall(VecDestroy(&v0));
154 : /* Norm of the previous residual */
155 2 : PetscCall(SNESGetFunctionNorm(snes,&norm0));
156 : /* Make the tolerance smaller than the last residual
157 : SNES will converge right away if the initial is setup correctly */
158 2 : PetscCall(SNESSetTolerances(snes,norm0*1.2,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT));
159 2 : PetscCall(EPSSolve(eps));
160 : /* Number of Newton iterations supposes to be zero */
161 2 : PetscCall(SNESGetIterationNumber(snes,&nits));
162 2 : if (nits) PetscCall(PetscPrintf(comm," Number of Newton iterations %" PetscInt_FMT " should be zero \n",nits));
163 : }
164 :
165 : /*
166 : Optional: Get some information from the solver and display it
167 : */
168 19 : PetscCall(EPSGetType(eps,&type));
169 19 : PetscCall(EPSGetTolerances(eps,&tol,NULL));
170 19 : PetscCall(EPSPowerGetNonlinear(eps,&nonlin));
171 19 : PetscCall(EPSPowerGetUpdate(eps,&update));
172 24 : PetscCall(PetscPrintf(comm," Solution method: %s%s\n\n",type,nonlin?(update?" (nonlinear with monolithic update)":" (nonlinear)"):""));
173 19 : PetscCall(EPSGetDimensions(eps,&nev,NULL,NULL));
174 19 : PetscCall(PetscPrintf(comm," Number of requested eigenvalues: %" PetscInt_FMT "\n",nev));
175 :
176 : /* print eigenvalue and error */
177 19 : PetscCall(EPSGetConverged(eps,&nconv));
178 19 : if (nconv>0) {
179 19 : PetscScalar k;
180 19 : PetscReal na,nb;
181 19 : Vec a,b,eigen;
182 19 : PetscCall(DMCreateGlobalVector(dm,&a));
183 19 : PetscCall(VecDuplicate(a,&b));
184 19 : PetscCall(VecDuplicate(a,&eigen));
185 19 : PetscCall(EPSGetEigenpair(eps,0,&k,NULL,eigen,NULL));
186 19 : PetscCall(FormFunctionA(snes,eigen,a,&user));
187 19 : PetscCall(FormFunctionB(snes,eigen,b,&user));
188 19 : PetscCall(VecAXPY(a,-k,b));
189 19 : PetscCall(VecNorm(a,NORM_2,&na));
190 19 : if (use_custom_norm) PetscCall(FormNorm(snes,b,&nb,&user));
191 13 : else PetscCall(VecNorm(b,NORM_2,&nb));
192 19 : relerr = na/(nb*PetscAbsScalar(k));
193 19 : if (relerr<10*tol) PetscCall(PetscPrintf(comm,"k: %g, relative error below tol\n",(double)PetscRealPart(k)));
194 0 : else PetscCall(PetscPrintf(comm,"k: %g, relative error: %g\n",(double)PetscRealPart(k),(double)relerr));
195 19 : PetscCall(VecDestroy(&a));
196 19 : PetscCall(VecDestroy(&b));
197 19 : PetscCall(VecDestroy(&eigen));
198 0 : } else PetscCall(PetscPrintf(comm,"Solver did not converge\n"));
199 :
200 19 : PetscCall(MatDestroy(&A));
201 19 : PetscCall(MatDestroy(&B));
202 19 : if (use_shell_matrix) PetscCall(MatDestroy(&P));
203 19 : PetscCall(DMDestroy(&dm));
204 19 : PetscCall(EPSDestroy(&eps));
205 19 : PetscCall(ISDestroy(&user.bdis));
206 19 : PetscCall(SlepcFinalize());
207 : return 0;
208 : }
209 :
210 : /* <|u|u, v> */
211 172544 : static void f0_u(PetscInt dim, PetscInt Nf, PetscInt NfAux,
212 : const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
213 : const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
214 : PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[])
215 : {
216 172544 : PetscScalar cof = PetscAbsScalar(u[0]);
217 :
218 172544 : f0[0] = cof*u[0];
219 172544 : }
220 :
221 : /* <|\nabla u| \nabla u, \nabla v> */
222 1378304 : static void f1_u(PetscInt dim, PetscInt Nf, PetscInt NfAux,
223 : const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
224 : const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
225 : PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f1[])
226 : {
227 1378304 : PetscInt d;
228 1378304 : PetscScalar cof = 0;
229 4134912 : for (d = 0; d < dim; ++d) cof += u_x[d]*u_x[d];
230 :
231 1378304 : cof = PetscSqrtScalar(cof);
232 :
233 4134912 : for (d = 0; d < dim; ++d) f1[d] = u_x[d]*cof;
234 1378304 : }
235 :
236 : /* approximate Jacobian for <|u|u, v> */
237 0 : static void g0_uu(PetscInt dim, PetscInt Nf, PetscInt NfAux,
238 : const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
239 : const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
240 : PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g0[])
241 : {
242 0 : g0[0] = 1.0*PetscAbsScalar(u[0]);
243 0 : }
244 :
245 : /* approximate Jacobian for <|\nabla u| \nabla u, \nabla v> */
246 219392 : static void g3_uu(PetscInt dim, PetscInt Nf, PetscInt NfAux,
247 : const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
248 : const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
249 : PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g3[])
250 : {
251 219392 : PetscInt d;
252 :
253 658176 : for (d = 0; d < dim; ++d) g3[d*dim+d] = 1.0;
254 219392 : }
255 :
256 19 : PetscErrorCode SetupDiscretization(DM dm)
257 : {
258 19 : PetscFE fe;
259 19 : MPI_Comm comm;
260 :
261 19 : PetscFunctionBeginUser;
262 : /* Create finite element */
263 19 : PetscCall(PetscObjectGetComm((PetscObject)dm,&comm));
264 19 : PetscCall(PetscFECreateDefault(comm,2,1,PETSC_FALSE,NULL,-1,&fe));
265 19 : PetscCall(PetscObjectSetName((PetscObject)fe,"u"));
266 19 : PetscCall(DMSetField(dm,0,NULL,(PetscObject)fe));
267 19 : PetscCall(DMCreateDS(dm));
268 19 : PetscCall(PetscFEDestroy(&fe));
269 19 : PetscFunctionReturn(PETSC_SUCCESS);
270 : }
271 :
272 19 : PetscErrorCode CreateSquareMesh(MPI_Comm comm,DM *dm)
273 : {
274 19 : PetscInt cells[] = {8,8};
275 19 : PetscInt dim = 2;
276 19 : DM pdm;
277 19 : PetscMPIInt size;
278 :
279 19 : PetscFunctionBegin;
280 19 : PetscCall(DMPlexCreateBoxMesh(comm,dim,PETSC_FALSE,cells,NULL,NULL,NULL,PETSC_TRUE,dm));
281 19 : PetscCall(DMSetFromOptions(*dm));
282 19 : PetscCall(DMSetUp(*dm));
283 19 : PetscCallMPI(MPI_Comm_size(comm,&size));
284 19 : if (size > 1) {
285 0 : PetscCall(DMPlexDistribute(*dm,0,NULL,&pdm));
286 0 : PetscCall(DMDestroy(dm));
287 0 : *dm = pdm;
288 : }
289 19 : PetscFunctionReturn(PETSC_SUCCESS);
290 : }
291 :
292 19 : PetscErrorCode BoundaryGlobalIndex(DM dm,const char labelname[],IS *bdis)
293 : {
294 19 : IS bdpoints;
295 19 : PetscInt nindices,*indices,numDof,offset,npoints,i,j;
296 19 : const PetscInt *bdpoints_indices;
297 19 : DMLabel bdmarker;
298 19 : PetscSection gsection;
299 :
300 19 : PetscFunctionBegin;
301 19 : PetscCall(DMGetGlobalSection(dm,&gsection));
302 19 : PetscCall(DMGetLabel(dm,labelname,&bdmarker));
303 19 : PetscCall(DMLabelGetStratumIS(bdmarker,1,&bdpoints));
304 19 : PetscCall(ISGetLocalSize(bdpoints,&npoints));
305 19 : PetscCall(ISGetIndices(bdpoints,&bdpoints_indices));
306 : nindices = 0;
307 1235 : for (i=0;i<npoints;i++) {
308 1216 : PetscCall(PetscSectionGetDof(gsection,bdpoints_indices[i],&numDof));
309 1216 : if (numDof<=0) continue;
310 608 : nindices += numDof;
311 : }
312 19 : PetscCall(PetscCalloc1(nindices,&indices));
313 : nindices = 0;
314 1235 : for (i=0;i<npoints;i++) {
315 1216 : PetscCall(PetscSectionGetDof(gsection,bdpoints_indices[i],&numDof));
316 1216 : if (numDof<=0) continue;
317 608 : PetscCall(PetscSectionGetOffset(gsection,bdpoints_indices[i],&offset));
318 1216 : for (j=0;j<numDof;j++) indices[nindices++] = offset+j;
319 : }
320 19 : PetscCall(ISRestoreIndices(bdpoints,&bdpoints_indices));
321 19 : PetscCall(ISDestroy(&bdpoints));
322 19 : PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)dm),nindices,indices,PETSC_OWN_POINTER,bdis));
323 19 : PetscFunctionReturn(PETSC_SUCCESS);
324 : }
325 :
326 857 : static PetscErrorCode FormJacobian(SNES snes,Vec X,Mat A,Mat B,void *ctx)
327 : {
328 857 : DM dm;
329 857 : Vec Xloc;
330 :
331 857 : PetscFunctionBegin;
332 857 : PetscCall(SNESGetDM(snes,&dm));
333 857 : PetscCall(DMGetLocalVector(dm,&Xloc));
334 857 : PetscCall(VecZeroEntries(Xloc));
335 857 : PetscCall(DMGlobalToLocalBegin(dm,X,INSERT_VALUES,Xloc));
336 857 : PetscCall(DMGlobalToLocalEnd(dm,X,INSERT_VALUES,Xloc));
337 857 : CHKMEMQ;
338 857 : PetscCall(DMPlexSNESComputeJacobianFEM(dm,Xloc,A,B,ctx));
339 857 : if (A!=B) {
340 456 : PetscCall(MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY));
341 456 : PetscCall(MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY));
342 : }
343 857 : CHKMEMQ;
344 857 : PetscCall(DMRestoreLocalVector(dm,&Xloc));
345 857 : PetscFunctionReturn(PETSC_SUCCESS);
346 : }
347 :
348 857 : PetscErrorCode FormJacobianA(SNES snes,Vec X,Mat A,Mat B,void *ctx)
349 : {
350 857 : DM dm;
351 857 : PetscDS prob;
352 857 : PetscWeakForm wf;
353 857 : AppCtx *userctx = (AppCtx *)ctx;
354 :
355 857 : PetscFunctionBegin;
356 857 : PetscCall(MatSetOption(B,MAT_KEEP_NONZERO_PATTERN,PETSC_TRUE));
357 857 : PetscCall(SNESGetDM(snes,&dm));
358 857 : PetscCall(DMGetDS(dm,&prob));
359 857 : PetscCall(PetscDSGetWeakForm(prob, &wf));
360 857 : PetscCall(PetscWeakFormClearIndex(wf, NULL, 0, 0, 0, PETSC_WF_G3, 0));
361 857 : PetscCall(PetscWeakFormSetIndexJacobian(wf, NULL, 0, 0, 0, 0, 0, NULL, 0, NULL, 0, NULL, 0, g3_uu));
362 857 : PetscCall(FormJacobian(snes,X,A,B,ctx));
363 857 : PetscCall(MatZeroRowsIS(B,userctx->bdis,1.0,NULL,NULL));
364 857 : PetscFunctionReturn(PETSC_SUCCESS);
365 : }
366 :
367 0 : PetscErrorCode FormJacobianB(SNES snes,Vec X,Mat A,Mat B,void *ctx)
368 : {
369 0 : DM dm;
370 0 : PetscDS prob;
371 0 : PetscWeakForm wf;
372 0 : AppCtx *userctx = (AppCtx *)ctx;
373 :
374 0 : PetscFunctionBegin;
375 0 : PetscCall(MatSetOption(B,MAT_KEEP_NONZERO_PATTERN,PETSC_TRUE));
376 0 : PetscCall(SNESGetDM(snes,&dm));
377 0 : PetscCall(DMGetDS(dm,&prob));
378 0 : PetscCall(PetscDSGetWeakForm(prob, &wf));
379 0 : PetscCall(PetscWeakFormClearIndex(wf, NULL, 0, 0, 0, PETSC_WF_G3, 0));
380 0 : PetscCall(PetscWeakFormSetIndexJacobian(wf, NULL, 0, 0, 0, 0, 0, g0_uu, 0, NULL, 0, NULL, 0, NULL));
381 0 : PetscCall(FormJacobian(snes,X,A,B,ctx));
382 0 : PetscCall(MatZeroRowsIS(B,userctx->bdis,0.0,NULL,NULL));
383 0 : PetscFunctionReturn(PETSC_SUCCESS);
384 : }
385 :
386 222 : PetscErrorCode FormFunctionAB(SNES snes,Vec x,Vec Ax,Vec Bx,void *ctx)
387 : {
388 222 : PetscFunctionBegin;
389 : /*
390 : * In real applications, users should have a generic formFunctionAB which
391 : * forms Ax and Bx simultaneously for an more efficient calculation.
392 : * In this example, we just call FormFunctionA+FormFunctionB to mimic how
393 : * to use FormFunctionAB
394 : */
395 222 : PetscCall(FormFunctionA(snes,x,Ax,ctx));
396 222 : PetscCall(FormFunctionB(snes,x,Bx,ctx));
397 222 : PetscFunctionReturn(PETSC_SUCCESS);
398 : }
399 :
400 6058 : static PetscErrorCode FormFunction(SNES snes,Vec X,Vec F,void *ctx)
401 : {
402 6058 : DM dm;
403 6058 : Vec Xloc,Floc;
404 :
405 6058 : PetscFunctionBegin;
406 6058 : PetscCall(SNESGetDM(snes,&dm));
407 6058 : PetscCall(DMGetLocalVector(dm,&Xloc));
408 6058 : PetscCall(DMGetLocalVector(dm,&Floc));
409 6058 : PetscCall(VecZeroEntries(Xloc));
410 6058 : PetscCall(VecZeroEntries(Floc));
411 6058 : PetscCall(DMGlobalToLocalBegin(dm,X,INSERT_VALUES,Xloc));
412 6058 : PetscCall(DMGlobalToLocalEnd(dm,X,INSERT_VALUES,Xloc));
413 6058 : CHKMEMQ;
414 6058 : PetscCall(DMPlexSNESComputeResidualFEM(dm,Xloc,Floc,ctx));
415 6058 : CHKMEMQ;
416 6058 : PetscCall(VecZeroEntries(F));
417 6058 : PetscCall(DMLocalToGlobalBegin(dm,Floc,ADD_VALUES,F));
418 6058 : PetscCall(DMLocalToGlobalEnd(dm,Floc,ADD_VALUES,F));
419 6058 : PetscCall(DMRestoreLocalVector(dm,&Xloc));
420 6058 : PetscCall(DMRestoreLocalVector(dm,&Floc));
421 6058 : PetscFunctionReturn(PETSC_SUCCESS);
422 : }
423 :
424 5384 : PetscErrorCode FormFunctionA(SNES snes,Vec X,Vec F,void *ctx)
425 : {
426 5384 : DM dm;
427 5384 : PetscDS prob;
428 5384 : PetscWeakForm wf;
429 5384 : PetscInt nindices,iStart,iEnd,i;
430 5384 : AppCtx *userctx = (AppCtx *)ctx;
431 5384 : PetscScalar *array,value;
432 5384 : const PetscInt *indices;
433 5384 : PetscInt vecstate;
434 :
435 5384 : PetscFunctionBegin;
436 5384 : PetscCall(SNESGetDM(snes,&dm));
437 5384 : PetscCall(DMGetDS(dm,&prob));
438 : /* hook functions */
439 5384 : PetscCall(PetscDSGetWeakForm(prob, &wf));
440 5384 : PetscCall(PetscWeakFormClearIndex(wf, NULL, 0, 0, 0, PETSC_WF_F0, 0));
441 5384 : PetscCall(PetscWeakFormSetIndexResidual(wf, NULL, 0, 0, 0, 0, NULL, 0, f1_u));
442 5384 : PetscCall(FormFunction(snes,X,F,ctx));
443 : /* Boundary condition */
444 5384 : PetscCall(VecLockGet(X,&vecstate));
445 5384 : if (vecstate>0) PetscCall(VecLockReadPop(X));
446 5384 : PetscCall(VecGetOwnershipRange(X,&iStart,&iEnd));
447 5384 : PetscCall(VecGetArray(X,&array));
448 5384 : PetscCall(ISGetLocalSize(userctx->bdis,&nindices));
449 5384 : PetscCall(ISGetIndices(userctx->bdis,&indices));
450 177672 : for (i=0;i<nindices;i++) {
451 172288 : value = array[indices[i]-iStart] - 0.0;
452 172288 : PetscCall(VecSetValue(F,indices[i],value,INSERT_VALUES));
453 : }
454 5384 : PetscCall(ISRestoreIndices(userctx->bdis,&indices));
455 5384 : PetscCall(VecRestoreArray(X,&array));
456 5384 : if (vecstate>0) PetscCall(VecLockReadPush(X));
457 5384 : PetscCall(VecAssemblyBegin(F));
458 5384 : PetscCall(VecAssemblyEnd(F));
459 5384 : PetscFunctionReturn(PETSC_SUCCESS);
460 : }
461 :
462 214 : PetscErrorCode FormNorm(SNES snes,Vec Bx,PetscReal *norm,void *ctx)
463 : {
464 214 : PetscFunctionBegin;
465 214 : PetscCall(VecNorm(Bx,NORM_2,norm));
466 214 : PetscFunctionReturn(PETSC_SUCCESS);
467 : }
468 :
469 0 : PetscErrorCode MatMult_A(Mat A,Vec x,Vec y)
470 : {
471 0 : AppCtx *userctx;
472 :
473 0 : PetscFunctionBegin;
474 0 : PetscCall(MatShellGetContext(A,&userctx));
475 0 : PetscCall(FormFunctionA(userctx->snes,x,y,userctx));
476 0 : PetscFunctionReturn(PETSC_SUCCESS);
477 : }
478 :
479 674 : PetscErrorCode FormFunctionB(SNES snes,Vec X,Vec F,void *ctx)
480 : {
481 674 : DM dm;
482 674 : PetscDS prob;
483 674 : PetscWeakForm wf;
484 674 : PetscInt nindices,iStart,iEnd,i;
485 674 : AppCtx *userctx = (AppCtx *)ctx;
486 674 : PetscScalar value;
487 674 : const PetscInt *indices;
488 :
489 674 : PetscFunctionBegin;
490 674 : PetscCall(SNESGetDM(snes,&dm));
491 674 : PetscCall(DMGetDS(dm,&prob));
492 : /* hook functions */
493 674 : PetscCall(PetscDSGetWeakForm(prob, &wf));
494 674 : PetscCall(PetscWeakFormClearIndex(wf, NULL, 0, 0, 0, PETSC_WF_F1, 0));
495 674 : PetscCall(PetscWeakFormSetIndexResidual(wf, NULL, 0, 0, 0, 0, f0_u, 0, NULL));
496 674 : PetscCall(FormFunction(snes,X,F,ctx));
497 : /* Boundary condition */
498 674 : PetscCall(VecGetOwnershipRange(F,&iStart,&iEnd));
499 674 : PetscCall(ISGetLocalSize(userctx->bdis,&nindices));
500 674 : PetscCall(ISGetIndices(userctx->bdis,&indices));
501 22242 : for (i=0;i<nindices;i++) {
502 21568 : value = 0.0;
503 21568 : PetscCall(VecSetValue(F,indices[i],value,INSERT_VALUES));
504 : }
505 674 : PetscCall(ISRestoreIndices(userctx->bdis,&indices));
506 674 : PetscCall(VecAssemblyBegin(F));
507 674 : PetscCall(VecAssemblyEnd(F));
508 674 : PetscFunctionReturn(PETSC_SUCCESS);
509 : }
510 :
511 0 : PetscErrorCode MatMult_B(Mat B,Vec x,Vec y)
512 : {
513 0 : AppCtx *userctx;
514 :
515 0 : PetscFunctionBegin;
516 0 : PetscCall(MatShellGetContext(B,&userctx));
517 0 : PetscCall(FormFunctionB(userctx->snes,x,y,userctx));
518 0 : PetscFunctionReturn(PETSC_SUCCESS);
519 : }
520 :
521 : /*TEST
522 :
523 : testset:
524 : requires: double
525 : args: -petscspace_degree 1 -petscspace_poly_tensor -checkfunctionlist 0
526 : output_file: output/ex34_1.out
527 : test:
528 : suffix: 1
529 : test:
530 : suffix: 2
531 : args: -eps_power_update -form_function_ab {{0 1}}
532 : filter: sed -e "s/ with monolithic update//"
533 : test:
534 : suffix: 3
535 : args: -use_shell_matrix -eps_power_snes_mf_operator 1
536 : test:
537 : suffix: 4
538 : args: -use_shell_matrix -eps_power_update -init_eps_power_snes_mf_operator 1 -eps_power_snes_mf_operator 1 -form_function_ab {{0 1}}
539 : filter: sed -e "s/ with monolithic update//"
540 : test:
541 : suffix: 5
542 : args: -use_shell_matrix -eps_power_update -init_eps_power_snes_mf_operator 1 -eps_power_snes_mf_operator 1 -form_function_ab {{0 1}} -test_init_sol 1
543 : filter: sed -e "s/ with monolithic update//"
544 :
545 : test:
546 : suffix: 6
547 : args: -use_shell_matrix -eps_power_update -init_eps_power_snes_mf_operator 1 -eps_power_snes_mf_operator 1 -form_function_ab {{0 1}} -eps_monitor_all
548 : output_file: output/ex34_6.out
549 : filter: sed -e "s/\([+-].*i\)//g" -e "1,3s/[0-9]//g" -e "/[45] EPS/d"
550 : test:
551 : suffix: 7
552 : args: -use_custom_norm -sign_normalization 1 -eps_power_snes_mf_operator 1
553 : test:
554 : suffix: 8
555 : args: -use_custom_norm -sign_normalization 1 -eps_power_update -form_function_ab {{0 1}} -eps_power_snes_mf_operator 1 -init_eps_power_snes_mf_operator 1
556 : filter: sed -e "s/ with monolithic update//"
557 : test:
558 : suffix: 9
559 : requires: !complex
560 : args: -use_custom_norm {{0 1}} -sign_normalization 0 -eps_power_snes_mf_operator 1
561 : test:
562 : suffix: 10
563 : requires: !complex
564 : args: -use_custom_norm {{0 1}} -sign_normalization 0 -eps_power_update -form_function_ab {{0 1}} -eps_power_snes_mf_operator 1 -init_eps_power_snes_mf_operator 1
565 : filter: sed -e "s/ with monolithic update//"
566 : test:
567 : suffix: 11
568 : requires: complex
569 : args: -use_custom_norm {{0 1}} -sign_normalization 0 -eps_power_snes_type nrichardson -eps_power_snes_atol 1e-12
570 : test:
571 : suffix: 12
572 : requires: complex
573 : args: -use_custom_norm {{0 1}} -sign_normalization 0 -eps_power_update -init_eps_power_snes_type nrichardson -init_eps_max_it 2 -eps_power_snes_mf_operator 1
574 : filter: sed -e "s/ with monolithic update//"
575 : TEST*/
|
