examples/SparseMatrix_petsc_Hamiltonian_real_time_evolution

examples/SparseMatrix_petsc_Hamiltonian_real_time_evolution/src/main.cpp

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