examples/ShellMatrix_petsc_Lindbladian

examples/ShellMatrix_petsc_Lindbladian/src/main.cpp

  1#include <iostream>
  2#include <mpi.h>
  3#include <complex>
  4
  5/* slepceps includes eigen value routines */
  6#include <slepc.h> 
  7
  8/* Including Eigen */
  9#include <Eigen/Dense>
 10#include <Eigen/Eigenvalues>
 11
 12/* Time measurement */
 13#include <chrono>
 14
 15/* Maximal system size */
 16#define MaxSites 64
 17
 18/* Hilbert space dimension */
 19#define Q 2
 20
 21/* ScalarType */
 22#define ScalarType std::complex<double>
 23
 24/* Monitoring function */
 25PetscErrorCode monitor(EPS eps,PetscInt its,PetscInt nconv,PetscScalar *eigr,PetscScalar *eigi,PetscReal *errest,PetscInt nest,void *mctx){
 26
 27	PetscMPIInt myrank;
 28	PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &myrank));
 29
 30	if (myrank == 0){
 31		std::cout << "[myrank: " << myrank << "] - [monitor] - iteration = " << its << " - nconv = " << nconv << " - errors: ";
 32		for(PetscInt c_for = 0; c_for < nconv+1; c_for++){
 33			std::cout << errest[c_for] << " ";
 34		}
 35		std::cout << std::endl;
 36	}
 37	PetscErrorCode ierr = PETSC_SUCCESS;
 38	return ierr;
 39}
 40
 41#include "DanceQ.h"
 42using namespace danceq;
 43
 44int main(int argc, char *argv[]) {
 45
 46
 47
 48	PetscInitialize(&argc, &argv, NULL, NULL);
 49
 50	/* System size */
 51	PetscInt L = 6;
 52
 53	/* Anisotropy */
 54	PetscReal Jz = 1.;
 55
 56	/* Number of eigenvalues */
 57	PetscInt nev = 2;
 58
 59	/* Maximal dimension of subspace */
 60	PetscInt ncv = 100;
 61
 62	/* Coupling strength */
 63	PetscReal gamma = 1;
 64
 65	/* Spectral shift improves the convergence significantly */
 66	PetscReal shift = 100.;
 67
 68	/* Tolerance for the solver */
 69	PetscReal tol = 1e-8;
 70
 71	PetscOptionsGetInt(NULL, NULL, "-L", &L, NULL);
 72	PetscOptionsGetReal(NULL, NULL, "-Jz", &Jz, NULL);
 73	PetscOptionsGetInt(NULL, NULL, "-nev", &nev, NULL);
 74	PetscOptionsGetInt(NULL, NULL, "-ncv", &ncv, NULL);
 75	PetscOptionsGetReal(NULL, NULL, "-shift", &shift, NULL);
 76	PetscOptionsGetReal(NULL, NULL, "-tol", &tol, NULL);
 77	PetscOptionsGetReal(NULL, NULL, "-gamma", &gamma, NULL);
 78
 79	/* Number of MPI ranks */
 80	PetscMPIInt world_size;
 81
 82	/* Rank number */
 83	PetscMPIInt myrank;
 84
 85	PetscCall(SlepcInitialize(&argc, &argv, (char *)0, (char *)0));
 86	PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &world_size));
 87	PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &myrank));
 88
 89	if(myrank == 0){
 90		std::cout << "[myrank: " << myrank << "] - [main] - MPI processes: world_size = " << world_size << std::endl;
 91		std::cout << "[myrank: " << myrank << "] - [main] - system size: L = " << L << std::endl;
 92		std::cout << "[myrank: " << myrank << "] - [main] - interaction strength: Jz = " << Jz << std::endl;
 93		std::cout << "[myrank: " << myrank << "] - [main] - coupling strength: gamma = " << gamma << std::endl;
 94		std::cout << "[myrank: " << myrank << "] - [main] - number of requested eigen values: nev = " << nev << std::endl;
 95		std::cout << "[myrank: " << myrank << "] - [main] - maximal subspace dimension: ncv = " << ncv << std::endl;
 96		std::cout << "[myrank: " << myrank << "] - [main] - spectral shift to improve convergence: shift = " << shift << std::endl;
 97		std::cout << "[myrank: " << myrank << "] - [main] - tolerance: tol = " << tol << std::endl;
 98	}
 99
