examples/SparseMatrix_petsc_Lindbladian

examples/SparseMatrix_petsc_Lindbladian/src/main.cpp

  1#include <iostream>
  2#include <mpi.h>
  3
  4/* slepcmfn includes matrix exponential */
  5#include <slepc.h> 
  6#include <slepcconf.h>
  7
  8/* Including Eigen */
  9#include <Eigen/Dense>
 10#include <Eigen/Eigenvalues>
 11
 12/* Maximal system size */
 13#define MaxSites 64
 14
 15/* Hilbert space dimension */
 16#define Q 2
 17
 18/* ScalarType */
 19#define ScalarType std::complex<double>
 20
 21#include "DanceQ.h"
 22using namespace danceq;
 23
 24int main(int argc, char *argv[]) {
 25
 26	/* Initializes Petsc and MPI */
 27	PetscInitialize(&argc, &argv, NULL, NULL);
 28
 29	/* System size */
 30	PetscInt L = 8;
 31
 32	/* Number of steps */
 33	PetscInt steps = 10;
 34
 35	/* Size of step */
 36	PetscReal delta_t = 0.5;
 37
 38	/* Field strength */
 39	PetscReal gamma = 0.5;
 40
 41	/* Input using Petsc */
 42	PetscOptionsGetInt(NULL, NULL, "-L", &L, NULL);
 43	PetscOptionsGetInt(NULL, NULL, "-steps", &steps, NULL);
 44	PetscOptionsGetReal(NULL, NULL, "-delta_t", &delta_t, NULL);
 45	PetscOptionsGetReal(NULL, NULL, "-gamma", &gamma, NULL);
 46
 47	/* Number of MPI ranks */
 48	PetscMPIInt world_size;
 49
 50	/* Rank number */
 51	PetscMPIInt myrank;
 52
 53	PetscCall(SlepcInitialize(&argc, &argv, (char *)0, (char *)0));
 54	PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &world_size));
 55	PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &myrank));
 56
 57
 58	/* Creating fully imbalance density matrix rho_init = |Q-1;0;Q-1;0;...><Q-1;0;Q-1;0;...> as Petsc Vec */
 59	State rho_init;
 60	rho_init.set_range(State::template_struct::template_Q-1,0UL,2UL*L);
 61	for(uint64_t i_for = 0UL; i_for < L; i_for++){
 62		if(i_for % 2 == 1){
 63			rho_init.set_state(0,i_for);
 64			rho_init.set_state(0,i_for+L);
 65		}
 66	}
 67
 68	/* Retrieve particle number */
 69	uint64_t n = rho_init.get_particle_number(0,L);
 70
 71	if(myrank == 0){
 72		std::cout << "[myrank: " << myrank << "] - [main] - world_size = " << world_size << std::endl;
 73		std::cout << "[myrank: " << myrank << "] - [main] - delta_t = " << delta_t << std::endl;
 74		std::cout << "[myrank: " << myrank << "] - [main] - steps = " << steps << std::endl;
 75		std::cout << "[myrank: " << myrank << "] - [main] - L = " << L << std::endl;
 76		std::cout << "[myrank: " << myrank << "] - [main] - n = " << n << std::endl;
 77		std::cout << "[myrank: " << myrank << "] - [main] - gamma = " << gamma << std::endl;
 78		std::cout << "[myrank: " << myrank << "] - [main] - initial density matrix: " << rho_init.get_string(0UL,L) << " " << rho_init.get_string(L,2UL*L) << std::endl;
 79	}
 80
 81	/* Lindbladian */
 82	Lindbladian_U1 Lind(L,n); 
 83
 84	/* Defining the SU(2) symmetry Heisenberg model with PBC */
 85	double J = 1.;
 86	for(uint64_t i = 0; i < L; i++){
 87		Lind.add_operator(J, {i,(i+1)%L}, {"Sx","Sx"});
 88		Lind.add_operator(J, {i,(i+1)%L}, {"Sy","Sy"});
 89		Lind.add_operator(J, {i,(i+1)%L}, {"Sz","Sz"});
 90		Lind.add_jump_operator(gamma, {i}, {"Sz"});
 91	}
 92
 93	/* Hamiltonian class for the observable */
 94	Hamiltonian_U1 obs(L,n); 
 95
 96	/* Defining the observable for measurements */
 97	for(uint64_t i = 0; i < L; i++){
 98		obs.add_operator(std::pow(-1.,i), {L-i-1}, {"Sz"});
 99	}
