In [1]:
#!/usr/bin/python
# interpolate scalar gradient onto nedelec space
import petsc4py
import sys
petsc4py.init(sys.argv)
from petsc4py import PETSc
from dolfin import *
import mshr
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 MHDmatrixPrecondSetup as PrecondSetup
import NSprecondSetup
import MHDprec as MHDpreconditioner
import memory_profiler
import gc
import MHDmulti
import MHDmatrixSetup as MHDsetup
import Lshaped
#@profile
m = 2
set_log_active(False)
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'
MU[0]= 1e0
for xx in xrange(1,m):
print xx
level[xx-1] = xx + 1
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)
# domain = mshr.Rectangle(Point(0., 0.), Point(1., 2.)) + mshr.Rectangle(Point(1., 0.), Point(2., 1.))
# mesh = mshr.generate_mesh(domain, nn)
mesh, boundaries, domains = Lshaped.Domain(nn)
# set_log_level(WARNING)
order = 2
parameters['reorder_dofs_serial'] = False
Velocity = VectorFunctionSpace(mesh, "CG", order)
Pressure = FunctionSpace(mesh, "CG", order-1)
Magnetic = FunctionSpace(mesh, "N1curl", order-1)
Lagrange = FunctionSpace(mesh, "CG", order-1)
W = MixedFunctionSpace([Velocity, Pressure, Magnetic,Lagrange])
# W = Velocity*Pressure*Magnetic*Lagrange
Velocitydim[xx-1] = Velocity.dim()
Pressuredim[xx-1] = Pressure.dim()
Magneticdim[xx-1] = Magnetic.dim()
Lagrangedim[xx-1] = Lagrange.dim()
Wdim[xx-1] = W.dim()
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
FSpaces = [Velocity,Pressure,Magnetic,Lagrange]
kappa = 1.0
Mu_m =1e4
MU = 1.0
N = FacetNormal(mesh)
ds = Measure('ds', domain=mesh, subdomain_data=boundaries)
dx = Measure('dx', domain=mesh, subdomain_data=domains)
# g = inner(p0*N - MU*grad(u0)*N,v)*dx
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]
u0, p0, b0, r0, F_NS, F_M, gradu0 = Lshaped.Solution(mesh, params)
u_k, p_k = Lshaped.Stokes(Velocity, Pressure, F_NS, u0, p0, gradu0, params, boundaries, domains)
b_k, r_k = Lshaped.Maxwell(Magnetic, Lagrange, F_M, b0, r0, params,boundaries, domains)
# print u_k.vector().array()
# print p_k.vector().array()
# print b_k.vector().array()
# print r_k.vector().array()
# bcu = DirichletBC(W.sub(0),u0, boundaries, 1)
# bcb = DirichletBC(W.sub(2),b0, boundaries, 1)
# bcr = DirichletBC(W.sub(3),r0, boundaries, 1)
# print bcu.get_boundary_values()
# print bcb.get_boundary_values()
# print bcr.get_boundary_values()
# ssss
(u, p, b, r) = TrialFunctions(W)
(v, q, c, s) = TestFunctions(W)
n = FacetNormal(mesh)
# u_k = Function(Velocity)
# p_k = Function(Pressure)
# b_k = Function(Magnetic)
# r_k = Function(Lagrange)
m11 = params[1]*params[0]*inner(curl(b),curl(c))*dx(0)
m21 = inner(c,grad(r))*dx(0)
m12 = inner(b,grad(s))*dx(0)
a11 = params[2]*inner(grad(v), grad(u))*dx(0) + inner((grad(u)*u_k),v)*dx(0) #+(1./2)*div(u_k)*inner(u,v)*dx(0) - (1./2)*inner(u_k,n)*inner(u,v)*ds(0)
a12 = -div(v)*p*dx(0)
a21 = -div(u)*q*dx(0)
CoupleT = params[0]*(v[0]*b_k[1]-v[1]*b_k[0])*curl(b)*dx(0)
Couple = -params[0]*(u[0]*b_k[1]-u[1]*b_k[0])*curl(c)*dx(0)
a = m11 + m12 + m21 + a11 + a21 + a12 + Couple + CoupleT
Lns = inner(v, F_NS)*dx(0) + inner(p0*n - MU*gradu0*n,v)*ds(2)
Lmaxwell = inner(c, F_M)*dx(0)
m11 = params[1]*params[0]*inner(curl(b_k),curl(c))*dx(0)
m21 = inner(c,grad(r_k))*dx(0)
m12 = inner(b_k,grad(s))*dx(0)
a11 = params[2]*inner(grad(v), grad(u_k))*dx(0) + inner((grad(u_k)*u_k),v)*dx(0) #+(1./2)*div(u_k)*inner(u_k,v)*dx(0) - (1./2)*inner(u_k,n)*inner(u_k,v)*ds(0)
a12 = -div(v)*p_k*dx(0)
a21 = -div(u_k)*q*dx(0)
CoupleT = params[0]*(v[0]*b_k[1]-v[1]*b_k[0])*curl(b_k)*dx(0)
Couple = -params[0]*(u_k[0]*b_k[1]-u_k[1]*b_k[0])*curl(c)*dx(0)
L = Lns + Lmaxwell - (m11 + m12 + m21 + a11 + a21 + a12 + Couple + CoupleT)
MO.PrintStr("Seting up initial guess matricies",2,"=","\n\n","\n")
BCtime = time.time()
BC = MHDsetup.BoundaryIndices(mesh)
MO.StrTimePrint("BC index function, time: ", time.time()-BCtime)
Hiptmairtol = 1e-6
HiptmairMatrices = PrecondSetup.MagneticSetup(Magnetic, Lagrange, b0, r0, Hiptmairtol, params)
MO.PrintStr("Setting up MHD initial guess",5,"+","\n\n","\n\n")
# u_k,p_k,b_k,r_k = common.InitialGuess(FSpaces,[u0,p0,b0,r0],[F_NS,F_M],params,HiptmairMatrices,1e-10,Neumann=None,options ="New")
ones = Function(Pressure)
ones.vector()[:]=(0*ones.vector().array()+1)
# pConst = - assemble(p_k*dx)/assemble(ones*dx)
p_k.vector()[:] += - assemble(p_k*dx)/assemble(ones*dx)
x = Iter.u_prev(u_k,p_k,b_k,r_k)
KSPlinearfluids, MatrixLinearFluids = PrecondSetup.FluidLinearSetup(Pressure, MU)
kspFp, Fp = PrecondSetup.FluidNonLinearSetup(Pressure, MU, u_k)
#plot(b_k)
# ns,maxwell,CoupleTerm,Lmaxwell,Lns = forms.MHD2D(mesh, W,F_M,F_NS, u_k,b_k,params,IterType,"CG",Saddle,Stokes)
# RHSform = forms.PicardRHS(mesh, W, u_k, p_k, b_k, r_k, params,"CG",Saddle,Stokes)
# bcu = DirichletBC(W.sub(0),Expression(("0.0","0.0")), boundaries, 1)
# bcb = DirichletBC(W.sub(2),Expression(("0.0","0.0")), boundaries, 1)
# bcr = DirichletBC(W.sub(3),Expression(("0.0")), boundaries, 1)
# bcs = [bcu,bcb,bcr]
IS = MO.IndexSet(W, 'Blocks')
parameters['linear_algebra_backend'] = 'uBLAS'
eps = 1.0 # error measure ||u-u_k||
tol = 1.0E-4 # tolerance
iter = 0 # iteration counter
maxiter = 10 # max no of iterations allowed
SolutionTime = 0
outer = 0
# parameters['linear_algebra_backend'] = 'uBLAS'
# FSpaces = [Velocity,Magnetic,Pressure,Lagrange]
u_is = PETSc.IS().createGeneral(range(Velocity.dim()))
NS_is = PETSc.IS().createGeneral(range(Velocity.dim()+Pressure.dim()))
M_is = PETSc.IS().createGeneral(range(Velocity.dim()+Pressure.dim(),W.dim()))
OuterTol = 1e-5
InnerTol = 1e-5
NSits =0
Mits =0
TotalStart =time.time()
SolutionTime = 0
# while eps > tol and iter < maxiter:
# iter += 1
# MO.PrintStr("Iter "+str(iter),7,"=","\n\n","\n\n")
# # if iter == 1:
# # bcu = DirichletBC(W.sub(0),u0, boundaries, 1)
# # bcb = DirichletBC(W.sub(2),b0, boundaries, 1)
# # bcr = DirichletBC(W.sub(3),r0, boundaries, 1)
# # bcs = [bcu,bcb,bcr]
# # else:
# bcu = DirichletBC(W.sub(0),Expression(("0.0","0.0")), boundaries, 1)
# bcb = DirichletBC(W.sub(2),Expression(("0.0","0.0")), boundaries, 1)
# bcr = DirichletBC(W.sub(3),Expression("0.0"), boundaries, 1)
# bcs = [bcu,bcb,bcr]
# # if iter == 1:
# # , L
# A, b = assemble_system(a, L, bcs)
# # AA = assemble(a)
# # bb = assemble(L)
# # for bc in bcs:
# # bc.apply(AA,bb)
# # print A.sparray().todense()
# # MO.StoreMatrix(A.sparray(),'name')
# A, b = CP.Assemble(A,b)
# u = b.duplicate()
# # print b.array
# # ssss
# # L = assemble(L)
# # print L.array()
# # for bc in bcs:
# # bc.apply(L)
# # print L.array()
# # MO.StrTimePrint("MHD total assemble, time: ", time.time()-AssembleTime)
# # u = b.duplicate()
# # kspFp, Fp = PrecondSetup.FluidNonLinearSetup(Pressure, MU, u_k)
# # print "Inititial guess norm: ", u.norm(PETSc.NormType.NORM_INFINITY)
# # #A,Q
# stime = time.time()
# u, mits,nsits = S.solve(A,b,u,params,W,'Direct',IterType,OuterTol,InnerTol,HiptmairMatrices,Hiptmairtol,KSPlinearfluids, Fp,1)
# Soltime = time.time()- stime
# MO.StrTimePrint("MHD solve, time: ", Soltime)
# Mits += mits
# NSits += nsits
# SolutionTime += Soltime
# u1, p1, b1, r1, eps= Iter.PicardToleranceDecouple(u,x,FSpaces,dim,"2",iter)
# p1.vector()[:] += - assemble(p1*dx)/assemble(ones*dx)
# u_k.assign(u1)
# p_k.assign(p1)
# b_k.assign(b1)
# r_k.assign(r1)
# uOld= np.concatenate((u_k.vector().array(),p_k.vector().array(),b_k.vector().array(),r_k.vector().array()), axis=0)
# x = IO.arrayToVec(uOld)
XX= np.concatenate((u_k.vector().array(),p_k.vector().array(),b_k.vector().array(),r_k.vector().array()), axis=0)
SolTime[xx-1] = SolutionTime/iter
NSave[xx-1] = (float(NSits)/iter)
Mave[xx-1] = (float(Mits)/iter)
iterations[xx-1] = iter
TotalTime[xx-1] = time.time() - TotalStart
dim = [Velocity.dim(), Pressure.dim(), Magnetic.dim(),Lagrange.dim()]
ExactSolution = [u0,p0,b0,r0]
errL2u[xx-1], errH1u[xx-1], errL2p[xx-1], errL2b[xx-1], errCurlb[xx-1], errL2r[xx-1], errH1r[xx-1] = Iter.Errors(x,mesh,FSpaces,ExactSolution,order,dim, "CG")
if xx > 1:
l2uorder[xx-1] = np.abs(np.log2(errL2u[xx-2]/errL2u[xx-1]))
H1uorder[xx-1] = np.abs(np.log2(errH1u[xx-2]/errH1u[xx-1]))
l2porder[xx-1] = np.abs(np.log2(errL2p[xx-2]/errL2p[xx-1]))
l2border[xx-1] = np.abs(np.log2(errL2b[xx-2]/errL2b[xx-1]))
Curlborder[xx-1] = np.abs(np.log2(errCurlb[xx-2]/errCurlb[xx-1]))
l2rorder[xx-1] = np.abs(np.log2(errL2r[xx-2]/errL2r[xx-1]))
H1rorder[xx-1] = np.abs(np.log2(errH1r[xx-2]/errH1r[xx-1]))
import pandas as pd
LatexTitles = ["l","DoFu","Dofp","V-L2","L2-order","V-H1","H1-order","P-L2","PL2-order"]
LatexValues = np.concatenate((level,Velocitydim,Pressuredim,errL2u,l2uorder,errH1u,H1uorder,errL2p,l2porder), axis=1)
LatexTable = pd.DataFrame(LatexValues, columns = LatexTitles)
pd.set_option('precision',3)
LatexTable = MO.PandasFormat(LatexTable,"V-L2","%2.4e")
LatexTable = MO.PandasFormat(LatexTable,'V-H1',"%2.4e")
LatexTable = MO.PandasFormat(LatexTable,"H1-order","%1.2f")
LatexTable = MO.PandasFormat(LatexTable,'L2-order',"%1.2f")
LatexTable = MO.PandasFormat(LatexTable,"P-L2","%2.4e")
LatexTable = MO.PandasFormat(LatexTable,'PL2-order',"%1.2f")
print LatexTable.to_latex()
print "\n\n Magnetic convergence"
MagneticTitles = ["l","B DoF","R DoF","B-L2","L2-order","B-Curl","HCurl-order"]
MagneticValues = np.concatenate((level,Magneticdim,Lagrangedim,errL2b,l2border,errCurlb,Curlborder),axis=1)
MagneticTable= pd.DataFrame(MagneticValues, columns = MagneticTitles)
pd.set_option('precision',3)
MagneticTable = MO.PandasFormat(MagneticTable,"B-Curl","%2.4e")
MagneticTable = MO.PandasFormat(MagneticTable,'B-L2',"%2.4e")
MagneticTable = MO.PandasFormat(MagneticTable,"L2-order","%1.2f")
MagneticTable = MO.PandasFormat(MagneticTable,'HCurl-order',"%1.2f")
print MagneticTable.to_latex()
print "\n\n Lagrange convergence"
LagrangeTitles = ["l","B DoF","R DoF","R-L2","L2-order","R-H1","H1-order"]
LagrangeValues = np.concatenate((level,Magneticdim,Lagrangedim,errL2r,l2rorder,errH1r,H1rorder),axis=1)
LagrangeTable= pd.DataFrame(LagrangeValues, columns = LagrangeTitles)
pd.set_option('precision',3)
LagrangeTable = MO.PandasFormat(LagrangeTable,"R-L2","%2.4e")
LagrangeTable = MO.PandasFormat(LagrangeTable,'R-H1',"%2.4e")
LagrangeTable = MO.PandasFormat(LagrangeTable,"L2-order","%1.2f")
LagrangeTable = MO.PandasFormat(LagrangeTable,'H1-order',"%1.2f")
print LagrangeTable.to_latex()
import pandas as pd
print "\n\n Iteration table"
if IterType == "Full":
IterTitles = ["l","DoF","AV solve Time","Total picard time","picard iterations","Av Outer its","Av Inner its",]
else:
IterTitles = ["l","DoF","AV solve Time","Total picard time","picard iterations","Av NS iters","Av M iters"]
IterValues = np.concatenate((level,Wdim,SolTime,TotalTime,iterations,Mave,NSave),axis=1)
IterTable= pd.DataFrame(IterValues, columns = IterTitles)
if IterType == "Full":
IterTable = MO.PandasFormat(IterTable,'Av Outer its',"%2.1f")
IterTable = MO.PandasFormat(IterTable,'Av Inner its',"%2.1f")
else:
IterTable = MO.PandasFormat(IterTable,'Av NS iters',"%2.1f")
IterTable = MO.PandasFormat(IterTable,'Av M iters',"%2.1f")
print IterTable
print " \n Outer Tol: ",OuterTol, "Inner Tol: ", InnerTol
# tableName = "2d_Lshaped_nu="+str(MU)+"_nu_m="+str(Mu_m)+"_kappa="+str(kappa)+"_l="+str(np.min(level))+"-"+str(np.max(level))+"Approx.tex"
# IterTable.to_latex(tableName)
# # # if (ShowResultPlots == 'yes'):
# plot(u_k)
# plot(interpolate(u0,Velocity))
#
# plot(p_k)
#
# plot(interpolate(p0,Pressure))
#
# plot(b_k)
# plot(interpolate(b0,Magnetic))
#
# plot(r_k)
# plot(interpolate(r0,Lagrange))
#
# interactive()
interactive()
In [3]:
import mshr
mshr.generate_mesh?
In [9]:
RectangleMesh(-1,0,1,1,4,2)
domain = mshr.Rectangle(Point(-1., -1.), Point(1., 1.)) - mshr.Rectangle(Point(0., -1.), Point(1., 0.) ) mesh = mshr.generate_mesh(domain, 4)
mesh.smooth()
In [19]:
mesh = RectangleMesh(-1,-1,1,1,4,4)
In [20]:
mesh
Out[20]:
In [57]:
print np.arange(1,6,2)
cell_f = CellFunction('size_t', mesh, 0)
for cell in cells(mesh):
v = cell.get_vertex_coordinates()
y = v[np.arange(0,6,2)]
x = v[np.arange(1,6,2)]
xone = np.ones(3)
xone[x <0] = 0
yone = np.ones(3)
yone[x <0] = 0
print x,xone
# for i in range(len(v)):
# print v[0:2]
v = cell.get_vertex_coordinates()