GCC Code Coverage Report

Directory: ./
File: src/eps/tutorials/ex60.c
Date: 2026-07-29 03:58:07
Exec Total Coverage
Lines: 159 169 94.1%
Functions: 8 8 100.0%
Branches: 400 708 56.5%
Line Branch Exec Source
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 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";
24
25 #include <slepceps.h>
26
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 self-consistent field (SCF) iteration.
30 */
31 typedef struct {
32 // DGP context
33
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
38
39 // Constants
40 PetscReal L; // Length of the semi-domain [-L, L]
41 PetscReal h; // Step size
42 PetscReal omega; // Angular frequency
43
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
48
49 // SCF context
50 EPS eps; // Solver for each H
51 Mat H; // Hamiltonian
52 Vec x; // Eigenvector
53 } DGPContext;
54
55 /*
56 Initialize the context for the 2D discrete Gross-Pitaevskii problem.
57
58 This function allocates the necessary memory and builds the base linear
59 Hamiltonian A0 as
60
61 A0 = h^2 * ( -0.5*L + V - omega*1i*L_z )
62
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.
68
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 20 PetscErrorCode DGPCreate(MPI_Comm comm,PetscInt N,PetscReal beta,PetscBool symm,DGPContext **ctx_out)
76 {
77 20 DGPContext *ctx;
78
79 20 PetscInt Istart,Iend,i,II,j;
80 20 PetscReal h2;
81 20 PetscReal x,y,V;
82 20 PetscScalar alpha;
83
84
1/2
✓ Branch 0 taken 1 times.
✗ Branch 1 not taken.
 20 PetscFunctionBeginUser;
85
4/6
✓ Branch 0 taken 1 times.
✓ Branch 1 taken 4 times.
✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 1 times.
 20 PetscCall(PetscNew(&ctx));
86
87 // Inputs
88 20 ctx->N=N;
89 20 ctx->beta=beta;
90 20 ctx->symm=symm;
91
92 // Constants
93 20 ctx->L=1.0;
94 20 ctx->h=2.0*ctx->L/(ctx->N+1.0);
95 20 ctx->omega=0.85;
96 20 h2=ctx->h*ctx->h;
97
98 // 1. Create A0
99
4/6
✓ Branch 0 taken 1 times.
✓ Branch 1 taken 4 times.
✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 1 times.
 20 PetscCall(MatCreate(comm,&ctx->A0));
100
4/6
✓ Branch 0 taken 1 times.
✓ Branch 1 taken 4 times.
✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 1 times.
 20 PetscCall(MatSetSizes(ctx->A0,PETSC_DECIDE,PETSC_DECIDE,N*N,N*N));
101
4/6
✓ Branch 0 taken 1 times.
✓ Branch 1 taken 4 times.
✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 1 times.
 20 PetscCall(MatSetFromOptions(ctx->A0));
102
103 // 2. Fill A0
104
4/6
✓ Branch 0 taken 1 times.
✓ Branch 1 taken 4 times.
✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 1 times.
 20 PetscCall(MatGetOwnershipRange(ctx->A0,&Istart,&Iend));
105
106
2/2
✓ Branch 0 taken 5 times.
✓ Branch 1 taken 5 times.
 200 for (II=Istart;II<Iend;II++) {
107 180 i=II/ctx->N; j=II-i*ctx->N;
108
109 // Potential for mesh point with coordinates (x,y)
110 180 x=-ctx->L+(j+1)*ctx->h;
111 180 y=-ctx->L+(i+1)*ctx->h;
112
1/2
✓ Branch 0 taken 5 times.
✗ Branch 1 not taken.
 180 if (ctx->symm) V=0.5*(x*x+y*y);
113 ✗ else V=0.5*(x*x+100.0*y*y);
114
115 180 alpha = -h2*ctx->omega/(2.0*ctx->h)*PETSC_i;
116
4/6
✓ Branch 0 taken 1 times.
✓ Branch 1 taken 4 times.
✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 1 times.
 180 PetscCall(MatSetValue(ctx->A0,II,II,0.5*4.0+h2*V,INSERT_VALUES));
117
6/8
✓ Branch 0 taken 5 times.
✓ Branch 1 taken 5 times.
✓ Branch 2 taken 1 times.
✓ Branch 3 taken 4 times.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✓ Branch 7 taken 1 times.
 180 if (i>0) PetscCall(MatSetValue(ctx->A0,II,II-ctx->N,-0.5+alpha*x,INSERT_VALUES));
118
6/8
✓ Branch 0 taken 5 times.
✓ Branch 1 taken 5 times.
✓ Branch 2 taken 1 times.
✓ Branch 3 taken 4 times.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✓ Branch 7 taken 1 times.
 180 if (i<ctx->N-1) PetscCall(MatSetValue(ctx->A0,II,II+ctx->N,-0.5-alpha*x,INSERT_VALUES));
119
6/8
✓ Branch 0 taken 5 times.
✓ Branch 1 taken 5 times.
✓ Branch 2 taken 1 times.
✓ Branch 3 taken 4 times.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✓ Branch 7 taken 1 times.
 180 if (j>0) PetscCall(MatSetValue(ctx->A0,II,II-1,-0.5-alpha*y,INSERT_VALUES));
120
6/8
✓ Branch 0 taken 5 times.
✓ Branch 1 taken 5 times.
✓ Branch 2 taken 1 times.
✓ Branch 3 taken 4 times.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✓ Branch 7 taken 1 times.
 180 if (j<ctx->N-1) PetscCall(MatSetValue(ctx->A0,II,II+1,-0.5+alpha*y,INSERT_VALUES));
121 }
122
123
4/6
✓ Branch 0 taken 1 times.
✓ Branch 1 taken 4 times.
✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 1 times.
 20 PetscCall(MatAssemblyBegin(ctx->A0,MAT_FINAL_ASSEMBLY));
124
4/6
✓ Branch 0 taken 1 times.
✓ Branch 1 taken 4 times.
✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 1 times.
 20 PetscCall(MatAssemblyEnd(ctx->A0,MAT_FINAL_ASSEMBLY));
125
126 // 3. Create vectors to store the condensate density
127
4/6
✓ Branch 0 taken 1 times.
✓ Branch 1 taken 4 times.
✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 1 times.
 20 PetscCall(MatCreateVecs(ctx->A0,&ctx->rho,&ctx->z));
128
129 20 *ctx_out=ctx;
130
5/12
✓ Branch 0 taken 1 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 1 times.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 1 times.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✓ Branch 9 taken 1 times.
✗ Branch 10 not taken.
✗ Branch 11 not taken.
 20 PetscFunctionReturn(PETSC_SUCCESS);
131 }
132
133 /*
134 Configure the SLEPc objects necessary for the SCF iteration.
135
136 This function prepares the Hamiltonian matrix by cloning the structure of A0,
137 initializes the vector to store the condensate state, and configures the eigenvalue
138 solver (EPS) by defining the problem type and adjusting its tolerance.
139
140 Arguments:
141 ctx - Pointer to the previously initialized DGP context
142 tol - Desired tolerance for the internal eigenvalue solver
143 */
144 20 PetscErrorCode DGPSetupSCF(DGPContext *ctx,PetscReal tol)
145 {
146 20 MPI_Comm comm;
147
148
1/2
✓ Branch 0 taken 1 times.
✗ Branch 1 not taken.
 20 PetscFunctionBeginUser;
149 20 comm = PetscObjectComm((PetscObject)ctx->A0);
150 /* Prepare the H matrix */
151
4/6
✓ Branch 0 taken 1 times.
✓ Branch 1 taken 4 times.
✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 1 times.
 20 PetscCall(MatDuplicate(ctx->A0,MAT_DO_NOT_COPY_VALUES,&ctx->H));
152
153 /* Configure the eigenvector x */
154
4/6
✓ Branch 0 taken 1 times.
✓ Branch 1 taken 4 times.
✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 1 times.
 20 PetscCall(MatCreateVecs(ctx->A0,&ctx->x,NULL));
155
156 /* Configure SLEPc EPS */
157
4/6
✓ Branch 0 taken 1 times.
✓ Branch 1 taken 4 times.
✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 1 times.
 20 PetscCall(EPSCreate(comm,&ctx->eps));
158
4/6
✓ Branch 0 taken 1 times.
✓ Branch 1 taken 4 times.
✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 1 times.
 20 PetscCall(EPSSetProblemType(ctx->eps,EPS_HEP));
159
4/6
✓ Branch 0 taken 1 times.
✓ Branch 1 taken 4 times.
✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 1 times.
 20 PetscCall(EPSSetWhichEigenpairs(ctx->eps,EPS_SMALLEST_REAL));
160
4/6
✓ Branch 0 taken 1 times.
✓ Branch 1 taken 4 times.
✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 1 times.
 20 PetscCall(EPSSetDimensions(ctx->eps,1,PETSC_DECIDE,PETSC_DECIDE));
161
4/6
✓ Branch 0 taken 1 times.
✓ Branch 1 taken 4 times.
✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 1 times.
 20 PetscCall(EPSSetTolerances(ctx->eps,tol,PETSC_DECIDE));
162
4/6
✓ Branch 0 taken 1 times.
✓ Branch 1 taken 4 times.
✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 1 times.
 20 PetscCall(EPSSetFromOptions(ctx->eps));
163
5/12
✓ Branch 0 taken 1 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 1 times.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 1 times.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✓ Branch 9 taken 1 times.
✗ Branch 10 not taken.
✗ Branch 11 not taken.
 4 PetscFunctionReturn(PETSC_SUCCESS);
164 }
165
166 /*
167 Generate the initial guess X0 by solving the base linear Hamiltonian.
168
169 This function solves the eigenvalue problem (A0 * x = E * x) and fills
170 the vector x with the lowest energy eigenvector obtained. This
171 initial guess improves the convergence of the SCF loop.
172
173 Arguments:
174 ctx - Pointer to the previously initialized DGP context
175 */
176 20 PetscErrorCode DGPGenerateInitialGuess(DGPContext *ctx)
177 {
178 20 PetscInt nconv;
179 20 MPI_Comm comm;
180
181
1/2
✓ Branch 0 taken 1 times.
✗ Branch 1 not taken.
 20 PetscFunctionBeginUser;
182 20 comm = PetscObjectComm((PetscObject)ctx->A0);
183 /* 1. Configure and solve the base linear problem A0 */
184
4/6
✓ Branch 0 taken 1 times.
✓ Branch 1 taken 4 times.
✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 1 times.
 20 PetscCall(EPSSetOperators(ctx->eps,ctx->A0,NULL));
185
4/6
✓ Branch 0 taken 1 times.
✓ Branch 1 taken 4 times.
✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 1 times.
 20 PetscCall(EPSSolve(ctx->eps));
186
187 /* Safety check: guarantee that SLEPc found the solution */
188
4/6
✓ Branch 0 taken 1 times.
✓ Branch 1 taken 4 times.
✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 1 times.
 20 PetscCall(EPSGetConverged(ctx->eps,&nconv));
189
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 5 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
 20 PetscCheck(nconv>=1,comm,PETSC_ERR_NOT_CONVERGED,"SLEPc could not find the initial linear ground state.");
190
191 /* 2. Extract the eigenvector and store it in x */
192
4/6
✓ Branch 0 taken 1 times.
✓ Branch 1 taken 4 times.
✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 1 times.
 20 PetscCall(EPSGetEigenvector(ctx->eps,0,ctx->x,NULL));
193
5/12
✓ Branch 0 taken 1 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 1 times.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 1 times.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✓ Branch 9 taken 1 times.
✗ Branch 10 not taken.
✗ Branch 11 not taken.
 4 PetscFunctionReturn(PETSC_SUCCESS);
194 }
195
196 /*
197 Compute the condensate density from the ground-state eigenvector.
198
199 This function computes the density rho(x) as the squared magnitude
200 of the components of the eigenvector X0. The formula used is: rho(x) = |x|.^2.
201
202 Arguments:
203 ctx - Pointer to the initialized DGP context
204 rho_out - Pre-created vector where the computed density will be stored
205 */
206 410 PetscErrorCode DGPCalculateDensity(DGPContext *ctx,Vec rho_out)
207 {
208 410 PetscInt i,n_loc;
209 410 const PetscScalar *x_local;
210 410 PetscScalar *rho_local;
211
212
1/2
✓ Branch 0 taken 1 times.
✗ Branch 1 not taken.
 410 PetscFunctionBeginUser;
213
4/6
✓ Branch 0 taken 1 times.
✓ Branch 1 taken 4 times.
✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 1 times.
 410 PetscCall(VecSet(rho_out,0.0));
214
4/6
✓ Branch 0 taken 1 times.
✓ Branch 1 taken 4 times.
✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 1 times.
 410 PetscCall(VecGetLocalSize(rho_out,&n_loc));
215
4/6
✓ Branch 0 taken 1 times.
✓ Branch 1 taken 4 times.
✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 1 times.
 410 PetscCall(VecGetArray(rho_out,&rho_local));
216
217
4/6
✓ Branch 0 taken 1 times.
✓ Branch 1 taken 4 times.
✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 1 times.
 410 PetscCall(VecGetArrayRead(ctx->x,&x_local));
218
219
2/2
✓ Branch 0 taken 5 times.
✓ Branch 1 taken 5 times.
 4100 for (i=0;i<n_loc;i++) {
220 /* Compute the squared magnitude: rho = Real^2 + Imag^2 */
221 3690 PetscReal mod=PetscAbsScalar(x_local[i]);
222 3690 rho_local[i]=mod*mod;
223 }
224
225
4/6
✓ Branch 0 taken 1 times.
✓ Branch 1 taken 4 times.
✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 1 times.
 410 PetscCall(VecRestoreArrayRead(ctx->x,&x_local));
226
4/6
✓ Branch 0 taken 1 times.
✓ Branch 1 taken 4 times.
✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 1 times.
 410 PetscCall(VecRestoreArray(rho_out,&rho_local));
227
5/12
✓ Branch 0 taken 1 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 1 times.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 1 times.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✓ Branch 9 taken 1 times.
✗ Branch 10 not taken.
✗ Branch 11 not taken.
 82 PetscFunctionReturn(PETSC_SUCCESS);
228 }
229
230 /*
231 Build the Hamiltonian from an input density.
232
233 This function assembles the updated Hamiltonian matrix stored in the context
234 using the formula: H = A0 + beta * Diag(rho_in).
235
236 Arguments:
237 ctx - Pointer to the initialized DGP context (ctx->H is updated)
238 rho_in - Vector with the input density
239 */
240 390 PetscErrorCode DGPBuildHamiltonian(DGPContext *ctx,Vec rho_in)
241 {
242
1/2
✓ Branch 0 taken 1 times.
✗ Branch 1 not taken.
 390 PetscFunctionBeginUser;
243 // 1. Build H = A0 + beta * Diag(rho_in)
244
4/6
✓ Branch 0 taken 1 times.
✓ Branch 1 taken 4 times.
✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 1 times.
 390 PetscCall(MatCopy(ctx->A0,ctx->H,SAME_NONZERO_PATTERN));
245
4/6
✓ Branch 0 taken 1 times.
✓ Branch 1 taken 4 times.
✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 1 times.
 390 PetscCall(VecCopy(rho_in,ctx->z));
246
247
4/6
✓ Branch 0 taken 1 times.
✓ Branch 1 taken 4 times.
✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 1 times.
 390 PetscCall(VecScale(ctx->z,ctx->beta));
248
4/6
✓ Branch 0 taken 1 times.
✓ Branch 1 taken 4 times.
✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 1 times.
 390 PetscCall(MatDiagonalSet(ctx->H,ctx->z,ADD_VALUES));
249
5/12
✓ Branch 0 taken 1 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 1 times.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 1 times.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✓ Branch 9 taken 1 times.
✗ Branch 10 not taken.
✗ Branch 11 not taken.
 78 PetscFunctionReturn(PETSC_SUCCESS);
250 }
251
252 /*
253 Evaluation function (Callback) for the SNES nonlinear solver in DGP.
254
255 In each SNES iteration, this function receives a proposed density (rho_in),
256 builds the Hamiltonian, solves the eigenvalue equation, computes the
257 resulting density, and returns the residual F = rho_out - rho_in.
258
259 Arguments:
260 snes - The nonlinear solver context
261 rho_in - Condensate density proposed by SNES
262 F - Vector where the computed residual will be stored
263 ctx_void - Pointer to the user's DGP context (DGPContext)
264 */
265 390 PetscErrorCode DGPIterationSCF(SNES snes,Vec rho_in,Vec F,void *ctx_void)
266 {
267 390 DGPContext *ctx=(DGPContext*)ctx_void;
268 390 PetscInt nconv;
269 390 MPI_Comm comm;
270
271
1/2
✓ Branch 0 taken 1 times.
✗ Branch 1 not taken.
 390 PetscFunctionBeginUser;
272 390 comm = PetscObjectComm((PetscObject)ctx->A0);
273 /* 1. Build the physics (Hamiltonian) with the guess proposed by SNES */
274
4/6
✓ Branch 0 taken 1 times.
✓ Branch 1 taken 4 times.
✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 1 times.
 390 PetscCall(DGPBuildHamiltonian(ctx,rho_in));
275
276 /* 2. Solve the eigenvalue problem with the current H */
277
4/6
✓ Branch 0 taken 1 times.
✓ Branch 1 taken 4 times.
✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 1 times.
 390 PetscCall(EPSSetOperators(ctx->eps,ctx->H,NULL));
278
4/6
✓ Branch 0 taken 1 times.
✓ Branch 1 taken 4 times.
✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 1 times.
 390 PetscCall(EPSSolve(ctx->eps));
279
280 // Safety check: verify that SLEPc found the ground state
281
4/6
✓ Branch 0 taken 1 times.
✓ Branch 1 taken 4 times.
✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 1 times.
 390 PetscCall(EPSGetConverged(ctx->eps,&nconv));
282
1/8
✗ Branch 0 not taken.
✓ Branch 1 taken 5 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
 390 if (nconv<1) PetscCall(PetscPrintf(comm,"Warning: SLEPc did not find the ground state in this iteration.\n"));
283
284 /* 3. Extract the new eigenvector and store it directly in ctx->x */
285
4/6
✓ Branch 0 taken 1 times.
✓ Branch 1 taken 4 times.
✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 1 times.
 390 PetscCall(EPSGetEigenvector(ctx->eps,0,ctx->x,NULL));
286
287 /* 4. Compute the new density (rho_out) generated by this eigenvector */
288 // Store the result in ctx->rho
289
4/6
✓ Branch 0 taken 1 times.
✓ Branch 1 taken 4 times.
✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 1 times.
 390 PetscCall(DGPCalculateDensity(ctx,ctx->rho));
290
291 /* 5. Calculate the nonlinear residual: F = rho_out - rho_in */
292
4/6
✓ Branch 0 taken 1 times.
✓ Branch 1 taken 4 times.
✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 1 times.
 390 PetscCall(VecWAXPY(F,-1.0,ctx->rho,rho_in));
293
5/12
✓ Branch 0 taken 1 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 1 times.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 1 times.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✓ Branch 9 taken 1 times.
✗ Branch 10 not taken.
✗ Branch 11 not taken.
 78 PetscFunctionReturn(PETSC_SUCCESS);
294 }
295
296 /*
297 Free the memory associated with the DGP context.
298
299 This function destroys all internal PETSc objects created during
300 initialization and frees the memory of the main structure.
301
302 Arguments:
303 ctx - Pointer to the DGP context (set to NULL upon completion)
304 */
305 20 PetscErrorCode DGPDestroy(DGPContext **ctx)
306 {
307
1/2
✓ Branch 0 taken 1 times.
✗ Branch 1 not taken.
 20 PetscFunctionBeginUser;
308
2/14
✓ Branch 0 taken 4 times.
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
✗ Branch 13 not taken.
 20 if (!*ctx) PetscFunctionReturn(PETSC_SUCCESS);
309
310
4/6
✓ Branch 0 taken 1 times.
✓ Branch 1 taken 4 times.
✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 1 times.
 20 PetscCall(MatDestroy(&(*ctx)->A0));
311
4/6
✓ Branch 0 taken 1 times.
✓ Branch 1 taken 4 times.
✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 1 times.
 20 PetscCall(VecDestroy(&(*ctx)->rho));
312
4/6
✓ Branch 0 taken 1 times.
✓ Branch 1 taken 4 times.
✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 1 times.
 20 PetscCall(VecDestroy(&(*ctx)->z));
313
314
4/6
✓ Branch 0 taken 1 times.
✓ Branch 1 taken 4 times.
✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 1 times.
 20 PetscCall(MatDestroy(&(*ctx)->H));
315
4/6
✓ Branch 0 taken 1 times.
✓ Branch 1 taken 4 times.
✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 1 times.
 20 PetscCall(VecDestroy(&(*ctx)->x));
316
4/6
✓ Branch 0 taken 1 times.
✓ Branch 1 taken 4 times.
✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 1 times.
 20 PetscCall(EPSDestroy(&(*ctx)->eps));
317
318
5/8
✓ Branch 0 taken 5 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 1 times.
✓ Branch 3 taken 4 times.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✓ Branch 7 taken 1 times.
 20 PetscCall(PetscFree(*ctx));
319
5/12
✓ Branch 0 taken 1 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 1 times.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 1 times.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✓ Branch 9 taken 1 times.
✗ Branch 10 not taken.
✗ Branch 11 not taken.
 4 PetscFunctionReturn(PETSC_SUCCESS);
320 }
321
322 20 int main(int argc,char **argv)
323 {
324 20 DGPContext *ctx;
325 20 Vec H_diag,rho_guess;
326 20 PetscScalar diagonal_sum,kr;
327 20 PetscInt n=3,maxit=100,its;
328 20 PetscInt nconv;
329 20 PetscReal beta=2.0,rtol=SLEPC_DEFAULT_TOL,stol=1e-12;
330 20 SNES snes;
331 20 Vec F; // Vector to store the residual (F = rho_out - rho_in)
332 20 PetscBool symm=PETSC_TRUE,verbose=PETSC_FALSE;
333 20 SNESConvergedReason reason;
334 20 SNESLineSearch linesearch;
335
336
1/2
✓ Branch 0 taken 1 times.
✗ Branch 1 not taken.
 20 PetscFunctionBeginUser;
337
4/6
✓ Branch 0 taken 1 times.
✓ Branch 1 taken 4 times.
✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 1 times.
 20 PetscCall(SlepcInitialize(&argc,&argv,NULL,help));
338
339 20 PetscCheck(PetscDefined(USE_COMPLEX),PETSC_COMM_WORLD,PETSC_ERR_SUP,"This example requires complex scalars");
340
341
4/6
✓ Branch 0 taken 1 times.
✓ Branch 1 taken 4 times.
✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 1 times.
 20 PetscCall(PetscPrintf(PETSC_COMM_WORLD,"--- DGP SNES SCF ---\n"));
342
343 /* 1. Read DGP-specific command line options */
344
4/6
✓ Branch 0 taken 1 times.
✓ Branch 1 taken 4 times.
✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 1 times.
 20 PetscCall(PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL));
