Actual source code: ex60.c
1: /*
2: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3: SLEPc - Scalable Library for Eigenvalue Problem Computations
4: Copyright (c) 2002-, Universitat Politecnica de Valencia, Spain
6: This file is part of SLEPc.
7: SLEPc is distributed under a 2-clause BSD license (see LICENSE).
8: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
9: */
11: static char help[] = "2-D discrete Gross-Pitaevskii model solved as a Nonlinear Eigenvalue Problem with eigenvector dependence (NEPv).\n\n"
12: "Implemented by embedding an EPS inside a nonlinear SNES loop.\n\n"
13: "The Gross-Pitaevskii (GP) equation models the ground state of a rotating Bose-Einstein Condensate (BEC).\n"
14: "The goal is to self-consistently find the macroscopic wavefunction (eigenvector) and its corresponding\n"
15: "chemical potential or ground-state energy (eigenvalue) under a magnetic trap.\n\n"
16: "The nonlinear Hamiltonian matrix is defined as H(x) = A0 + beta * Diag(rho(x)), where:\n"
17: " A0 = base linear Hamiltonian containing 2-D kinetic energy, a harmonic trap potential, and angular momentum rotation.\n"
18: " rho = condensate density vector computed from the ground-state wavefunction x as rho = |x|^2.\n\n"
19: "The command line options are:\n"
20: " -n <n>, where <n> = number of grid points per dimension.\n"
21: " -beta <beta>, where <beta> = real scaling parameter controlling the nonlinearity and strength of the atomic interactions.\n"
22: " -symm <bool>: trap potential profile (true = symmetric, false = asymmetric).\n"
23: " -verbose, to print the converged Hamiltonian diagonal and its sum.\n\n";
25: #include <slepceps.h>
27: /*
28: Context structure to store the necessary data regarding the discrete Gross-Pitaevskii (DGP) problem
29: as well as the objects needed for the nonlinear iteration loop.
30: */
31: typedef struct {
32: // DGP context
34: // Inputs
35: PetscInt N; // Number of points per dimension
36: PetscReal beta; // Parameter controlling the nonlinearity
37: PetscBool symm; // True: Symmetric potential, False: Asymmetric potential
39: // Constants
40: PetscReal L; // Length of the semi-domain [-L, L]
41: PetscReal h; // Step size
42: PetscReal omega; // Angular frequency
44: // Working vectors and matrices
45: Vec rho; // Vector to store the condensate density
46: Vec z; // Intermediate vector to store the result z = beta * rho
47: Mat A0; // Base linear Hamiltonian
49: // Nonlinear iteration context
50: EPS eps; // Solver for each H
51: Mat H; // Hamiltonian
52: Vec x; // Eigenvector
53: } DGPContext;
55: /*
56: Initialize the context for the 2D discrete Gross-Pitaevskii problem.
58: This function allocates the necessary memory and builds the base linear
59: Hamiltonian A0 as
61: A0 = h^2 * ( -0.5*L + V - omega*1i*L_z )
63: where
64: L is the discrete 2-D Laplacian,
65: L_z is the discrete 2-D angular momentum operator (without scaling), and
66: V is a diagonal matrix with the potential for each coordinate (x,y):
67: (x.^2+y.^2)/2 for the symmetric case, (x.^2+100*y.^2)/2 for the non-symmetric case.
69: Arguments:
70: N - Grid size per dimension
71: beta - Parameter controlling the repulsive interaction of the condensate
72: symm - Defines the shape of the magnetic trap (True = symmetric, False = asymmetric)
73: ctx_out - Pointer to the memory address where the created context will be stored
74: */
75: PetscErrorCode DGPCreateContext(MPI_Comm comm,PetscInt N,PetscReal beta,PetscBool symm,DGPContext **ctx_out)
76: {
77: DGPContext *ctx;
79: PetscInt Istart,Iend,i,II,j;
80: PetscReal h2;
81: PetscReal x,y,V;
82: PetscScalar alpha;
84: PetscFunctionBeginUser;
85: PetscCall(PetscNew(&ctx));
87: // Inputs
88: ctx->N=N;
89: ctx->beta=beta;
90: ctx->symm=symm;
92: // Constants
93: ctx->L=1.0;
94: ctx->h=2.0*ctx->L/(ctx->N+1.0);
95: ctx->omega=0.85;
96: h2=ctx->h*ctx->h;
97: // Scaling factor for h^2 * (-omega * 1i * L_z) with L_z = (x or y)/(2h)
98: alpha = -h2*ctx->omega/(2.0*ctx->h)*PETSC_i;
100: // 1. Create A0
101: PetscCall(MatCreate(comm,&ctx->A0));
102: PetscCall(MatSetSizes(ctx->A0,PETSC_DECIDE,PETSC_DECIDE,N*N,N*N));
103: PetscCall(MatSetFromOptions(ctx->A0));
105: // 2. Fill A0
106: PetscCall(MatGetOwnershipRange(ctx->A0,&Istart,&Iend));
108: for (II=Istart;II<Iend;II++) {
109: i=II/ctx->N; j=II%ctx->N;
111: // Potential for mesh point with coordinates (x,y)
112: x=-ctx->L+(j+1)*ctx->h;
113: y=-ctx->L+(i+1)*ctx->h;
114: if (ctx->symm) V=0.5*(x*x+y*y);
115: else V=0.5*(x*x+100.0*y*y);
117: // Diagonal: h^2*(-0.5*L + V - omega*1i*L_z) with L_diag = -4/h^2 and L_z_diag = 0
118: PetscCall(MatSetValue(ctx->A0,II,II,0.5*4.0+h2*V,INSERT_VALUES));
120: // Off-diagonal: h^2*(-0.5*L + V - omega*1i*L_z) with L_offdiag = 1/h^2, V_offdiag = 0, and L_z_offdiag = +/- (x or y)/(2h)
121: if (i>0) PetscCall(MatSetValue(ctx->A0,II,II-ctx->N,-0.5+alpha*x,INSERT_VALUES));
122: if (i<ctx->N-1) PetscCall(MatSetValue(ctx->A0,II,II+ctx->N,-0.5-alpha*x,INSERT_VALUES));
123: if (j>0) PetscCall(MatSetValue(ctx->A0,II,II-1,-0.5-alpha*y,INSERT_VALUES));
124: if (j<ctx->N-1) PetscCall(MatSetValue(ctx->A0,II,II+1,-0.5+alpha*y,INSERT_VALUES));
125: }
127: PetscCall(MatAssemblyBegin(ctx->A0,MAT_FINAL_ASSEMBLY));
128: PetscCall(MatAssemblyEnd(ctx->A0,MAT_FINAL_ASSEMBLY));
130: // 3. Create vectors to store the condensate density
131: PetscCall(MatCreateVecs(ctx->A0,&ctx->rho,&ctx->z));
133: *ctx_out=ctx;
134: PetscFunctionReturn(PETSC_SUCCESS);
135: }
137: /*
138: Configure the SLEPc objects necessary for the nonlinear iteration.
140: This function prepares the Hamiltonian matrix by cloning the structure of A0,
141: initializes the vector to store the condensate state, and configures the eigenvalue
142: solver (EPS) by defining the problem type and adjusting its tolerance.
144: Arguments:
145: ctx - Pointer to the previously initialized DGP context
146: tol - Desired tolerance for the internal eigenvalue solver
147: */
148: PetscErrorCode DGPSetUpIterationWorkspace(DGPContext *ctx,PetscReal tol)
149: {
150: MPI_Comm comm;
152: PetscFunctionBeginUser;
153: comm = PetscObjectComm((PetscObject)ctx->A0);
154: /* Prepare the H matrix */
155: PetscCall(MatDuplicate(ctx->A0,MAT_DO_NOT_COPY_VALUES,&ctx->H));
157: /* Configure the eigenvector x */
158: PetscCall(MatCreateVecs(ctx->A0,&ctx->x,NULL));
160: /* Configure SLEPc EPS */
161: PetscCall(EPSCreate(comm,&ctx->eps));
162: PetscCall(EPSSetProblemType(ctx->eps,EPS_HEP));
163: PetscCall(EPSSetWhichEigenpairs(ctx->eps,EPS_SMALLEST_REAL));
164: PetscCall(EPSSetDimensions(ctx->eps,1,PETSC_DECIDE,PETSC_DECIDE));
165: PetscCall(EPSSetTolerances(ctx->eps,tol,PETSC_DECIDE));
166: PetscCall(EPSSetFromOptions(ctx->eps));
167: PetscFunctionReturn(PETSC_SUCCESS);
168: }
170: /*
171: Generate the initial guess X0 by solving the base linear Hamiltonian.
173: This function solves the eigenvalue problem (A0 * x = E * x) and fills
174: the vector x with the lowest energy eigenvector obtained. This
175: initial guess improves the convergence of the nonlinear iteration loop.
177: Arguments:
178: ctx - Pointer to the previously initialized DGP context
179: */
180: PetscErrorCode DGPGenerateInitialGuess(DGPContext *ctx)
181: {
182: PetscInt nconv;
183: MPI_Comm comm;
185: PetscFunctionBeginUser;
186: comm = PetscObjectComm((PetscObject)ctx->A0);
187: /* 1. Configure and solve the base linear problem A0 */
188: PetscCall(EPSSetOperators(ctx->eps,ctx->A0,NULL));
189: PetscCall(EPSSolve(ctx->eps));
191: /* Safety check: guarantee that SLEPc found the solution */
192: PetscCall(EPSGetConverged(ctx->eps,&nconv));
193: PetscCheck(nconv>=1,comm,PETSC_ERR_NOT_CONVERGED,"SLEPc could not find the initial linear ground state.");
195: /* 2. Extract the eigenvector and store it in x */
196: PetscCall(EPSGetEigenvector(ctx->eps,0,ctx->x,NULL));
197: PetscFunctionReturn(PETSC_SUCCESS);
198: }
200: /*
201: Compute the eigenvector-dependent term from the ground-state eigenvector.
203: In the context of the Gross-Pitaevskii model, this function computes the condensate
204: density rho(x) as the squared magnitude of the components of the eigenvector X0.
205: The formula used is: rho(x) = |x|.^2.
207: Arguments:
208: ctx - Pointer to the initialized DGP context
209: rho_out - Pre-created vector where the computed term will be stored
210: */
211: PetscErrorCode DGPComputeEigvecDependentTerm(DGPContext *ctx,Vec rho_out)
212: {
213: PetscInt i,n_loc;
214: const PetscScalar *x_local;
215: PetscScalar *rho_local;
217: PetscFunctionBeginUser;
218: PetscCall(VecSet(rho_out,0.0));
219: PetscCall(VecGetLocalSize(rho_out,&n_loc));
220: PetscCall(VecGetArray(rho_out,&rho_local));
222: PetscCall(VecGetArrayRead(ctx->x,&x_local));
224: for (i=0;i<n_loc;i++) {
225: /* Compute the squared magnitude: rho = Real^2 + Imag^2 */
226: rho_local[i]=PetscRealPart(x_local[i])*PetscRealPart(x_local[i])+PetscImaginaryPart(x_local[i])*PetscImaginaryPart(x_local[i]);
227: }
229: PetscCall(VecRestoreArrayRead(ctx->x,&x_local));
230: PetscCall(VecRestoreArray(rho_out,&rho_local));
231: PetscFunctionReturn(PETSC_SUCCESS);
232: }
234: /*
235: Build the Hamiltonian from an input eigenvector-dependent term.
237: This function assembles the updated Hamiltonian matrix stored in the context
238: using the formula: H = A0 + beta * Diag(rho_in).
240: Arguments:
241: ctx - Pointer to the initialized DGP context (ctx->H is updated)
242: rho_in - Vector with the proposed input eigenvector-dependent term
243: */
244: PetscErrorCode DGPBuildHamiltonian(DGPContext *ctx,Vec rho_in)
245: {
246: PetscFunctionBeginUser;
247: // 1. Build H = A0 + beta * Diag(rho_in)
248: PetscCall(MatCopy(ctx->A0,ctx->H,SAME_NONZERO_PATTERN));
249: PetscCall(VecCopy(rho_in,ctx->z));
251: PetscCall(VecScale(ctx->z,ctx->beta));
252: PetscCall(MatDiagonalSet(ctx->H,ctx->z,ADD_VALUES));
253: PetscFunctionReturn(PETSC_SUCCESS);
254: }
256: /*
257: Evaluation function (Callback) for the SNES nonlinear solver in DGP.
259: In each SNES iteration, this function receives a proposed eigenvector-dependent
260: term (rho_in), builds the Hamiltonian, solves the eigenvalue equation, computes the
261: resulting term, and returns the residual F = rho_out - rho_in.
263: Arguments:
264: snes - The nonlinear solver context
265: rho_in - Input state (eigenvector-dependent term) proposed by SNES
266: F - Vector where the computed residual will be stored
267: ctx_void - Pointer to the user's DGP context (DGPContext)
268: */
269: PetscErrorCode DGPNonlinearIteration(SNES snes,Vec rho_in,Vec F,void *ctx_void)
270: {
271: DGPContext *ctx=(DGPContext*)ctx_void;
272: PetscInt nconv;
273: MPI_Comm comm;
275: PetscFunctionBeginUser;
276: comm = PetscObjectComm((PetscObject)ctx->A0);
277: /* 1. Build the physics (Hamiltonian) with the guess proposed by SNES */
278: PetscCall(DGPBuildHamiltonian(ctx,rho_in));
280: /* 2. Solve the eigenvalue problem with the current H, using the previous state as initial guess */
281: PetscCall(EPSSetOperators(ctx->eps,ctx->H,NULL));
282: PetscCall(EPSSetInitialSpace(ctx->eps,1,&ctx->x));
283: PetscCall(EPSSolve(ctx->eps));
285: // Safety check: verify that SLEPc found the ground state
286: PetscCall(EPSGetConverged(ctx->eps,&nconv));
287: PetscCheck(nconv>=1,comm,PETSC_ERR_NOT_CONVERGED,"SLEPc did not find the ground state in this SNES iteration");
289: /* 3. Extract the new eigenvector and store it directly in ctx->x */
290: PetscCall(EPSGetEigenvector(ctx->eps,0,ctx->x,NULL));
292: /* 4. Compute the new eigenvector-dependent term (rho_out) generated by this eigenvector */
293: // Store the result in ctx->rho
294: PetscCall(DGPComputeEigvecDependentTerm(ctx,ctx->rho));
296: /* 5. Calculate the nonlinear residual: F = rho_in - rho_out */
297: PetscCall(VecWAXPY(F,-1.0,ctx->rho,rho_in));
298: PetscFunctionReturn(PETSC_SUCCESS);
299: }
301: /*
302: Free the memory associated with the DGP context.
304: This function destroys all internal PETSc objects created during
305: initialization and frees the memory of the main structure.
307: Arguments:
308: ctx - Pointer to the DGP context (set to NULL upon completion)
309: */
310: PetscErrorCode DGPDestroyContext(DGPContext **ctx)
311: {
312: PetscFunctionBeginUser;
313: if (!*ctx) PetscFunctionReturn(PETSC_SUCCESS);
315: PetscCall(MatDestroy(&(*ctx)->A0));
316: PetscCall(VecDestroy(&(*ctx)->rho));
317: PetscCall(VecDestroy(&(*ctx)->z));
319: PetscCall(MatDestroy(&(*ctx)->H));
320: PetscCall(VecDestroy(&(*ctx)->x));
321: PetscCall(EPSDestroy(&(*ctx)->eps));
323: PetscCall(PetscFree(*ctx));
324: PetscFunctionReturn(PETSC_SUCCESS);
325: }
327: int main(int argc,char **argv)
328: {
329: DGPContext *ctx;
330: Vec H_diag,rho_guess;
331: PetscScalar diagonal_sum,kr;
332: PetscInt n=3,maxit=100,its;
333: PetscInt nconv;
334: PetscReal beta=2.0,rtol=SLEPC_DEFAULT_TOL,stol=1e-12;
335: SNES snes;
336: Vec F; // Vector to store the residual (F = rho_out - rho_in)
337: PetscBool symm=PETSC_TRUE,verbose=PETSC_FALSE;
338: SNESConvergedReason reason;
339: SNESLineSearch linesearch;
341: PetscFunctionBeginUser;
342: PetscCall(SlepcInitialize(&argc,&argv,NULL,help));
344: PetscCheck(PetscDefined(USE_COMPLEX),PETSC_COMM_WORLD,PETSC_ERR_SUP,"This example requires complex scalars");
346: PetscCall(PetscPrintf(PETSC_COMM_WORLD,"--- DGP SNES NEPv ---\n"));
348: /* 1. Read DGP-specific command line options */
349: PetscCall(PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL));
350: PetscCall(PetscOptionsGetReal(NULL,NULL,"-beta",&beta,NULL));
351: PetscCall(PetscOptionsGetBool(NULL,NULL,"-symm",&symm,NULL));
353: /* 2. Initialization of the DGP context and base linear Hamiltonian A0 */
354: PetscCall(DGPCreateContext(PETSC_COMM_WORLD,n,beta,symm,&ctx));
356: /* 3. Configure objects for the nonlinear iteration (H matrix, EPS, Vec x) */
357: PetscCall(DGPSetUpIterationWorkspace(ctx,rtol*0.1));
359: /* 4. Generate the initial guess by solving the linear problem A0 */
360: PetscCall(DGPGenerateInitialGuess(ctx));
362: /* 5. Prepare the initial density (rho_guess) from the linear ground state */
363: PetscCall(VecDuplicate(ctx->rho,&rho_guess));
364: PetscCall(DGPComputeEigvecDependentTerm(ctx,rho_guess));
366: /* 6. Prepare the residual vector F by cloning rho */
367: PetscCall(VecDuplicate(ctx->rho,&F));
369: /* 7. Configure the nonlinear solver engine (SNES) with the DGP callback */
370: PetscCall(SNESCreate(PETSC_COMM_WORLD,&snes));
371: PetscCall(SNESSetFunction(snes,F,DGPNonlinearIteration,ctx));
372: PetscCall(SNESSetTolerances(snes,PETSC_DETERMINE,rtol,stol,maxit,PETSC_DETERMINE));
374: // Default -> NRICHARDSON with step lambda = 1.0 to simulate basic SCF
375: PetscCall(SNESSetType(snes,SNESNRICHARDSON));
376: PetscCall(SNESGetLineSearch(snes,&linesearch));
377: PetscCall(SNESLineSearchSetType(linesearch,SNESLINESEARCHNONE));
378: PetscCall(SNESLineSearchSetDamping(linesearch,1.0));
380: PetscCall(SNESSetFromOptions(snes));
382: PetscCall(SNESGetTolerances(snes,NULL,&rtol,NULL,&maxit,NULL));
383: PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Current parameters: n=%" PetscInt_FMT ", beta=%g, symm=%s, rtol=%g, maxit=%" PetscInt_FMT "\n",n,(double)beta,symm ? "True" : "False",(double)rtol,maxit));
385: /* 8. Nonlinear iteration loop */
386: PetscCall(PetscPrintf(PETSC_COMM_WORLD,"\nStarting nonlinear iterations...\n"));
387: PetscCall(SNESSolve(snes,NULL,rho_guess));
389: PetscCall(SNESGetConvergedReason(snes,&reason));
390: PetscCall(SNESGetIterationNumber(snes,&its));
392: /* 9. Analysis of results */
393: if (reason>0) {
394: PetscCall(PetscPrintf(PETSC_COMM_WORLD,"--> CONVERGENCE REACHED in %" PetscInt_FMT " iterations (Reason: %s).\n",its,SNESConvergedReasons[reason]));
396: /* Extract the eigenvalue of the converged ground state */
397: PetscCall(EPSGetConverged(ctx->eps,&nconv));
398: if (nconv>0) {
399: PetscCall(EPSGetEigenvalue(ctx->eps,0,&kr,NULL));
400: PetscCall(PetscPrintf(PETSC_COMM_WORLD,"\nEigenvalue found = %10.6f\n",(double)PetscRealPart(kr)));
401: }
403: PetscCall(PetscOptionsHasName(NULL,NULL,"-verbose",&verbose));
404: if (verbose) {
405: PetscCall(MatCreateVecs(ctx->H,NULL,&H_diag));
406: PetscCall(MatGetDiagonal(ctx->H,H_diag));
407: PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Diagonal of the final Hamiltonian:\n"));
408: PetscCall(VecView(H_diag,PETSC_VIEWER_STDOUT_WORLD));
410: PetscCall(VecSum(H_diag,&diagonal_sum));
411: PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Sum of the diagonal (Approx total energy): %g\n",(double)PetscRealPart(diagonal_sum)));
412: PetscCall(VecDestroy(&H_diag));
413: }
415: } else {
416: PetscCall(PetscPrintf(PETSC_COMM_WORLD,"--> ERROR: SNES did not converge (Reason: %s)\n",SNESConvergedReasons[reason]));
417: }
419: /* 10. Clean up memory */
420: PetscCall(VecDestroy(&F));
421: PetscCall(VecDestroy(&rho_guess));
422: PetscCall(SNESDestroy(&snes));
423: PetscCall(DGPDestroyContext(&ctx));
425: PetscCall(SlepcFinalize());
426: return 0;
427: }
429: /*TEST
431: build:
432: requires: complex
434: testset:
435: filter: sed -e "s/rtol=1e-05/rtol=1e-08/" -e "s/rtol=1e-16/rtol=1e-08/" -e "s/[0-9]\{1,\} iterations/23 iterations/" -e "s/CONVERGED_SNORM_RELATIVE/CONVERGED_FNORM_RELATIVE/"
436: output_file: output/ex60_1.out
437: test:
438: suffix: 1
439: test:
440: suffix: 2
441: args: -snes_type anderson -snes_anderson_m 7
442: test:
443: suffix: 3
444: args: -snes_type ngmres -npc_snes_type nrichardson -snes_npc_side right
445: test:
446: suffix: 4
447: args: -snes_type composite -snes_composite_type additiveoptimal -snes_composite_sneses anderson,nrichardson -sub_0_snes_anderson_m 7 -sub_0_snes_anderson_beta 0.4
449: TEST*/