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 17 : int main(int argc,char **argv)
48 : {
49 17 : DM dm;
50 17 : MPI_Comm comm;
51 17 : AppCtx user;
52 17 : EPS eps; /* eigenproblem solver context */
53 17 : ST st;
54 17 : EPSType type;
55 17 : Mat A,B,P;
56 17 : Vec v0;
57 17 : PetscContainer container;
58 17 : PetscInt nev,nconv,m,n,M,N;
59 17 : PetscBool nonlin,flg=PETSC_FALSE,update;
60 17 : SNES snes;
61 17 : PetscReal tol,relerr;
62 17 : PetscBool use_shell_matrix=PETSC_FALSE,test_init_sol=PETSC_FALSE,use_custom_norm=PETSC_FALSE,sign_normalization=PETSC_TRUE;
63 :
64 17 : PetscFunctionBeginUser;
65 17 : PetscCall(SlepcInitialize(&argc,&argv,NULL,help));
66 17 : comm = PETSC_COMM_WORLD;
67 : /* Create a quadrilateral mesh on domain (0,1)x(0,1) */
68 17 : PetscCall(CreateSquareMesh(comm,&dm));
69 : /* Setup basis function */
70 17 : PetscCall(SetupDiscretization(dm));
71 17 : PetscCall(BoundaryGlobalIndex(dm,"marker",&user.bdis));
72 : /* Check if we are going to use shell matrices */
73 17 : PetscCall(PetscOptionsGetBool(NULL,NULL,"-use_shell_matrix",&use_shell_matrix,NULL));
74 17 : 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 10 : PetscCall(DMCreateMatrix(dm,&A));
84 10 : PetscCall(MatDuplicate(A,MAT_COPY_VALUES,&B));
85 : }
86 : /* Check whether we should use a custom normalization */
87 17 : 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 17 : 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 17 : PetscCall(PetscObjectComposeFunction((PetscObject)A,"formFunction",FormFunctionA));
95 17 : PetscCall(PetscOptionsGetBool(NULL,NULL,"-form_function_ab",&flg,NULL));
96 17 : if (flg) PetscCall(PetscObjectComposeFunction((PetscObject)A,"formFunctionAB",FormFunctionAB));
97 17 : PetscCall(PetscObjectComposeFunction((PetscObject)A,"formJacobian",FormJacobianA));
98 17 : PetscCall(PetscObjectComposeFunction((PetscObject)B,"formFunction",FormFunctionB));
99 17 : if (use_custom_norm) PetscCall(PetscObjectComposeFunction((PetscObject)B,"formNorm",FormNorm));
100 17 : PetscCall(PetscContainerCreate(comm,&container));
101 17 : PetscCall(PetscContainerSetPointer(container,&user));
102 17 : PetscCall(PetscObjectCompose((PetscObject)A,"formFunctionCtx",(PetscObject)container));
103 17 : PetscCall(PetscObjectCompose((PetscObject)A,"formJacobianCtx",(PetscObject)container));
104 17 : PetscCall(PetscObjectCompose((PetscObject)B,"formFunctionCtx",(PetscObject)container));
105 17 : if (use_custom_norm) PetscCall(PetscObjectCompose((PetscObject)B,"formNormCtx",(PetscObject)container));
106 17 : PetscCall(PetscContainerDestroy(&container));
107 :
108 : /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
109 : Create the eigensolver and set various options
110 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
111 :
112 17 : PetscCall(EPSCreate(comm,&eps));
113 17 : PetscCall(EPSSetOperators(eps,A,B));
114 17 : PetscCall(EPSSetProblemType(eps,EPS_GNHEP));
115 17 : user.eps = eps;
116 : /*
117 : Use nonlinear inverse iteration
118 : */
119 17 : PetscCall(EPSSetType(eps,EPSPOWER));
120 17 : PetscCall(EPSPowerSetNonlinear(eps,PETSC_TRUE));
121 : /* Set the Bx sign normalization (or not) */
122 17 : PetscCall(EPSPowerSetSignNormalization(eps,sign_normalization));
123 : /*
124 : Attach DM to SNES
125 : */
126 17 : PetscCall(EPSPowerGetSNES(eps,&snes));
127 17 : user.snes = snes;
128 17 : PetscCall(SNESSetDM(snes,dm));
129 17 : PetscCall(EPSSetFromOptions(eps));
130 :
131 : /* Set a preconditioning matrix to ST */
132 17 : 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 17 : PetscCall(EPSSolve(eps));
142 :
143 17 : PetscCall(EPSGetConverged(eps,&nconv));
144 17 : PetscCall(PetscOptionsGetBool(NULL,NULL,"-test_init_sol",&test_init_sol,NULL));
145 17 : 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_CURRENT,PETSC_CURRENT,PETSC_CURRENT,PETSC_CURRENT));
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 17 : PetscCall(EPSGetType(eps,&type));
169 17 : PetscCall(EPSGetTolerances(eps,&tol,NULL));
170 17 : PetscCall(EPSPowerGetNonlinear(eps,&nonlin));
171 17 : PetscCall(EPSPowerGetUpdate(eps,&update));
172 22 : PetscCall(PetscPrintf(comm," Solution method: %s%s\n\n",type,nonlin?(update?" (nonlinear with monolithic update)":" (nonlinear)"):""));
173 17 : PetscCall(EPSGetDimensions(eps,&nev,NULL,NULL));
174 17 : PetscCall(PetscPrintf(comm," Number of requested eigenvalues: %" PetscInt_FMT "\n",nev));
175 :
176 : /* print eigenvalue and error */
177 17 : PetscCall(EPSGetConverged(eps,&nconv));
178 17 : if (nconv>0) {
179 17 : PetscScalar k;
180 17 : PetscReal na,nb;
181 17 : Vec a,b,eigen;
182 17 : PetscCall(DMCreateGlobalVector(dm,&a));
183 17 : PetscCall(VecDuplicate(a,&b));
184 17 : PetscCall(VecDuplicate(a,&eigen));
185 17 : PetscCall(EPSGetEigenpair(eps,0,&k,NULL,eigen,NULL));
186 17 : PetscCall(FormFunctionA(snes,eigen,a,&user));
187 17 : PetscCall(FormFunctionB(snes,eigen,b,&user));
188 17 : PetscCall(VecAXPY(a,-k,b));
189 17 : PetscCall(VecNorm(a,NORM_2,&na));
190 17 : if (use_custom_norm) PetscCall(FormNorm(snes,b,&nb,&user));
191 12 : else PetscCall(VecNorm(b,NORM_2,&nb));
192 17 : relerr = na/(nb*PetscAbsScalar(k));
193 17 : 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 17 : PetscCall(VecDestroy(&a));
196 17 : PetscCall(VecDestroy(&b));
197 17 : PetscCall(VecDestroy(&eigen));
198 0 : } else PetscCall(PetscPrintf(comm,"Solver did not converge\n"));
199 :
200 17 : PetscCall(MatDestroy(&A));
201 17 : PetscCall(MatDestroy(&B));
202 17 : if (use_shell_matrix) PetscCall(MatDestroy(&P));
203 17 : PetscCall(DMDestroy(&dm));
204 17 : PetscCall(EPSDestroy(&eps));
205 17 : PetscCall(ISDestroy(&user.bdis));
206 17 : PetscCall(SlepcFinalize());
207 : return 0;
208 : }
209 :
210 : /* <|u|u, v> */
211 234240 : 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 234240 : PetscScalar cof = PetscAbsScalar(u[0]);
217 :
218 234240 : f0[0] = cof*u[0];
219 234240 : }
220 :
221 : /* <|\nabla u| \nabla u, \nabla v> */
222 5744384 : 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 5744384 : PetscInt d;
228 5744384 : PetscScalar cof = 0;
229 17233152 : for (d = 0; d < dim; ++d) cof += u_x[d]*u_x[d];
230 :
231 5744384 : cof = PetscSqrtScalar(cof);
232 :
233 17233152 : for (d = 0; d < dim; ++d) f1[d] = u_x[d]*cof;
234 5744384 : }
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 196352 : 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 196352 : PetscInt d;
252 :
253 589056 : for (d = 0; d < dim; ++d) g3[d*dim+d] = 1.0;
254 196352 : }
255 :
256 17 : PetscErrorCode SetupDiscretization(DM dm)
257 : {
258 17 : PetscFE fe;
259 17 : MPI_Comm comm;
260 :
261 17 : PetscFunctionBeginUser;
262 : /* Create finite element */
263 17 : PetscCall(PetscObjectGetComm((PetscObject)dm,&comm));
264 17 : PetscCall(PetscFECreateDefault(comm,2,1,PETSC_FALSE,NULL,-1,&fe));
265 17 : PetscCall(PetscObjectSetName((PetscObject)fe,"u"));
266 17 : PetscCall(DMSetField(dm,0,NULL,(PetscObject)fe));
267 17 : PetscCall(DMCreateDS(dm));
268 17 : PetscCall(PetscFEDestroy(&fe));
269 17 : PetscFunctionReturn(PETSC_SUCCESS);
270 : }
271 :
272 17 : PetscErrorCode CreateSquareMesh(MPI_Comm comm,DM *dm)
273 : {
274 17 : PetscInt cells[] = {8,8};
275 17 : PetscInt dim = 2;
276 17 : DM pdm;
277 17 : PetscMPIInt size;
278 :
279 17 : PetscFunctionBegin;
280 17 : PetscCall(DMPlexCreateBoxMesh(comm,dim,PETSC_FALSE,cells,NULL,NULL,NULL,PETSC_TRUE,0,PETSC_TRUE,dm));
281 17 : PetscCall(DMSetFromOptions(*dm));
282 17 : PetscCall(DMSetUp(*dm));
283 17 : PetscCallMPI(MPI_Comm_size(comm,&size));
284 17 : if (size > 1) {
285 0 : PetscCall(DMPlexDistribute(*dm,0,NULL,&pdm));
286 0 : PetscCall(DMDestroy(dm));
287 0 : *dm = pdm;
288 : }
289 17 : PetscFunctionReturn(PETSC_SUCCESS);
290 : }
291 :
292 17 : PetscErrorCode BoundaryGlobalIndex(DM dm,const char labelname[],IS *bdis)
293 : {
294 17 : IS bdpoints;
295 17 : PetscInt nindices,*indices,numDof,offset,npoints,i,j;
296 17 : const PetscInt *bdpoints_indices;
297 17 : DMLabel bdmarker;
298 17 : PetscSection gsection;
299 :
300 17 : PetscFunctionBegin;
301 17 : PetscCall(DMGetGlobalSection(dm,&gsection));
302 17 : PetscCall(DMGetLabel(dm,labelname,&bdmarker));
303 17 : PetscCall(DMLabelGetStratumIS(bdmarker,1,&bdpoints));
304 17 : PetscCall(ISGetLocalSize(bdpoints,&npoints));
305 17 : PetscCall(ISGetIndices(bdpoints,&bdpoints_indices));
306 : nindices = 0;
307 1105 : for (i=0;i<npoints;i++) {
308 1088 : PetscCall(PetscSectionGetDof(gsection,bdpoints_indices[i],&numDof));
309 1088 : if (numDof<=0) continue;
310 544 : nindices += numDof;
311 : }
312 17 : PetscCall(PetscCalloc1(nindices,&indices));
313 : nindices = 0;
314 1105 : for (i=0;i<npoints;i++) {
315 1088 : PetscCall(PetscSectionGetDof(gsection,bdpoints_indices[i],&numDof));
316 1088 : if (numDof<=0) continue;
317 544 : PetscCall(PetscSectionGetOffset(gsection,bdpoints_indices[i],&offset));
318 1088 : for (j=0;j<numDof;j++) indices[nindices++] = offset+j;
319 : }
320 17 : PetscCall(ISRestoreIndices(bdpoints,&bdpoints_indices));
321 17 : PetscCall(ISDestroy(&bdpoints));
322 17 : PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)dm),nindices,indices,PETSC_OWN_POINTER,bdis));
323 17 : PetscFunctionReturn(PETSC_SUCCESS);
324 : }
325 :
326 767 : static PetscErrorCode FormJacobian(SNES snes,Vec X,Mat A,Mat B,void *ctx)
327 : {
328 767 : DM dm;
329 767 : Vec Xloc;
330 :
331 767 : PetscFunctionBegin;
332 767 : PetscCall(SNESGetDM(snes,&dm));
333 767 : PetscCall(DMGetLocalVector(dm,&Xloc));
334 767 : PetscCall(VecZeroEntries(Xloc));
335 767 : PetscCall(DMGlobalToLocalBegin(dm,X,INSERT_VALUES,Xloc));
336 767 : PetscCall(DMGlobalToLocalEnd(dm,X,INSERT_VALUES,Xloc));
337 767 : CHKMEMQ;
338 767 : PetscCall(DMPlexSNESComputeJacobianFEM(dm,Xloc,A,B,ctx));
339 767 : if (A!=B) {
340 350 : PetscCall(MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY));
341 350 : PetscCall(MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY));
342 : }
343 767 : CHKMEMQ;
344 767 : PetscCall(DMRestoreLocalVector(dm,&Xloc));
345 767 : PetscFunctionReturn(PETSC_SUCCESS);
346 : }
347 :
348 767 : PetscErrorCode FormJacobianA(SNES snes,Vec X,Mat A,Mat B,void *ctx)
349 : {
350 767 : DM dm;
351 767 : PetscDS prob;
352 767 : PetscWeakForm wf;
353 767 : AppCtx *userctx = (AppCtx *)ctx;
354 :
355 767 : PetscFunctionBegin;
356 767 : PetscCall(MatSetOption(B,MAT_KEEP_NONZERO_PATTERN,PETSC_TRUE));
357 767 : PetscCall(SNESGetDM(snes,&dm));
358 767 : PetscCall(DMGetDS(dm,&prob));
359 767 : PetscCall(PetscDSGetWeakForm(prob, &wf));
360 767 : PetscCall(PetscWeakFormClearIndex(wf, NULL, 0, 0, 0, PETSC_WF_G3, 0));
361 767 : PetscCall(PetscWeakFormSetIndexJacobian(wf, NULL, 0, 0, 0, 0, 0, NULL, 0, NULL, 0, NULL, 0, g3_uu));
362 767 : PetscCall(FormJacobian(snes,X,A,B,ctx));
363 767 : PetscCall(MatZeroRowsIS(B,userctx->bdis,1.0,NULL,NULL));
364 767 : 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 237 : PetscErrorCode FormFunctionAB(SNES snes,Vec x,Vec Ax,Vec Bx,void *ctx)
387 : {
388 237 : 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 237 : PetscCall(FormFunctionA(snes,x,Ax,ctx));
396 237 : PetscCall(FormFunctionB(snes,x,Bx,ctx));
397 237 : PetscFunctionReturn(PETSC_SUCCESS);
398 : }
399 :
400 23354 : static PetscErrorCode FormFunction(SNES snes,Vec X,Vec F,void *ctx)
401 : {
402 23354 : DM dm;
403 23354 : Vec Xloc,Floc;
404 :
405 23354 : PetscFunctionBegin;
406 23354 : PetscCall(SNESGetDM(snes,&dm));
407 23354 : PetscCall(DMGetLocalVector(dm,&Xloc));
408 23354 : PetscCall(DMGetLocalVector(dm,&Floc));
409 23354 : PetscCall(VecZeroEntries(Xloc));
410 23354 : PetscCall(VecZeroEntries(Floc));
411 23354 : PetscCall(DMGlobalToLocalBegin(dm,X,INSERT_VALUES,Xloc));
412 23354 : PetscCall(DMGlobalToLocalEnd(dm,X,INSERT_VALUES,Xloc));
413 23354 : CHKMEMQ;
414 23354 : PetscCall(DMPlexSNESComputeResidualFEM(dm,Xloc,Floc,ctx));
415 23354 : CHKMEMQ;
416 23354 : PetscCall(VecZeroEntries(F));
417 23354 : PetscCall(DMLocalToGlobalBegin(dm,Floc,ADD_VALUES,F));
418 23354 : PetscCall(DMLocalToGlobalEnd(dm,Floc,ADD_VALUES,F));
419 23354 : PetscCall(DMRestoreLocalVector(dm,&Xloc));
420 23354 : PetscCall(DMRestoreLocalVector(dm,&Floc));
421 23354 : PetscFunctionReturn(PETSC_SUCCESS);
422 : }
423 :
424 22439 : PetscErrorCode FormFunctionA(SNES snes,Vec X,Vec F,void *ctx)
425 : {
426 22439 : DM dm;
427 22439 : PetscDS prob;
428 22439 : PetscWeakForm wf;
429 22439 : PetscInt nindices,iStart,iEnd,i;
430 22439 : AppCtx *userctx = (AppCtx *)ctx;
431 22439 : PetscScalar *array,value;
432 22439 : const PetscInt *indices;
433 22439 : PetscInt vecstate;
434 :
435 22439 : PetscFunctionBegin;
436 22439 : PetscCall(SNESGetDM(snes,&dm));
437 22439 : PetscCall(DMGetDS(dm,&prob));
438 : /* hook functions */
439 22439 : PetscCall(PetscDSGetWeakForm(prob, &wf));
440 22439 : PetscCall(PetscWeakFormClearIndex(wf, NULL, 0, 0, 0, PETSC_WF_F0, 0));
441 22439 : PetscCall(PetscWeakFormSetIndexResidual(wf, NULL, 0, 0, 0, 0, NULL, 0, f1_u));
442 22439 : PetscCall(FormFunction(snes,X,F,ctx));
443 : /* Boundary condition */
444 22439 : PetscCall(VecLockGet(X,&vecstate));
445 22439 : if (vecstate>0) PetscCall(VecLockReadPop(X));
446 22439 : PetscCall(VecGetOwnershipRange(X,&iStart,&iEnd));
447 22439 : PetscCall(VecGetArray(X,&array));
448 22439 : PetscCall(ISGetLocalSize(userctx->bdis,&nindices));
449 22439 : PetscCall(ISGetIndices(userctx->bdis,&indices));
450 740487 : for (i=0;i<nindices;i++) {
451 718048 : value = array[indices[i]-iStart] - 0.0;
452 718048 : PetscCall(VecSetValue(F,indices[i],value,INSERT_VALUES));
453 : }
454 22439 : PetscCall(ISRestoreIndices(userctx->bdis,&indices));
455 22439 : PetscCall(VecRestoreArray(X,&array));
456 22439 : if (vecstate>0) PetscCall(VecLockReadPush(X));
457 22439 : PetscCall(VecAssemblyBegin(F));
458 22439 : PetscCall(VecAssemblyEnd(F));
459 22439 : PetscFunctionReturn(PETSC_SUCCESS);
460 : }
461 :
462 293 : PetscErrorCode FormNorm(SNES snes,Vec Bx,PetscReal *norm,void *ctx)
463 : {
464 293 : PetscFunctionBegin;
465 293 : PetscCall(VecNorm(Bx,NORM_2,norm));
466 293 : 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 915 : PetscErrorCode FormFunctionB(SNES snes,Vec X,Vec F,void *ctx)
480 : {
481 915 : DM dm;
482 915 : PetscDS prob;
483 915 : PetscWeakForm wf;
484 915 : PetscInt nindices,iStart,iEnd,i;
485 915 : AppCtx *userctx = (AppCtx *)ctx;
486 915 : PetscScalar value;
487 915 : const PetscInt *indices;
488 :
489 915 : PetscFunctionBegin;
490 915 : PetscCall(SNESGetDM(snes,&dm));
491 915 : PetscCall(DMGetDS(dm,&prob));
492 : /* hook functions */
493 915 : PetscCall(PetscDSGetWeakForm(prob, &wf));
494 915 : PetscCall(PetscWeakFormClearIndex(wf, NULL, 0, 0, 0, PETSC_WF_F1, 0));
495 915 : PetscCall(PetscWeakFormSetIndexResidual(wf, NULL, 0, 0, 0, 0, f0_u, 0, NULL));
496 915 : PetscCall(FormFunction(snes,X,F,ctx));
497 : /* Boundary condition */
498 915 : PetscCall(VecGetOwnershipRange(F,&iStart,&iEnd));
499 915 : PetscCall(ISGetLocalSize(userctx->bdis,&nindices));
500 915 : PetscCall(ISGetIndices(userctx->bdis,&indices));
501 30195 : for (i=0;i<nindices;i++) {
502 29280 : value = 0.0;
503 29280 : PetscCall(VecSetValue(F,indices[i],value,INSERT_VALUES));
504 : }
505 915 : PetscCall(ISRestoreIndices(userctx->bdis,&indices));
506 915 : PetscCall(VecAssemblyBegin(F));
507 915 : PetscCall(VecAssemblyEnd(F));
508 915 : 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*/
|