345
4/6
✓ Branch 0 taken 1 times.
✓ Branch 1 taken 4 times.
✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 1 times.
 20 PetscCall(PetscOptionsGetReal(NULL,NULL,"-beta",&beta,NULL));
346
4/6
✓ Branch 0 taken 1 times.
✓ Branch 1 taken 4 times.
✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 1 times.
 20 PetscCall(PetscOptionsGetBool(NULL,NULL,"-symm",&symm,NULL));
347
348 /* 2. Initialization of the DGP context and base linear Hamiltonian A0 */
349
4/6
✓ Branch 0 taken 1 times.
✓ Branch 1 taken 4 times.
✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 1 times.
 20 PetscCall(DGPCreate(PETSC_COMM_WORLD,n,beta,symm,&ctx));
350
351 /* 3. Configure objects for the SCF (H matrix, EPS, Vec x) */
352
4/6
✓ Branch 0 taken 1 times.
✓ Branch 1 taken 4 times.
✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 1 times.
 20 PetscCall(DGPSetupSCF(ctx,rtol*0.1));
353
354 /* 4. Generate the initial guess by solving the linear problem A0 */
355
4/6
✓ Branch 0 taken 1 times.
✓ Branch 1 taken 4 times.
✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 1 times.
 20 PetscCall(DGPGenerateInitialGuess(ctx));
