# Memory leak using trilinos

```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)
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):

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 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 ]
```