HI All,
I think I run into some memory leak issues using trilinos as my solver both
for serial and parallel case. The total memory goes up from initial around
300 MB to the my system memory limit (like 29 GB), and program stops. I
have attached my example code below.
Somethings I am not sure if I implement it correctly:
1. in the example code, I basically solve multiple equations at the same
time. I use list, and append my CellVarible together
2. Similar, I use list and append my equation together
3. Create coupled equation by looping through the list element, and solve
this coupled equation.
4. During the solution step, I loop through the list element, and do
updateOld()
The results of the numerical calculation looks correct to me, I just have
this memory issue that prevent me from going longer time steps because I
run out of memory.
I wonder is there any workaround/fix to this problem ?
I am using the build from Anaconda repository, default trillions solver,
and my OS is Ubuntu 16.04.
from fipy import *
import numpy as np
import matplotlib as matplotlib_plot
import scipy.interpolate
import resource
def make2dList(rows, cols):
a=[]
for row in xrange(rows): a += [[0]*cols]
return a
mesh = Gmsh2D('''
lc = 1e-2;
Point(1) = {-1, 0, 0, lc};
Point(2) = {0, 0, 0, lc};
Point(3) = {1, 0, 0, lc};
Point(4) = {0,-0.1,0,lc};
Ellipse(5) = {1,2,2,4};
Ellipse(6) = {4,2,2,3};
Line(7) = {3,2};
Line(8) = {2,1};
Circle(9) = {1,2,3};
Line Loop(9) = {6,7, 8,5};
Line Loop(10) = {-6,-5,9};
Plane Surface (11) = {9};
Plane Surface (12) = {10};
Physical Surface("top") = {11};
Physical Surface("bottom") = {12};
Physical Line("interface") = {5, 6};
''' % locals())
Np = 3
phi =[];
for species_i in range(Np-1):
species = CellVariable(mesh=mesh, hasOld = True, value=1.0/Np)
phi.append(species)
xFace,yFace=mesh.faceCenters
Epsilon = 0.1
velocityX_Top = (1-Epsilon*Epsilon)/(Epsilon*Epsilon)*(yFace + Epsilon*np.sqrt(1-xFace*xFace))
velocityY_Top = (1-Epsilon*Epsilon)/Epsilon*(xFace*yFace)/(np.sqrt(1-xFace*xFace))
velocityX_Bottom = yFace
velocityY_Bottom = -xFace
velocityVector = FaceVariable(mesh=mesh, rank=1)
topFaces = numerix.unique(mesh.cellFaceIDs[..., mesh.physicalCells["top"].value].flatten().filled())
bottomFaces = numerix.unique(mesh.cellFaceIDs[..., mesh.physicalCells["bottom"].value].flatten().filled())
velocityVector[0, topFaces] = velocityX_Top[..., topFaces]
velocityVector[1, topFaces] = velocityY_Top[..., topFaces]
velocityVector[0, bottomFaces] = velocityX_Bottom[..., bottomFaces]
velocityVector[1, bottomFaces] = velocityY_Bottom[..., bottomFaces]
local_shear_rate = FaceVariable(mesh=mesh, rank=1)
local_shear_rate[0, topFaces] = 0.
local_shear_rate[1, topFaces] = 1.0
local_shear_rate[0, bottomFaces] = 0.
local_shear_rate[1, bottomFaces] = 0.
# Velocity Field Plot
#viewer_velocity = MatplotlibVectorViewer(vars=velocityVector,title="velocity_vector")
#viewer_velocity.plot()
Pe = 5.
D = FaceVariable(mesh=mesh, rank=1)
D[0,topFaces] = 0.
D[1,topFaces] = Epsilon*Epsilon/Pe
D[0,bottomFaces] = 0.
D[1,bottomFaces] = 0.
Lambda = make2dList(Np, Np)
radius = np.linspace(0.001,0.003,num=Np)
beta = 0.26
R = 0.075
c = (1-pow(Epsilon,2))/(R*pow(Epsilon,3))
for i in range(Np):
for j in range(Np):
Lambda[i][j] = c*beta*2.0*np.minimum(radius[i],radius[j])*np.log(radius[i]/radius[j])
coupledEqn = 0
eq = []
Segregation_Term = [0]*(Np-1)
Last_Species = 1.0
for species_i in range(Np-1):
print species_i
for species_j in range(Np -1):
Segregation_Term[species_i] = Segregation_Term[species_i] + Lambda[species_i][species_j]*phi[species_j]
Last_Species = Last_Species - phi[species_j]
Segregation_Term[species_i] = Segregation_Term[species_i] + Lambda[species_i][Np-1]*Last_Species
eq.append(TransientTerm(var=phi[species_i]) \
+ VanLeerConvectionTerm(coeff = velocityVector + Epsilon*local_shear_rate*Segregation_Term[species_i].faceValue, var=phi[species_i])\
== DiffusionTerm(coeff=D,var=phi[species_i]))
phi[species_i].faceGrad.constrain([1.0/Epsilon*local_shear_rate*Pe*(phi[species_i]*Segregation_Term[species_i]).faceValue], where = mesh.facesTop | mesh.physicalFaces["interface"]) #boundary conditions at the top and bottom
Last_Species = 1.0
for species_i in range(Np-1):
coupledEqn = coupledEqn & eq[species_i]
timeStepDuration = 1/2000.
steps = 1100
k = 0
for step in range(steps):
res = 1000
while res > 1e-6:
res = coupledEqn.sweep(dt=timeStepDuration)
print "MaxRSS:", resource.getrusage(resource.RUSAGE_SELF).ru_maxrss
print res
for species_i in range(Np-1):
phi[species_i].updateOld()
print step
_______________________________________________
fipy mailing list
fipy@nist.gov
http://www.ctcms.nist.gov/fipy
[ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
