-
Notifications
You must be signed in to change notification settings - Fork 13
/
transientTerms.jl
109 lines (82 loc) · 3.09 KB
/
transientTerms.jl
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
# ===============================
# Written by AAE
# TU Delft, Spring 2014
# simulkade.com
# ===============================
function transientTerm(phi_old::CellValue, dt::Real, alfa::Real)
transientTerm(phi_old, dt, alfa.*ones(Float64,tuple(phi_old.domain.dims...)))
end
function transientTerm(phi_old::CellValue, dt::Real)
transientTerm(phi_old, dt, 1.0)
end
function transientTerm(phi_old::CellValue, dt::Real, alfa::CellValue)
d = phi_old.domain.dimension
if d==1 || d==1.5
transientTerm1D(phi_old, dt, alfa.value[2:end-1])
elseif d==2 || d==2.5 || d==2.8
transientTerm2D(phi_old, dt, alfa.value[2:end-1,2:end-1])
elseif d==3 || d==3.2
transientTerm3D(phi_old, dt, alfa[2:end-1,2:end-1,2:end-1])
end
end
function transientTerm(phi_old::CellValue, dt::Real, alfa::Array{T}) where T<:Real
d = phi_old.domain.dimension
if d==1 || d==1.5
transientTerm1D(phi_old, dt, alfa)
elseif d==2 || d==2.5 || d==2.8
transientTerm2D(phi_old, dt, alfa)
elseif d==3 || d==3.2
transientTerm3D(phi_old, dt, alfa)
end
end
function transientTerm1D(phi_old::CellValue,
dt::Real, alfa::Array{T}) where T<:Real
# returns the matrix and RHS for a d(phi)/dt term
# extract data from the mesh structure
Nx = phi_old.domain.dims[1]
G = [1:Nx+2;]
# rearrange the matrix of k and build the sparse matrix for internal cells
row_index = reshape(G[2:Nx+1],Nx) # main diagonal (only internal cells)
AP_diag = reshape(alfa/dt,Nx)
M = sparse(row_index, row_index, AP_diag, Nx+2, Nx+2)
# define the RHS Vector
RHS = zeros(Nx+2)
# assign the values of the RHS vector
RHS[row_index] .= reshape(alfa.*phi_old.value[2:Nx+1]/dt,Nx)
(M, RHS)
end
function transientTerm2D(phi_old::CellValue,
dt::Real, alfa::Array{T}) where T<:Real
# returns the matrix and RHS for a d(phi)/dt term
# extract data from the mesh structure
Nx = phi_old.domain.dims[1]
Ny = phi_old.domain.dims[2]
G=reshape([1:(Nx+2)*(Ny+2);], Nx+2, Ny+2)
# rearrange the matrix of k and build the sparse matrix for internal cells
row_index = reshape(G[2:Nx+1,2:Ny+1],Nx*Ny) # main diagonal (only internal cells)
AP_diag = reshape(alfa/dt,Nx*Ny)
M = sparse(row_index, row_index, AP_diag, (Nx+2)*(Ny+2), (Nx+2)*(Ny+2))
# define the RHS Vector
RHS = zeros((Nx+2)*(Ny+2))
# assign the values of the RHS vector
RHS[row_index] .= reshape(alfa.*phi_old.value[2:Nx+1,2:Ny+1]/dt,Nx*Ny)
(M, RHS)
end
function transientTerm3D(phi_old::CellValue,
dt::Real, alfa::Array{T}) where T<:Real
# returns the matrix and RHS for a d(phi)/dt term
# extract data from the mesh structure
Nx = phi_old.domain.dims[1]
Ny = phi_old.domain.dims[2]
Nz = phi_old.domain.dims[3]
G=reshape([1:(Nx+2)*(Ny+2)*(Nz+2);], Nx+2, Ny+2, Nz+2)
# rearrange the matrix of k and build the sparse matrix for internal cells
row_index = reshape(G[2:Nx+1,2:Ny+1,2:Nz+1],Nx*Ny*Nz) # main diagonal (only internal cells)
AP_diag = reshape(alfa/dt,Nx*Ny*Nz)
M = sparse(row_index, row_index, AP_diag, (Nx+2)*(Ny+2)*(Nz+2), (Nx+2)*(Ny+2)*(Nz+2))
# define the RHS Vector
RHS = zeros((Nx+2)*(Ny+2)*(Nz+2))
# assign the values of the RHS vector
RHS[row_index] .= reshape(alfa.*phi_old.value[2:Nx+1,2:Ny+1,2:Nz+1]./dt,Nx*Ny*Nz)
(M, RHS)
end