In [1]:
import petsc4py
import sys
petsc4py.init(sys.argv)
from petsc4py import PETSc
from dolfin import *
Print = PETSc.Sys.Print
# from MatrixOperations import *
import numpy as np
import PETScIO as IO
import common
import scipy
import scipy.io
import time
import BiLinear as forms
import IterOperations as Iter
import MatrixOperations as MO
import CheckPetsc4py as CP
import ExactSol
import Solver as S
import NSprecondSetup
import memory_profiler
import gc
from scipy.sparse import coo_matrix, csr_matrix, spdiags, bmat, linalg
#@profile
m = 2
errL2u =np.zeros((m-1,1))
errH1u =np.zeros((m-1,1))
errL2p =np.zeros((m-1,1))
errL2b =np.zeros((m-1,1))
errCurlb =np.zeros((m-1,1))
errL2r =np.zeros((m-1,1))
errH1r =np.zeros((m-1,1))
l2uorder = np.zeros((m-1,1))
H1uorder =np.zeros((m-1,1))
l2porder = np.zeros((m-1,1))
l2border = np.zeros((m-1,1))
Curlborder =np.zeros((m-1,1))
l2rorder = np.zeros((m-1,1))
H1rorder = np.zeros((m-1,1))
NN = np.zeros((m-1,1))
DoF = np.zeros((m-1,1))
Velocitydim = np.zeros((m-1,1))
Magneticdim = np.zeros((m-1,1))
Pressuredim = np.zeros((m-1,1))
Lagrangedim = np.zeros((m-1,1))
Wdim = np.zeros((m-1,1))
iterations = np.zeros((m-1,1))
SolTime = np.zeros((m-1,1))
udiv = np.zeros((m-1,1))
MU = np.zeros((m-1,1))
level = np.zeros((m-1,1))
NSave = np.zeros((m-1,1))
Mave = np.zeros((m-1,1))
TotalTime = np.zeros((m-1,1))
nn = 2
dim = 2
ShowResultPlots = 'yes'
split = 'Linear'
parameters["form_compiler"]["no-evaluate_basis_derivatives"] = False
MU[0]= 1e0
for xx in xrange(1,m):
print xx
level[xx-1] = xx + 0
nn = 2**(level[xx-1])
# Create mesh and define function space
nn = int(nn)
NN[xx-1] = nn/2
# parameters["form_compiler"]["quadrature_degree"] = 6
# parameters = CP.ParameterSetup()
mesh = UnitSquareMesh(nn,nn)
# mesh = RectangleMesh(0,0,2*np.pi,2*np.pi,nn,nn)
order = 2
parameters['reorder_dofs_serial'] = False
Velocity = VectorFunctionSpace(mesh, "CG", order)
Pressure = FunctionSpace(mesh, "CG", order)
Magnetic = FunctionSpace(mesh, "N1curl", order)
Lagrange = FunctionSpace(mesh, "CG", order)
W = MixedFunctionSpace([Velocity, Pressure, Magnetic,Lagrange])
print "\n\nW: ",Wdim[xx-1],"Velocity: ",Velocitydim[xx-1],"Pressure: ",Pressuredim[xx-1],"Magnetic: ",Magneticdim[xx-1],"Lagrange: ",Lagrangedim[xx-1],"\n\n"
dim = [Velocity.dim(), Pressure.dim(), Magnetic.dim(), Lagrange.dim()]
def boundary(x, on_boundary):
return on_boundary
u0, p0,b0, r0, Laplacian, Advection, gradPres,CurlCurl, gradR, NS_Couple, M_Couple = ExactSol.MHD2D(4,1)
bcu = DirichletBC(Velocity,u0, boundary)
bcb = DirichletBC(Magnetic,Expression(('0','0')), boundary)
bcr = DirichletBC(Lagrange,Expression(('0')), boundary)
# bc = [u0,p0,b0,r0]
bcs = [bcu,bcb,bcr]
FSpaces = [Velocity,Pressure,Magnetic,Lagrange]
kappa = 10.0
Mu_m =10.0
MU = 1.0/1
IterType = 'Full'
Split = "No"
Saddle = "No"
Stokes = "No"
SetupType = 'python-class'
F_NS = -MU*Laplacian+Advection+gradPres-kappa*NS_Couple
if kappa == 0:
F_M = Mu_m*CurlCurl+gradR -kappa*M_Couple
else:
F_M = Mu_m*kappa*CurlCurl+gradR -kappa*M_Couple
params = [kappa,Mu_m,MU]
MO.PrintStr("Seting up initial guess matricies",2,"=","\n\n","\n")
#C = HiptmairMatrices[0]
#Px = HiptmairMatrices[1][0]
#Py = HiptmairMatrices[1][1]
In [26]:
row = np.zeros(6)
col = np.zeros(6)
dataX = np.zeros(6)
dataY = np.zeros(6)
E = Edge(mesh,0)
VertexNodes = E.entities(0)
V1 = Vertex(mesh,VertexNodes[0])
V2 = Vertex(mesh,VertexNodes[1])
Xtang = V2.x(0)-V1.x(0)
Ytang = V2.x(1)-V1.x(1)
row[0] = E.global_index()
col[0] = VertexNodes[0]
col[1] = VertexNodes[1]
col[2] = 0+mesh.num_vertices()
dataX[2:5] = (1./6)*Xtang
dataY[2:5] = (1./6)*Ytang
row[1] = E.global_index()
col[3] = VertexNodes[0]
col[4] = VertexNodes[1]
col[5] = 0+mesh.num_vertices()
dataX[2:5] = (1./6)*Xtang
dataY[2:5] = (1./6)*Ytang
# for i in range(2):
# print i
# for j in range(3):
0+mesh.num_vertices()
print mesh.num_edges()
In [80]:
row = np.zeros(2*mesh.num_edges())
col = np.zeros(12*mesh.num_edges())
dataX = np.zeros(12*mesh.num_edges())
dataY = np.zeros(12*mesh.num_edges())
EdgeDOFmap = Magnetic.dofmap()
k = 0
i = 0
for l in range(mesh.num_cells()):
C = Cell(mesh,l)
CellDOF = EdgeDOFmap.cell_dofs(l)
print CellDOF
# print len(C.entities(1))
jj = 0
for j in range(len(C.entities(1))):
# print j
E = Edge(mesh,C.entities(1)[j])
VertexNodes = E.entities(0)
V1 = Vertex(mesh,VertexNodes[0])
V2 = Vertex(mesh,VertexNodes[1])
Xtang = V2.x(0)-V1.x(0)
Ytang = V2.x(1)-V1.x(1)
row[k] = CellDOF[jj]
jj = jj + 1
kk = i
col[i] = VertexNodes[0]
i = i+1
col[i] = VertexNodes[1]
i = i+1
col[i] = j+mesh.num_vertices()
dataX[kk:i+1] = (1./6)*Xtang
dataY[kk:i+1] = (1./6)*Ytang
k = k+1
row[k] = CellDOF[jj]
jj = jj + 1
kk = i
col[i] = VertexNodes[0]
i = i+1
col[i] = VertexNodes[1]
i = i+1
col[i] = j+mesh.num_vertices()
dataX[kk:i+1] = (1./6)*Xtang
dataY[kk:i+1] = (1./6)*Ytang
In [75]:
col
Out[75]:
In [58]:
C = Cell(mesh,0)
In [62]:
C.entities(2)
Out[62]:
In [5]:
import petsc4py, sys
petsc4py.init(sys.argv)
from petsc4py import PETSc
from dolfin import *
# import matplotlib.pylab as plt
import PETScIO as IO
import numpy as np
import scipy.sparse as sparse
import CheckPetsc4py as CP
import scipy.sparse.linalg as sparselin
import scipy as sp
import scipy.io
mm = 2
for x in xrange(1,mm):
nn = 2**(x)
print nn
mesh = UnitSquareMesh(nn,nn)
print "Num edges ", mesh.num_edges()
print "Num vertices", mesh.num_vertices()
print "Num faces", mesh.num_faces()
order = 2
parameters['reorder_dofs_serial'] = False
V = FunctionSpace(mesh, "N1curl", order)
Q = FunctionSpace(mesh, "CG", order)
(v) = TrialFunction(V)
(u) = TestFunction(V)
m = inner(u,v)*dx
# parameters['linear_algebra_backend'] = 'PETSc'
M = assemble(m)
M = CP.Assemble(M)
ksp = PETSc.KSP().create()
ksp.setOperators(M)
x = M.getVecLeft()
ksp.setFromOptions()
ksp.setType(ksp.Type.PREONLY)
# ksp.setTolerances(1e-1)
ksp.pc.setType(ksp.pc.Type.LU)
OptDB = PETSc.Options()
ksp.setFromOptions()
C = sparse.csr_matrix((V.dim(),Q.dim()))
tic()
for i in range(0,Q.dim()):
uOut = Function(V)
uu = Function(Q)
x = M.getVecRight()
zero = np.zeros((Q.dim(),1))[:,0]
zero[i] = 1
uu.vector()[:] = zero
L = assemble(inner(u, grad(uu))*dx)
# bcb.apply(L)
rhs = IO.arrayToVec(L.array())
ksp.solve(rhs,x)
# x = project(grad(uu),V)
P = x.array
uOut.vector()[:] = P
low_values_indices = np.abs(P) < 1e-3
P[low_values_indices] = 0
low_values_indices = np.abs(P) > 1e-3
P[low_values_indices] = 1/3
#P=np.around(P)
pn = P.nonzero()[0]
for j in range(0,len(pn)):
C[pn[j],i] = P[pn[j]]
del uu
print toc()
Out[5]:
In [6]:
C
Out[6]:
In [ ]:
C