[firedrake] Memory Allocation Problem

19 Jul 2016

      Dear Firedrakers,

I'm doing a Monte Carlo simulation in which the samples are functions 
with random numbers at each node on my mesh. I then create another 
function whose values are taken from combinations of the random numbers 
of the first function. My problem is that the code always runs out of 
memory, after a certain number (7700 or so) iterations!

I wonder if anyone can help, as I don't know anyway to get around this 
problem!

Thanks very much,

Tom Bendall

Error message:
Traceback (most recent call last):
   File "gibbs_distribution_orthogonal.py", line 236, in <module>
     main()
   File "gibbs_distribution_orthogonal.py", line 185, in main
     Q0.dat.data[:] = Q0_data
   File "<decorator-gen-13>", line 2, in data
   File 
"/scratch/tmb15/firedrake/local/lib/python2.7/site-packages/pyop2/versioning.py", 
line 115, in modifies
     retval = method(self, *args, **kwargs)
   File 
"/scratch/tmb15/firedrake/local/lib/python2.7/site-packages/pyop2/base.py", 
line 1760, in data
     _trace.evaluate(set([self]), set([self]))
   File 
"/scratch/tmb15/firedrake/local/lib/python2.7/site-packages/pyop2/base.py", 
line 170, in evaluate
     comp._run()
   File 
"/scratch/tmb15/firedrake/local/lib/python2.7/site-packages/pyop2/base.py", 
line 4073, in _run
     return self.compute()
   File 
"/scratch/tmb15/firedrake/local/lib/python2.7/site-packages/pyop2/base.py", 
line 4114, in compute
     fun = self._jitmodule
   File 
"/scratch/tmb15/firedrake/local/lib/python2.7/site-packages/pyop2/utils.py", 
line 64, in __get__
     obj.__dict__[self.__name__] = result = self.fget(obj)
   File 
"/scratch/tmb15/firedrake/local/lib/python2.7/site-packages/pyop2/sequential.py", 
line 151, in _jitmodule
     direct=self.is_direct, iterate=self.iteration_region)
   File 
"/scratch/tmb15/firedrake/local/lib/python2.7/site-packages/pyop2/caching.py", 
line 203, in __new__
     obj = make_obj()
   File 
"/scratch/tmb15/firedrake/local/lib/python2.7/site-packages/pyop2/caching.py", 
line 193, in make_obj
     obj.__init__(*args, **kwargs)
   File 
"/scratch/tmb15/firedrake/local/lib/python2.7/site-packages/pyop2/host.py", 
line 700, in __init__
     self.compile()
   File 
"/scratch/tmb15/firedrake/local/lib/python2.7/site-packages/pyop2/host.py", 
line 803, in compile
     compiler=compiler.get('name'))
   File 
"/scratch/tmb15/firedrake/local/lib/python2.7/site-packages/pyop2/compilation.py", 
line 304, in load
     dll = compiler.get_so(src, extension)
   File 
"/scratch/tmb15/firedrake/local/lib/python2.7/site-packages/pyop2/compilation.py", 
line 199, in get_so
     return ctypes.CDLL(soname)
   File "/usr/lib/python2.7/ctypes/__init__.py", line 365, in __init__
     self._handle = _dlopen(self._name, mode)
OSError: /tmp/pyop2-cache-uid791676/14b1f4fc69cf8da1847643ac027a1af6.so: 
cannot apply additional memory protection after relocation: Cannot 
allocate memory

Code:
from firedrake import *
parameters["pyop2_options"]["log_level"] = "WARNING"
import numpy
from datetime import datetime

def main():

     startTime = datetime.now()

     # Variables to be changed
     n0 = 30
     Froude = 1.0
     Beta = 0.0
     setup = 'SV'
     datadumpfreq = 100
     MC = 1000000 # number of monte carlo samples

     # Set up size of mesh and grid-spacing
     dx0 = 1.0/n0
     deltaX_over_deltaT = 1

     # Mesh is declared for firedrake
     m = PeriodicIntervalMesh(n0, 1.0)
     mesh = ExtrudedMesh(m, layers=n0, layer_height=dx0)

     # Type of finite element method is declared
     Vdg = FunctionSpace(mesh,"DG",1)
     Vcg = FunctionSpace(mesh,"CG",1)

### COEFFICIENTS ###

     #Some coefficients
     F = Constant(Froude) #Rotational Froude number
     beta = Constant(Beta) #beta plane coefficient
     #Aiming for Courant number < 0.3, assume |u|=1
     # 0.3 > 1*dt/dx => dt < dx*0.3
     Dt = dx0 / deltaX_over_deltaT
     dt = Constant(Dt)

### INITIAL CONDITIONS ###

         # Declare initial state of relative PV

     if (setup == 'SV'):
         number_of_vortices = 1
         vortex_centre = []
         vortex_coord_string = []
         vortex_string = []
         R = "(1./15)" # Radius of vortex
         pI = "(3.14159265359)"
         normalising_constant = 8.0
         BC = 0.0
         # assign vortex centres
         vortex_centre.append([0.5, 0.5])
         for i in range(number_of_vortices):
vortex_coord_string.append("pow(pow(x[0]-"+str(vortex_centre[i][0])+",2)+pow(x[1]-"+str(vortex_centre[i][1])+",2),0.5)")
         for i in range(number_of_vortices):
vortex_string.append("exp(-pow(("+str(vortex_coord_string[i])+")/"+str(R)+",4))")
         # making string for function
         function_string = 
str(normalising_constant)+"*"+str(vortex_string[0])
     elif (setup == 'UF'):
         U = 0.2
         BC = 0.2
         function_string = str(U)+"*x[1]"
     elif (setup == 'BI'):
         U = 0.2
         BC = 0.0
         function_string = 
