Hello,

I am trying to run what should be a fairly simple simulation. My code is given below:

```
from netgen.csg import *
# Defining solid
sphere=Sphere(Pnt(0,0,0),0.55).bc('FSI')
pX=Plane(Pnt(0.5,0,0),Vec(1,0,0))
nX=Plane(Pnt(-0.5,0,0),Vec(-1,0,0))
pY=Plane(Pnt(0,0.5,0),Vec(0,1,0))
nY=Plane(Pnt(0,-0.5,0),Vec(0,-1,0))
pZ=Plane(Pnt(0,0,0.5),Vec(0,0,1))
nZ=Plane(Pnt(0,0,-0.5),Vec(0,0,-1))
Zielinski=sphere*pX*nX*pY*nY*pZ*nZ
# Making geometry from solid
geo=CSGeometry()
geo.Add(Zielinski.mat("domain"), bcmod=[(pX,"xplus"),(nX,"xminus"),(pY,"yplus"),(nY,"yminus"),(pZ,"zplus"),(nZ,"zminus")])
# Marking periodic surfaces
geo.PeriodicSurfaces(pX,nX)
geo.PeriodicSurfaces(pY,nY)
geo.PeriodicSurfaces(pZ,nZ)
# Meshing
from ngsolve.comp import Mesh
netgenMesh=geo.GenerateMesh(maxh=0.05)
mesh = Mesh(netgenMesh) # Converting netgen mesh to ngsolve mesh
Draw (mesh)
k = 3
velocityFiniteElementSpace = Periodic(VectorH1(mesh, order=k, dirichlet="FSI")) #discrete velocity space
pressureFiniteElementSpace = Periodic(H1(mesh, order=(k-1), dirichlet="FSI")) #discrete pressure space
productSpace = FESpace([velocityFiniteElementSpace,pressureFiniteElementSpace])
(trialFunction1,trialFunction2), (testFunction1,testFunction2) = productSpace.TnT()
e = CoefficientFunction( (1,0,0) ) # conservative source flux in COMSOL
LHS = BilinearForm(productSpace)
LHS += (InnerProduct(grad(trialFunction1),grad(testFunction1))-div(trialFunction1)*testFunction2-div(testFunction1)*trialFunction2)*dx
LHS.Assemble()
RHS = LinearForm(productSpace)
RHS += SymbolicLFI(e*testFunction1)
RHS.Assemble()
gfu = GridFunction(productSpace)
inv = LHS.mat.Inverse(freedofs=productSpace.FreeDofs(), inverse="umfpack")
res = RHS.vec.CreateVector()
res.data = RHS.vec - LHS.mat*gfu.vec
gfu.vec.data += inv * res
Redraw()
velocity = CoefficientFunction(gfu.components[0:2][0])
```

This is a modified version of the Stokes part of the Navier-Stokes tutorial (3.2 Incompressible Navier-Stokes equations — NGS-Py 6.2.2302 documentation). I am not trying to simulate a flow - this is a corollary for an acoustic problem.

The script runs, but if I increase the order (k) beyond 3 or if I decrease the element size any further, my workstation (which has 96GB RAM) runs out of memory and NGSolve closes with the message “Killed”. I am able to run the same simulation with the same geometry on OpenFOAM without it requiring nearly as much RAM. I know that OpenFOAM uses FVM and so this is an apples to oranges comparison, but surely, the simulation shouldn’t be taking so much RAM. Am I doing something wrong?