356
357 /* 5. Prepare the initial density (rho_guess) from the linear ground state */
358
4/6
✓ Branch 0 taken 1 times.
✓ Branch 1 taken 4 times.
✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 1 times.
 20 PetscCall(VecDuplicate(ctx->rho,&rho_guess));
359
4/6
✓ Branch 0 taken 1 times.
✓ Branch 1 taken 4 times.
✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 1 times.
 20 PetscCall(DGPCalculateDensity(ctx,rho_guess));
360
361 /* 6. Prepare the residual vector F by cloning rho */
362
4/6
✓ Branch 0 taken 1 times.
✓ Branch 1 taken 4 times.
✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 1 times.
 20 PetscCall(VecDuplicate(ctx->rho,&F));
363
364 /* 7. Configure the nonlinear solver engine (SNES) with the DGP callback */
365
4/6
✓ Branch 0 taken 1 times.
✓ Branch 1 taken 4 times.
✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 1 times.
 20 PetscCall(SNESCreate(PETSC_COMM_WORLD,&snes));
366
4/6
✓ Branch 0 taken 1 times.
✓ Branch 1 taken 4 times.
✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 1 times.
 20 PetscCall(SNESSetFunction(snes,F,DGPIterationSCF,ctx));
367
4/6
✓ Branch 0 taken 1 times.
✓ Branch 1 taken 4 times.
✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 1 times.
 20 PetscCall(SNESSetTolerances(snes,PETSC_DETERMINE,rtol,stol,maxit,PETSC_DETERMINE));