100
101	uint64_t dim = Lind.get_dim();
102	if(myrank == 0){
103		std::cout << "[myrank: " << myrank << "] - [main] - dim = " << dim << std::endl;
104		Lind.info();
105	}
106
107	/* Generating the matrices */
108	auto Lind_matrix = Lind.create_PetscSparseMatrix();
109	if(dim < 10){
110		PetscCall(MatView(Lind_matrix, PETSC_VIEWER_STDOUT_WORLD));
111	}
112
113
114	/* Retrieving the index of rho_init to create the initial vector */
115	PetscInt imb_index = static_cast<PetscInt>(Lind.get_index(rho_init));
116	if(myrank == 0){
117		std::cout << "[myrank: " << myrank << "] - [main] - rho_init = " << rho_init.get_string(0UL,2UL*L) << " with global index in the vectorized density matrix: " << imb_index << std::endl;
118	}
119
120	PetscScalar one = 1.;
121	Vec rho;
122	PetscCall(MatCreateVecs(Lind_matrix, &rho, PETSC_NULLPTR));
123	PetscCall(VecZeroEntries(rho));
124	PetscCall(VecSetValues(rho, 1, &imb_index, &one, INSERT_VALUES));
125	PetscCall(VecAssemblyBegin(rho));
126	PetscCall(VecAssemblyEnd(rho));
127	
128
129	/* Second Vec for computation */
130	Vec rhoPrime;
131	PetscCall(VecDuplicate(rho, &rhoPrime));
132
133	/* Setting up the matrix exponential with delta_t */
134	MFN mfn;
135	FN f;
136	PetscCall(MFNCreate(PETSC_COMM_WORLD,&mfn));
137	PetscCall(MFNSetOperator(mfn,Lind_matrix));
138	PetscCall(MFNGetFN(mfn,&f));
139	PetscCall(FNSetType(f,FNEXP));
140	PetscCall(FNSetScale(f,delta_t,1.0)); /* Setting the delta_t steps */
141	PetscCall(MFNSetTolerances(mfn,1e-8,PETSC_DEFAULT));
142	PetscCall(MFNSetFromOptions(mfn));
143
144	/* Loop of steps */
145	double t = 0.;
146	for(uint64_t i = 0; i < steps; i++){
147		t += delta_t;
148
149		/* Time evolution rhoPrime = e^{delta_t L} rho */
150		PetscCall(MFNSolve(mfn,rho,rhoPrime));
151
152		/* Setting rho = rhoPrime */
153		PetscCall(VecCopy(rhoPrime, rho));
154
155		/* Retrieving the density matrix as a Eigen matrix on rank 0 */
156		double S = von_Neumann_entropy_from_PetscVec(rho,0);
157
158		/* Retrieving the imbalance from the operator obs */
159		PetscScalar imbalance = obs.trace_over_density_matrix(rho);
160
161		/* Output */
162		if(myrank == 0){
163			std::cout << "[myrank: " << myrank << "] - [main] - t = " << t << ": S/Smax = " << S/std::log(dim)*2. << ", - imbalance = " << imbalance << std::endl;
164		}
165
166	}
167
168	/* Deleting Petsc objects */
169	PetscCall(MatDestroy(&Lind_matrix));
170	PetscCall(VecDestroy(&rho));
171	PetscCall(VecDestroy(&rhoPrime));
172
173	/* Cleaning up */
174	PetscCall(PetscFinalize());
175
176    return 0;
177}