Hello,I still have a large memory leak when using Trilinos. I am not sure where to start looking so I made an example code that produces my problem in hopes that someone can help me.
But! my example is cool. I implemented *Density Functional Theory in FiPy!***My code is slow, but runs in parallel and is simple (relative to most DFT codes). The example I have attached is just a lithium and hydrogen atom. The electrostatic boundary conditions are goofy but work well enough for demonstration purposes. If you set use_trilinos to True, the code will slowly use more memory. If not, it will try to use Pysparse.
Thanks, -Michael Waters
input.xyz
Description: Xmol XYZ data
#License: if you do something cool with this, please cite me, Michael Waters. lz = 3.0 # Ang lx = 4.0 ly = 3.0 dx = 0.03 dy = 0.03 dz = 0.03 # Ang nz = int(lz/dz) ny = int(ly/dy) nx = int(lx/dx) lz = nz*dz ly = ny*dy lx = nx*dx from mpi4py import MPI m4comm = MPI.COMM_WORLD if m4comm.Get_rank() == 0: print 'nz', nz ,'dz', dz, 'lz', lz print 'ny', ny ,'dy', dy, 'ly', ly print 'nx', nx ,'dx', dx, 'lx', lx number_of_electrons = 4 steps = 100 accuracy = 10.0**-5 phi_solver_iterations_per_step = 1000 save_infrequency = 5 initial_solver_iterations_per_step = 7 ### read_restart = True GC = True # does garbage collection at each self-consistent iteration use_trilinos = True ######################## Constants! qe = -1.0 # electron charge a0 = 0.52917721067 # codata 2014 in Ang hbar = 1.054571800 * 10**-34 #(SI) me0 = 9.10938356 * 10**-31 #(SI) eV = 1.6021766208 * 10**-19 #(SI) tcoeff = (hbar**2)/(2.0*me0)*(1.0/eV)*(10.0**10)**2 # = hbar^2/(2*Me) in units of [eV*Ang^2] from numpy import pi eu = 1.6021766208 * 10**-19 #electron charge (SI) epsilon0 = (8.854187817 * 10**-12) # in SI epsilon0 = epsilon0 * (1.0/eu)*(1.0/10.0**10) # conversion cc0 = 1.0/(4.0*pi*epsilon0) # about 14.4 eV*Ang usually used in atomistics ###########3 from fipy import __version__ as FiPy_version if m4comm.Get_rank() == 0: print 'FiPy Version:', FiPy_version ################ from fipy import * from numpy.random import rand from numpy import sqrt, savetxt, array, exp, arange, zeros, loadtxt, sign import os.path if GC: import gc if m4comm.Get_rank() == 0: print 'Creating mesh and cell variables' mesh = Grid3D(dx=dx, dy=dy, dz=dz, nx=nx, ny=ny, nz=nz,)# origin = (lx/2.0, ly/2.0, lz/2.0)) Vcores = CellVariable(name = 'Core Potential for Electrons', mesh=mesh, value = 0.0) rho_e = CellVariable(name = 'Electron Density', mesh=mesh, value = 0.0) phi_e = CellVariable(name = 'Hartree Potential', mesh=mesh, value = 0.0) ### choosing a solver if use_trilinos: from fipy.solvers.trilinos import LinearGMRESSolver as phi_method from fipy.solvers.trilinos.preconditioners import MultilevelSGSPreconditioner as phi_precon #from fipy.solvers import trilinos as trisol #print dir(trisol) #from fipy.solvers.trilinos import DefaultAsymmetricSolver as psi_method # this is just GMRES... #from fipy.solvers.trilinos import DefaultSolver as psi_method #from fipy.solvers.trilinos import DummySolver as psi_method #from fipy.solvers.trilinos import GeneralSolver as psi_method from fipy.solvers.trilinos import LinearBicgstabSolver as psi_method # works well with Jacobi, 20 iterations #from fipy.solvers.trilinos import LinearCGSSolver as psi_method # DNF with Jacobi, might need reduced minimum mixing paramter #from fipy.solvers.trilinos import LinearGMRESSolver as psi_method # works well with Jacobi, 32 iterations, optimal mixer might need to change iterations faster for it to keep up with the needs of the GMRES solver, looks like it could have done it in 21 #from fipy.solvers.trilinos import LinearLUSolver as psi_method # DNF with Jacobi #from fipy.solvers.trilinos import LinearPCGSolver as psi_method # standard, 23 iterations #from fipy.solvers.trilinos import TrilinosMLTest as psi_method # doesn't seem to be anything ##### choosing a preconditioner #from fipy.solvers.trilinos import preconditioners as precons #print dir(precons) #from fipy.solvers.trilinos.preconditioners import DomDecompPreconditioner as psi_precon # NotImplementedError #from fipy.solvers.trilinos.preconditioners import ICPreconditioner as psi_precon # segfaults... with permisions errors from fipy.solvers.trilinos.preconditioners import JacobiPreconditioner as psi_precon # standard. #from fipy.solvers.trilinos.preconditioners import MultilevelDDMLPreconditioner as psi_precon # bizzare warning about compilation..., DNF with LinearBicgstab #from fipy.solvers.trilinos.preconditioners import MultilevelDDPreconditioner as psi_precon # unstable #from fipy.solvers.trilinos.preconditioners import MultilevelNSSAPreconditioner as psi_precon # unstable #from fipy.solvers.trilinos.preconditioners import MultilevelSAPreconditioner as psi_precon # unstable #from fipy.solvers.trilinos.preconditioners import MultilevelSGSPreconditioner as psi_precon # unstable #from fipy.solvers.trilinos.preconditioners import MultilevelSolverSmootherPreconditioner as psi_precon # Segfaults... else: from fipy.solvers.pysparse import LinearPCGSolver as psi_method from fipy.solvers.pysparse import LinearPCGSolver as phi_method from fipy.solvers.pysparse import SsorPreconditioner as psi_precon from fipy.solvers.pysparse import SsorPreconditioner as phi_precon psi_solver = psi_method(iterations = initial_solver_iterations_per_step, precon = psi_precon() ) phi_solver = phi_method(iterations = phi_solver_iterations_per_step, precon = phi_precon() ) xc, yc, zc = mesh.cellCenters dv = dx*dy*dz simulation_volume = lx*ly*lz if m4comm.Get_rank() == 0: print 'simulation volume', simulation_volume ############## set initial applied potential if m4comm.Get_rank() == 0: print "Setting Core Potentials" atomic_data=loadtxt('input.xyz') for i in range(atomic_data.shape[0]): Z = atomic_data[i][0] x = lx/2.0 + atomic_data[i][1] y = ly/2.0 + atomic_data[i][2] z = lz/2.0 + atomic_data[i][3] #print x coeff = -Z/(4.0*pi*epsilon0) Vcores.setValue(coeff/sqrt( (xc-x)**2 + (yc-y)**2 + (zc-z)**2) + Vcores) if True: if m4comm.Get_rank() == 0: print 'saving core potentials...' savetxt('Vcores.txt', Vcores.globalValue) ########### write cell data fid = open('cell_data.txt','w') fid.write('numcells and dx in x direction\n') fid.write('%i\t%f\n'%(nx,dx)) fid.write('numcells and dy in y direction\n') fid.write('%i\t%f\n'%(ny,dy)) fid.write('numcells and dz in z direction\n') fid.write('%i\t%f\n'%(nz,dz)) fid.close() ################## Initialize all psi's psi_e_list = [] for n in range(number_of_electrons): psi_e_list.append( CellVariable(name = '$\psi_{e%i}$' %n , mesh=mesh, value = 1.0, hasOld = True)) ##### setting up charge density for n in range( number_of_electrons): rho_e = rho_e + qe*psi_e_list[n]*psi_e_list[n] ############## this sets up the electrostatic potential phi_e.equation = (0.0 == DiffusionTerm(coeff = epsilon0) + rho_e) ##################### some function definitions def inner_product(psi_a,psi_b): ans = ((psi_a * psi_b).cellVolumeAverage * simulation_volume).value return ans def normalize(psi): a = 1.0/sqrt( inner_product(psi,psi) ) # normalizing factor# #print 'normalizing factor',a psi.setValue( psi*a) def num_electrons(psi): return inner_product(psi,psi) def sandwich_H(psi_a, psi_b): # the hamiltonian is entered twice in this simulation, once as the equation to be solved and here as a function to be evaluated Hpsi = -tcoeff * psi_b.faceGrad.divergence + (Vcores + qe*phi_e)*psi_b return ((psi_a*Hpsi).cellVolumeAverage * simulation_volume).value def H_expect(psi): return sandwich_H(psi, psi) def subtract_projected_on_normed(psi_b,psi_a): coef=inner_product(psi_a,psi_b) psi_b.setValue(psi_b - psi_a*coef) def kinetic_energy_e(psi): Tpsi = -tcoeff* psi.faceGrad.divergence return ((psi*Tpsi).cellVolumeAverage * simulation_volume) ##### initialization of wave functions ############## if os.path.isfile('phi_e.txt') == False: read_restart = False if read_restart: if m4comm.Get_rank() == 0: print 'Reading restart' ## electrons for n in range(number_of_electrons): if os.path.isfile('psi_e_%i.txt'%n): if m4comm.Get_rank() == 0: print 'Reading electron wave function %i'%n input_data = loadtxt('psi_e_%i.txt'%n) psi_e_list[n].setValue(input_data) psi_e_list[n].updateOld() else: # does random initialization if m4comm.Get_rank() == 0: print 'Random initialization for missing electron wave functions %i'%n psi_e_list[n].setValue(rand(nx*ny*nz)-0.5) ## phi_e if m4comm.Get_rank() == 0: print 'Reading electron potential...' input_data = loadtxt('phi_e.txt') phi_e.setValue(input_data) else: if m4comm.Get_rank() == 0: print 'Guess initialization' # for all psi's for n in range(number_of_electrons): psi_e_list[n].setValue(rand(nx*ny*nz)-0.5) phi_e.setValue(rand(nx*ny*nz)) ############## orthonormalize functions then normalize if m4comm.Get_rank() == 0: print 'Orthonormalizing...' if number_of_electrons > 0: normalize(psi_e_list[0]) #normalize the lowest energy electron if number_of_electrons >1 :normalize(psi_e_list[1]) # only if there is second electron for n in range(2, number_of_electrons): for lower_wave in range(n%2,n,2): #since range goes from 0 to n-1, skipping opposite-spin electrons subtract_projected_on_normed(psi_e_list[n],psi_e_list[lower_wave]) normalize(psi_e_list[n]) ### checking normalization if False: for n in range(number_of_electrons): num = num_electrons(psi_e_list[n]) if m4comm.Get_rank() == 0: print n, num ####### initialize boundary conditions from fipy.boundaryConditions.constraint import Constraint if m4comm.Get_rank() == 0: print "applying BC's" constraint_list_e = [] for n in range(number_of_electrons): constraint_list_e.append( Constraint(0, mesh.exteriorFaces)) psi_e_list[n].constrain(constraint_list_e[n]) phi_e.constrain(Constraint(0.0, mesh.facesBack)) phi_e.constrain(Constraint(0.0, mesh.facesFront)) phi_e.faceGrad.constrain( 0.0, mesh.facesRight) phi_e.faceGrad.constrain( 0.0, mesh.facesLeft) phi_e.faceGrad.constrain( 0.0, mesh.facesTop) phi_e.faceGrad.constrain( 0.0, mesh.facesBottom) ############ initialize energies ... <psi_n|H|psi_n> = E_n energy_list_e = [] solver_iterations_per_step_e = [] rel_error_e = [] for n in range(number_of_electrons): e = H_expect(psi_e_list[n]) energy_list_e.append(Variable(value=e) ) solver_iterations_per_step_e.append(initial_solver_iterations_per_step) rel_error_e.append(1.0) Te = 0.0 # total kinetic energy for n in range(number_of_electrons): Te += kinetic_energy_e(psi_e_list[n]) toten_old = Te - inner_product(rho_e,Vcores) # total energy, rho is charge density and is always negative because electrons if m4comm.Get_rank() == 0: print '############### Starting Self-Consistent Iteration ###########' if m4comm.Get_rank() == 0: print 'Initial Total Electron Energy', toten_old if read_restart: mode = 'a' else: mode = 'w' fid = open('iteration_energies.log', mode) fid_rel_error = open('relative_errors.log', mode) #### HAMILTONIAN and Schroedinger equation! for n in range(0, number_of_electrons): psi_e_list[n].equation = (0 == DiffusionTerm(-tcoeff) + (Vcores + qe*phi_e - energy_list_e[n]) * psi_e_list[n] ) max_rel_error = 1.0 step = 1 while step<=steps and max_rel_error > accuracy: if m4comm.Get_rank() == 0: print 'Solving Phi e' phi_e_res = phi_e.equation.sweep(var = phi_e, solver = phi_solver ) if m4comm.Get_rank() == 0: print 'phi_e_res' , phi_e_res ##############Solving the ELECTRON wavefunctions with orthonormalization after each solve for n in range(0, number_of_electrons): if rel_error_e[n] < accuracy: if m4comm.Get_rank() == 0: print '--- Skipping electron %i because it is converged' %n else: psi_e_list[n].updateOld() ## update before solving for the next wavefunctions if m4comm.Get_rank() == 0: print '--- Working on electron %i' %n if m4comm.Get_rank() == 0: print 'Solving...' #### Solving it psi_solver.iterations = solver_iterations_per_step_e[n] psi_e_res = psi_e_list[n].equation.sweep(var = psi_e_list[n], solver = psi_solver ) if m4comm.Get_rank() == 0: print 'psi_e_res', psi_e_res ############## if m4comm.Get_rank() == 0: print 'Orthonormalizing...' normalize(psi_e_list[n]) for lower_wave in range(n%2,n,2): #only normalizes against the spin up or spin down electrons subtract_projected_on_normed(psi_e_list[n],psi_e_list[lower_wave]) normalize(psi_e_list[n]) Eold = float(energy_list_e[n].value) # evaluates the Variable Enew = H_expect(psi_e_list[n]) energy_list_e[n].setValue( Enew ) rel_error_e[n] = abs((Enew-Eold)/Enew) max_rel_error = rel_error_e[0] # finds the largest relative error for n in range(number_of_electrons): if rel_error_e[n] > max_rel_error: max_rel_error = rel_error_e[n] ########## Screen and log output, everything below here is not very iteresting. if m4comm.Get_rank() == 0: print '~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~' print '\nStep', step, Te = 0.0 for n in range(number_of_electrons): Te += kinetic_energy_e(psi_e_list[n]) Ecores = -inner_product(rho_e,Vcores) toten = Te + Ecores if m4comm.Get_rank() == 0: print 'Total Energy' , toten, de = toten - toten_old if m4comm.Get_rank() == 0: print 'dE', de, print 'dE/|toten|' , de/abs(toten), print 'max_rel_error', max_rel_error toten_old = toten if m4comm.Get_rank() == 0: print '\nElectron Energies', energy_list_e fid.write('\nStep ' + str(step)+ ' Total Energy ' + str(toten) +' dE '+ str(de) + ' dE/|toten| ' + str(de/abs(toten) ) + ' Max Rel Error ' +str(max_rel_error) )#+ '\nHole Energies '+ str(energy_list_h[0])) if m4comm.Get_rank() == 0: # this section saves file output def str_list(a): output = '' for thing in a: output = output + str(thing) + ' ' return output fid_rel_error.write(str_list(rel_error_e) + '\n') fid_rel_error.flush() ############# This output the energy eigenvalues if m4comm.Get_rank() == 0: e_energy_string = '' for n in range(number_of_electrons): e_energy_string = e_energy_string+' '+str(float(energy_list_e[n].value)) fid.write('\nElectron Energies ' + e_energy_string) ############# hartree_energy = 0.5*inner_product(rho_e,phi_e) if m4comm.Get_rank() == 0: fid.write('\nhartree_energy ' + str(hartree_energy) ) print 'hartree_energy ' + str(hartree_energy) fid.write('\nTe ' + str(Te) ) print 'Te ' + str(Te) fid.write('\nEcores ' + str(Ecores) ) print 'Ecores ' + str(Ecores) fid.flush() #### saving output### savetxt('Energy_eigenvalues_e.txt',array(energy_list_e)) if step%save_infrequency == 0: if m4comm.Get_rank() == 0: print 'Saving Cell data...' for n in range( number_of_electrons): savetxt('psi_e_%i.txt'%n,psi_e_list[n].globalValue) savetxt('phi_e.txt',phi_e.globalValue) savetxt('rho_e.txt',rho_e.globalValue) if GC: if m4comm.Get_rank() == 0: print 'Collecting Garbage...' gc.collect() step +=1 ## loop ends here! #### saving output### savetxt('Energy_eigenvalues_e.txt',array(energy_list_e)) if m4comm.Get_rank() == 0: print 'Finished!, Saving Cell data...' for n in range( number_of_electrons): savetxt('psi_e_%i.txt'%n,psi_e_list[n].globalValue) savetxt('phi_e.txt',phi_e.globalValue) savetxt('rho_e.txt',rho_e.globalValue)
_______________________________________________ fipy mailing list fipy@nist.gov http://www.ctcms.nist.gov/fipy [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]