368
369 // Default -> NRICHARDSON with step lambda = 1.0 to simulate basic SCF
370
4/6
✓ Branch 0 taken 1 times.
✓ Branch 1 taken 4 times.
✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 1 times.
 20 PetscCall(SNESSetType(snes,SNESNRICHARDSON));
371
4/6
✓ Branch 0 taken 1 times.
✓ Branch 1 taken 4 times.
✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 1 times.
 20 PetscCall(SNESGetLineSearch(snes,&linesearch));
372
4/6
✓ Branch 0 taken 1 times.
✓ Branch 1 taken 4 times.
✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 1 times.
 20 PetscCall(SNESLineSearchSetType(linesearch,SNESLINESEARCHNONE));
373
4/6
✓ Branch 0 taken 1 times.
✓ Branch 1 taken 4 times.
✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 1 times.
 20 PetscCall(SNESLineSearchSetDamping(linesearch,1.0));
374
375
4/6
✓ Branch 0 taken 1 times.
✓ Branch 1 taken 4 times.
✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 1 times.
 20 PetscCall(SNESSetFromOptions(snes));
376
377
4/6
✓ Branch 0 taken 1 times.
✓ Branch 1 taken 4 times.
✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 1 times.
 20 PetscCall(SNESGetTolerances(snes,NULL,&rtol,NULL,&maxit,NULL));
