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