100
101	/* Operator class */
102	Lindbladian_NoU1 Lind(L); 
103
104	/* Simply XXZ-Heisenberg model with periodic boundary conditions */
105	for(uint64_t i = 0; i < L; i++){
106		Lind.add_operator(.5, {i,(i+1)%L}, {"S+","S-"});
107		Lind.add_operator(.5, {i,(i+1)%L}, {"S-","S+"});
108		Lind.add_operator(Jz, {i,(i+1)%L}, {"Sz","Sz"});
109	}
110
111	/* Non-hermitian jump operators */
112	Lind.add_jump_operator(gamma, {0}, {"S+"});
113	Lind.add_jump_operator(gamma, {0}, {"S-"});
114
115
116	if(myrank == 0){
117		Lind.info();
118	}
119
120	/* A simple spectral shift improves the convergence significantly */
121	Lind.add_diagonal(std::complex<double>(shift,0.));
122
123	/* Creating the matrix shell for matrix-free for multiplication */
124	auto Lind_matrix = Lind.create_PetscShellMatrix();
125	
126	/* Setting up the eigen problem solver */
127	EPS eps;
128	PetscCall(EPSCreate(PETSC_COMM_WORLD, &eps));
129	PetscCall(EPSSetOperators(eps, Lind_matrix, NULL));
130	PetscCall(EPSSetProblemType(eps, EPS_NHEP));
131	PetscCall(EPSSetWhichEigenpairs(eps, EPS_LARGEST_REAL));
132	PetscCall(EPSSetType(eps,EPSKRYLOVSCHUR));
133	PetscCall(EPSSetDimensions(eps, nev, ncv, PETSC_DEFAULT));
134
135	/* Define tolerance, last argument refers to the maximal number of iterations */
136	PetscCall(EPSSetTolerances(eps,tol,10000));
137
138	/* Setting the monitoring function */
139	PetscCall(EPSMonitorSet(eps,&monitor,NULL,NULL));
140
141	/* Solve the system */
142	auto start = std::chrono::high_resolution_clock::now();
143	PetscCall(EPSSolve(eps));
144	auto end = std::chrono::high_resolution_clock::now();
145	double duration = std::chrono::duration_cast<std::chrono::seconds>(end-start).count();
146
147
148	/* Get information from the solver */
149	PetscInt nconv;
150	PetscCall(EPSGetConverged(eps, &nconv));
151	PetscInt iter_number;
152	EPSGetIterationNumber(eps,&iter_number);
153
154
155	if(myrank == 0){
156		for(uint64_t c_for = 0UL; c_for < static_cast<uint64_t>(nconv); c_for++){
157			PetscScalar E, E_imag;
158			PetscReal err;
159			EPSGetErrorEstimate(eps,static_cast<PetscInt>(c_for),&err);
160			PetscCall(EPSGetEigenvalue(eps,static_cast<PetscInt>(c_for), &E, &E_imag));
161
162			std::cout << "[myrank: " << myrank << "] - [main] - " << c_for << "th eigenvalue = " << E - shift << " with error = " << err << std::endl;
163		}	
164
165	}
166
167	if(myrank == 0){
168		std::cout << "[myrank: " << myrank << "] - [main] - MPI processes: world_size = " << world_size << std::endl;
169		std::cout << "[myrank: " << myrank << "] - [main] - system size: L = " << L << std::endl;
170		std::cout << "[myrank: " << myrank << "] - [main] - interaction strength: Jz = " << Jz << std::endl;
171		std::cout << "[myrank: " << myrank << "] - [main] - coupling strength: gamma = " << gamma << std::endl;
172		std::cout << "[myrank: " << myrank << "] - [main] - number of requested eigen values: nev = " << nev << std::endl;
173		std::cout << "[myrank: " << myrank << "] - [main] - maximal subspace dimension: ncv = " << ncv << std::endl;
174		std::cout << "[myrank: " << myrank << "] - [main] - spectral shift to improve convergence: shift = " << shift << std::endl;
175		std::cout << "[myrank: " << myrank << "] - [main] - tolerance: tol = " << tol << std::endl;
176		std::cout << "[myrank: " << myrank << "] - [main] - number of iteration = " << iter_number << std::endl;
177		std::cout << "[myrank: " << myrank << "] - [main] - time = " << duration << "s" <<  std::endl;
178		std::cout << "[myrank: " << myrank << "] - [main] - time per iteration = " << duration/iter_number << "s" <<  std::endl;
179	}
180
181
182	PetscCall(MatDestroy(&Lind_matrix));
183	PetscCall(PetscFinalize());
184
185    return 0;
186}