-
Notifications
You must be signed in to change notification settings - Fork 1
/
linear_stommel_sw_fenics.py
123 lines (94 loc) · 2.95 KB
/
linear_stommel_sw_fenics.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
from dolfin import *
from mshr import *
import numpy as np
#OPERATORS
zcross = lambda i: as_vector((-i[1],i[0]))
#BASIN & MESH
length = 1.0
width = 1.0#np.sqrt(2)
resolution = 25
geometry = Rectangle(Point(0.0, 0.0), Point(width, length))
mesh = generate_mesh(geometry, resolution)
#FUNCTION & VECTOR SPACES
DGv = FiniteElement("BDM", mesh.ufl_cell(), 2)
CG = FiniteElement("DG", mesh.ufl_cell(), 1)
G = FunctionSpace(mesh, MixedElement((DGv, CG)))
#TRIAL/TEST FUNCTIONS
(u, eta) = TrialFunctions(G)
(v, lmbda) = TestFunctions(G)
#SOLUTION SPACES
sol = Function(G)
#VARIABLES
r = Constant("0.1")
beta = Constant("1.0")
tau = Constant("0.001")
F = Constant("0.1")
Ed = Constant("0.0")
fcor = Expression("1.0+beta*x[1]", beta=beta, degree=2)
class Source(Expression):
def eval(self, values, x):
#values[0] = -tau*sin((pi*x[1])/(2*length))
values[0] = (length*tau/pi)*sin((pi*(x[1]-length/2.0))/length)
#Fwinds = Source(tau = tau, degree = 2)
Fwinds = Expression("(length*tau/pi)*sin((pi*(x[1]-length/2.0))/length)", tau = tau, length = length, degree=2)
# ==========================================================================
#CGs = FunctionSpace(mesh,"CG",3)
#tmp = project(betay, CGs)
#plot(tmp)
#interactive()
#BOUNDARY CONDITIONS
# Define function G such that G \cdot n = g
class BoundarySource(Expression):
def __init__(self, mesh, **kwargs):
self.mesh = mesh
def eval_cell(self, values, x, ufc_cell):
cell = Cell(self.mesh, ufc_cell.index)
n = cell.normal(ufc_cell.local_facet)
g = 0.0
values[0] = g*n[0]
values[1] = g*n[1]
def value_shape(self):
return (2,)
GS = BoundarySource(mesh, degree=2)
# Define essential boundary
def boundary(x):
return x[1] < DOLFIN_EPS or x[1] > 1.0 - DOLFIN_EPS or x[0] < DOLFIN_EPS or x[0] > 1.0 - DOLFIN_EPS
BCs = DirichletBC(G.sub(0), GS, boundary)
# ==========================================================================
F = inner(fcor*v, zcross(u))*dx + \
-div(v)*eta*dx \
+ r*inner(v, u)*dx + lmbda*div(u)*dx - \
Fwinds*v[0]*dx + F*eta*lmbda*dx
# Ed*inner(nabla_grad(lmbda), nabla_grad(eta))*dx + \
#inner(v, grad(eta))*dx + \
a, L = lhs(F), rhs(F)
solve(a==L, sol, BCs)
"""
#u, eta = sol.split()
# Assemble system
A = assemble(a)
b = assemble(L)
# Create Krylov solver
solver = PETScKrylovSolver("cg") #????
solver.set_operator(A)
# Create vector that spans the null space and normalize
null_vec = Vector(sol.vector())
G.dofmap().set(null_vec, 1.0)
null_vec *= 1.0/null_vec.norm("l2")
# Create null space basis object and attach to PETSc matrix
null_space = VectorSpaceBasis([null_vec])
as_backend_type(A).set_nullspace(null_space)
# Orthogonalize RHS vector b with respect to the null space (this
# gurantees a solution exists)
null_space.orthogonalize(b);
# Solve
solver.solve(sol.vector(), b)
"""
# Split
u, eta = sol.split()
File("u.pvd") << u
plot(u)
interactive()
File("eta.pvd") << eta
plot(eta)
interactive()