str(U)+"*x[1]*(x[1]-1)-2.0*"+str(U)+"/"+str(Froude)
     else:
         print "ERROR: SET-UP NOT CORRECT"

     Q_initial = Function(Vdg).interpolate(Expression(function_string))

     Z0 = assemble(0.5 * Q_initial * Q_initial * dx)
     P0 = assemble(Q_initial * dx)

     print "Z0 = ", Z0
     print "P0 = ", P0

     Basis = Function(Vdg)
     Basis.dat.data[:] = 0

     A = n0 * (3.0 ** 0.5)
     B = 4.0 * (n0 ** 2.0)

     f = Function(Vdg).interpolate(Expression(str(Beta)+"*x[1]"))

     sigma = (2 * Z0 / B - (P0 / A) ** 2.0 ) ** 0.5

### MAKE EMPTY FUNCTIONS AND LISTS ###

     # Declare functions
     Q0 = Function(Vdg) # The actual field in the finite element basis

     Q1 = Function(Vdg)
     Qmean = Function(Vdg)

     q0 = Function(Vdg).assign(Q0 - f)

     psi0 = Function(Vcg)
     psi1 = Function(Vcg)
     psimean = Function(Vcg)

     # Declare lists to store observables
     energy = []
     total_p_vorticity = []
     total_p_enstrophy = []

     # Put initial state of model into files
     qfile = File("monte_carlo_q_"+str(setup)+"_diss_ORTH_q.pvd")
     psifile = File("monte_carlo_psi_"+str(setup)+"_diss_ORTH_psi.pvd")

### DESCRIBE TIME THINGS ###

     # Assign how the model should step through time
     m = 0
     m_list = []
     videodump = 0
     datadump = 0

### PV INVERSION PROBLEM ###

     # set up PV inversion. This is the form of the psi-q equation to be 
solved by firedrake
     psi = TrialFunction(Vcg)
     phi = TestFunction(Vcg)

     Apsi = (F * psi * phi + inner(grad(psi),grad(phi)) ) * dx
     Lpsi = - (Q1 * phi - f * phi) * dx

     bc1 = [DirichletBC(Vcg, (-BC / Froude + Beta / Froude), 
"top"),DirichletBC(Vcg, 0.0, "bottom")]

     psi_problem = LinearVariationalProblem(Apsi, Lpsi, psi0, bcs = bc1)
     psi_solver = LinearVariationalSolver(psi_problem,
                                             solver_parameters={
             'ksp_type':'cg',
             'pc_type':'bjacobi'
             })

     Q0.dat.data[:] = 0
     Qmean.dat.data[:] = 0
     psimean.dat.data[:] = 0

     Q0_data = [] # A list to contain the coefficients in the finite 
element basis

     for i in range(len(Q0.dat.data)):
         Q0_data.append(0)

     E = 0
     P = 0
     Z = 0

     Q_orth = [] # A list to contain the coefficients in my orthogonal basis
     for i in range(4):
         Q_orth.append(0)

### OPEN DATA FILE ###

     with open("monte_carlo_"+str(setup)+"_diss_ORTH.dat", "w") as out_file:

         out_string = "time" + " " + "energy" + " " + 
"potential_vorticity" + " " + "potential_enstrophy" + "\n"
         out_file.write(out_string)

### RUN PROGRAM ###

         while(m < (MC - 1/2)):

             j = 0
             for k in range(len(Q0_data)/4):

                 # generate the random numbers
                 for i in range(4):
                     Q_orth[i] = sigma * numpy.random.randn() + P0 / A

                 # calculate the coefficients in the finite element basis
                 Q0_data[j] = n0 * (2 * (3 ** 0.5) * Q_orth[0] - (2 ** 
0.5) * Q_orth[2] - Q_orth[3])
                 Q0_data[j+1] = n0 * (2*(2 ** 0.5) * Q_orth[2] - Q_orth[3])
                 Q0_data[j+2] = 3 * n0 * Q_orth[3]
                 Q0_data[j+3] = n0 * (2 * (3 ** 0.5) * Q_orth[1] - (2 ** 
0.5) * Q_orth[2] - Q_orth[3])
                 j += 4

             # Transfer coefficients back to function
             Q0.dat.data[:] = Q0_data

             Qmean.assign((m * Qmean) + Q0)
             Qmean.assign(Qmean / (m + 1))

             Q1.assign(Q0)
             psi_solver.solve()

             psimean.assign((m * psimean) + psi0)
             psimean.assign(psimean / (m + 1))

             this_E = assemble(-F * psi0 * Q0 * dx)

             E = m * E + this_E
             E = E / (m + 1)

             this_P = assemble(Q0 * dx)
             this_Z = assemble(0.5 * Q0 * Q0 * dx)

             P = m * P + this_P
             P = P / (m + 1)

             Z = m * Z + this_Z
             Z = Z / (m + 1)

             m += 1

             # Add observables to lists
             datadump += 1
             if(datadump == datadumpfreq):
                 datadump -= datadumpfreq
                 energy.append(E)
                 total_p_vorticity.append(P)
                 total_p_enstrophy.append(Z)
                 m_list.append(m)
                 print m_list[-1]
                 print datetime.now() - startTime

                 out_string = str(m_list[-1]) + " " + str(energy[-1]) + 
" " + str(total_p_vorticity[-1]) + " " + str(total_p_enstrophy[-1]) + "\n"
                 out_file.write(out_string)

     qfile << Qmean
     psifile << psimean

### Run the function main defined in this file ###
if __name__ == "__main__":
    main()

[firedrake] Memory Allocation Problem

T. M. Bendall