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()