notebook.community

Edit and run 





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]



    


















    






DEBUG:FFC:Reusing form from cache.
DEBUG:FFC:Reusing form from cache.
DEBUG:FFC:Reusing form from cache.
DEBUG:FFC:Reusing form from cache.
DEBUG:FFC:Reusing form from cache.










    






1


W:   [ 0.] Velocity:   [ 0.] Pressure:   [ 0.] Magnetic:   [ 0.] Lagrange:   [ 0.] 



   >>>>>>>>>>>>>>>>>>>>>>>>>>
     MHD 2D Exact Solution:
   >>>>>>>>>>>>>>>>>>>>>>>>>>

 ----------------------
   NS Exact Solution:
 ----------------------
  u = ( x*y*exp(x + y) + x*exp(x + y) , -x*y*exp(x + y) - y*exp(x + y) )

  p = ( exp(y)*sin(x) )

 ---------------------------
   Maxwell Exact Solution:
 ---------------------------
  b = ( exp(x + y)*cos(x) , exp(x + y)*sin(x) - exp(x + y)*cos(x) )

  p = ( x*sin(2*pi*x)*sin(2*pi*y) )



=====================================
  Seting up initial guess matricies
=====================================


















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



    


















    






16

















In [92]:



    


row = np.zeros(Magnetic.dim())
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)

        print CellDOF[jj]
        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

        print CellDOF[jj]
        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()
        k = k+1
        dataX[kk:i+1] = (1./6)*Xtang
        dataY[kk:i+1] = (1./6)*Ytang



    


















    






[ 0  1  2  3  4  5 32 33]
0
1
2
3
4
5
[ 6  7  2  3  8  9 34 35]
6
7
2
3
8
9
[10 11 12 13 14 15 36 37]
10
11
12
13
14
15
[16 17 12 13  0  1 38 39]
16
17
12
13
0
1
[18 19 20 21  6  7 40 41]
18
19
20
21
6
7
[22 23 20 21 24 25 42 43]
22
23
20
21
24
25
[26 27 28 29 16 17 44 45]
26
27
28
29
16
17
[30 31 28 29 18 19 46 47]
30
31
28
29
18
19

















In [89]:



    


row



    


















    
Out[89]:








array([  0.,   1.,   2.,   3.,   4.,   5.,   6.,   7.,   2.,   3.,   8.,
         9.,  10.,  11.,  12.,  13.,  14.,  15.,  16.,  17.,  12.,  13.,
         0.,   1.,  18.,  19.,  20.,  21.,   6.,   7.,  22.,  23.,  20.,
        21.,  24.,  25.,  26.,  27.,  28.,  29.,  16.,  17.,  30.,  31.,
        28.,  29.,  18.,  19.])

















In [58]:



    


C = Cell(mesh,0)



    











In [62]:



    


C.entities(2)



    


















    
Out[62]:








array([], dtype=uint32)

















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



    


















    






DEBUG:FFC:Reusing form from cache.
DEBUG:FFC:Reusing form from cache.
DEBUG:FFC:Reusing form from cache.
DEBUG:FFC:Reusing form from cache.
DEBUG:FFC:Reusing form from cache.
DEBUG:FFC:Reusing form from cache.
DEBUG:FFC:Reusing form from cache.
DEBUG:FFC:Reusing form from cache.
DEBUG:FFC:Reusing form from cache.
DEBUG:FFC:Reusing form from cache.
DEBUG:FFC:Reusing form from cache.
DEBUG:FFC:Reusing form from cache.
DEBUG:FFC:Reusing form from cache.
DEBUG:FFC:Reusing form from cache.
DEBUG:FFC:Reusing form from cache.
DEBUG:FFC:Reusing form from cache.
DEBUG:FFC:Reusing form from cache.
DEBUG:FFC:Reusing form from cache.
DEBUG:FFC:Reusing form from cache.
DEBUG:FFC:Reusing form from cache.
DEBUG:FFC:Reusing form from cache.
DEBUG:FFC:Reusing form from cache.
DEBUG:FFC:Reusing form from cache.
DEBUG:FFC:Reusing form from cache.
DEBUG:FFC:Reusing form from cache.
DEBUG:FFC:Reusing form from cache.
DEBUG:FFC:Reusing form from cache.
DEBUG:FFC:Reusing form from cache.










    






2
Num edges    0
Num vertices 9
Num faces 8
0.0779409408569










    
Out[5]:








<dolfin.cpp.io.VTKPlotter; proxy of <Swig Object of type 'std::shared_ptr< dolfin::VTKPlotter > *' at 0x11a436660> >

















In [6]:



    


C



    


















    
Out[6]:








<48x25 sparse matrix of type '<type 'numpy.float64'>'
	with 160 stored elements in Compressed Sparse Row format>

















In [ ]:



    


C

Content source: wathen/PhD

Similar notebooks:

notebook.community | gallery | about