378
6/8
✓ Branch 0 taken 1 times.
✓ Branch 1 taken 4 times.
✓ Branch 2 taken 1 times.
✓ Branch 3 taken 4 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 1 times.
✗ Branch 6 not taken.
✓ Branch 7 taken 1 times.
 20 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));
379
380 /* 8. SCF loop */
381
4/6
✓ Branch 0 taken 1 times.
✓ Branch 1 taken 4 times.
✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 1 times.
 20 PetscCall(PetscPrintf(PETSC_COMM_WORLD,"\nStarting SCF iterations...\n"));
382
4/6
✓ Branch 0 taken 1 times.
✓ Branch 1 taken 4 times.
✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 1 times.
 20 PetscCall(SNESSolve(snes,NULL,rho_guess));
383
384
4/6
✓ Branch 0 taken 1 times.
✓ Branch 1 taken 4 times.
✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 1 times.
 20 PetscCall(SNESGetConvergedReason(snes,&reason));
385
4/6
✓ Branch 0 taken 1 times.
✓ Branch 1 taken 4 times.
✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 1 times.
 20 PetscCall(SNESGetIterationNumber(snes,&its));
386
387 /* 9. Analysis of results */
388
1/2
✓ Branch 0 taken 5 times.
✗ Branch 1 not taken.
 20 if (reason>0) {
389
4/6
✓ Branch 0 taken 1 times.
✓ Branch 1 taken 4 times.
✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 1 times.
 20 PetscCall(PetscPrintf(PETSC_COMM_WORLD,"--> CONVERGENCE REACHED in %" PetscInt_FMT " iterations (Reason: %s).\n",its,SNESConvergedReasons[reason]));
390
391 /* Extract the eigenvalue of the converged ground state */
392
4/6
✓ Branch 0 taken 1 times.
✓ Branch 1 taken 4 times.
✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 1 times.
 20 PetscCall(EPSGetConverged(ctx->eps,&nconv));
393
1/2
✓ Branch 0 taken 5 times.
✗ Branch 1 not taken.
 20 if (nconv>0) {
394
4/6
✓ Branch 0 taken 1 times.
✓ Branch 1 taken 4 times.
✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 1 times.
 20 PetscCall(EPSGetEigenvalue(ctx->eps,0,&kr,NULL));
395
4/6
✓ Branch 0 taken 1 times.
✓ Branch 1 taken 4 times.
✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 1 times.
 20 PetscCall(PetscPrintf(PETSC_COMM_WORLD,"\nEigenvalue found = %10.6f\n",(double)PetscRealPart(kr)));
396 }
397
398
4/6
✓ Branch 0 taken 1 times.
✓ Branch 1 taken 4 times.
✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 1 times.
 20 PetscCall(PetscOptionsHasName(NULL,NULL,"-verbose",&verbose));
399
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 5 times.
 20 if (verbose) {
400 ✗ PetscCall(MatCreateVecs(ctx->H,NULL,&H_diag));
401 ✗ PetscCall(MatGetDiagonal(ctx->H,H_diag));
402 ✗ PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Diagonal of the final Hamiltonian:\n"));
403 ✗ PetscCall(VecView(H_diag,PETSC_VIEWER_STDOUT_WORLD));
404
405 ✗ PetscCall(VecSum(H_diag,&diagonal_sum));
406 ✗ PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Sum of the diagonal (Approx total energy): %g\n",(double)PetscRealPart(diagonal_sum)));
407 ✗ PetscCall(VecDestroy(&H_diag));
408 }
409
410 } else {
411 ✗ PetscCall(SNESGetConvergedReason(snes,&reason));
412 ✗ PetscCall(PetscPrintf(PETSC_COMM_WORLD,"--> ERROR: The SCF did not converge (Reason: %s)\n",SNESConvergedReasons[reason]));
413 }
414
415 /* 10. Clean up memory */
416
4/6
✓ Branch 0 taken 1 times.
✓ Branch 1 taken 4 times.
✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 1 times.
 20 PetscCall(VecDestroy(&F));
417
4/6
✓ Branch 0 taken 1 times.
✓ Branch 1 taken 4 times.
✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 1 times.
 20 PetscCall(VecDestroy(&rho_guess));
418
4/6
✓ Branch 0 taken 1 times.
✓ Branch 1 taken 4 times.
✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 1 times.
 20 PetscCall(SNESDestroy(&snes));
419
4/6
✓ Branch 0 taken 1 times.
✓ Branch 1 taken 4 times.
✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 1 times.
 20 PetscCall(DGPDestroy(&ctx));
420
421
2/6
✓ Branch 0 taken 1 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
 20 PetscCall(SlepcFinalize());
422 return 0;
423 }
424
425 /*TEST
426
427 build:
428 requires: complex
429
430 testset:
431 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/"
432 output_file: output/ex60_1.out
433 test:
434 suffix: 1
435 test:
436 suffix: 2
437 args: -snes_type anderson -snes_anderson_m 7
438 test:
439 suffix: 3
440 args: -snes_type ngmres -npc_snes_type nrichardson -snes_npc_side right
441 test:
442 suffix: 4
443 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
444
445 TEST*/
446