-
Notifications
You must be signed in to change notification settings - Fork 0
/
exp_2_quadrate.jl
143 lines (102 loc) · 4.85 KB
/
exp_2_quadrate.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
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
# armijo_bas = 0.5
# armijo_sig = 0
# armijo_maxtry = 40
armijo_bas = 0.7
armijo_sig = 0
armijo_maxtry = 40
grad_bound = 1e-8
@everywhere const alpha = 0.001
@everywhere const beta = 0.001
# maxsteps = 1
maxsteps = 100000
save_every = 0
@everywhere time_regularization = false # geht nicht mit velocities_at interfaces
# velocities_at = "interfaces"
velocities_at = "centers"
transport_parallel = false # geht nicht gut, erst ab ca 500x500 Pixel sinnvoll
@everywhere interpolate_w_time = false
# das Verfahren mit Zeitregularisierung parallelisiert
# automatisch die Dimensionen, wenn mehr als ein Worker existiert
grad_parallel = true # betrifft nur die Verfahren ohne Zeitregularisierung
project_divfree = false # betrifft nur velocities_at = "interfaces"
#thr diese Optionen funktionieren nicht alle und die meisten sind sowieso unsinnvoll.
#poisson_solver = "multig" #geht nicht parallel
#poisson_solver = "gmres" #ungenau
poisson_solver = "lufact" #fur gegebene Probleme am besten. Eigentlich Cholesky-Faktorisierung fuer die interfaces und LU-Faktorisierung fuer center
#stokes_solver = "multig" #schlecht geeignet, langsam
#stokes_solver = "gmres" #ungenau
stokes_solver = "lufact"#fur gegebene Probleme am besten
timereg_solver = "multig"#fur gegebene Probleme am besten
#timereg_solver = "gmres" #ungenau
#timereg_solver = "lufact" #nur fuer sehr kleine Probleme benutzbar
#multigrid solver tolerance
@everywhere const mg_tol = 1e-1
@everywhere with_cfl_check = true
# Zeitregularisierung funktioniert nur mit Flussdiskretisierung an Zellmittelpunkten
# diese Zeile ist zu Sicherheit, damit man nichts falsch einstellt
# velocities_at = ~time_regularization ? velocities_at : "centers"
include("view.jl")
@everywhere const m = 60
@everywhere const n = 60
# @everywhere const m = 120
# @everywhere const n = 120
include("beispiele.jl")
# fuer die Konstruktion der Zeitregularisierungsmatrizen muss n_samples >=2 und n_zwischensamples >=3 sein!
@everywhere const n_samples = 5
@everywhere const auslassen = 7 # die Referenzsamples werden so gewählt, dass aus der Vorgabe werden immer `auslassen` Frames weggelassen werden
@everywhere const zwischen_ausgelassen = 3 # zwischen zwei ausgelassenen Frames sollen so viele Zwischenframes generiert werden.
# @everywhere const zwischen_ausgelassen = 12 # zwischen zwei ausgelassenen Frames sollen so viele Zwischenframes generiert werden.
# die Anzahl zwischen den Referenzframes zu generierenden Frames.
@everywhere const n_zwischensamples = auslassen + (auslassen+1) * zwischen_ausgelassen
# ...................... T, alle Frames/Zeitpunkte, also T-1 Zeitschritte von einem Frame auf den naechsten
@everywhere const T = (n_samples-1)*(n_zwischensamples+1) +1
@everywhere const T_vorgabe = auswahl_vorgabe(auslassen, n_samples)[end]
@everywhere const dt = 1/(T-1)
@everywhere const dx = 1/(max(m,n) -1)
iv1 = init_vorgabe(char_quadratR, m,n, T_vorgabe)
iv2 = init_vorgabe(char_quadratL, m,n, T_vorgabe)
I_vorgabe = iv1+iv2
# I_vorgabe = init_vorgabe(rot_circle_ex, 2*m,2*n, T_vorgabe)[m+1:2*m, n+1:2*n, :]
# I_vorgabe = init_vorgabe(transl_circle, m, n, T_vorgabe)
# I_vorgabe = init_vorgabe(deform_circle, m, n, T_vorgabe)
# I_vorgabe = init_vorgabe(__rot_circle_ex, m,n, T_vorgabe)
# I_vorgabe = readtaxi_alt()[3:192, 3:end-2, 1:T_vorgabe] # die dlm-dateien wurden am Rand mit Nullen aufgefuellt
# @everywhere rfac=0.3
# srand(1) #random seed setzen, damit fuer verschiedene Durchlaeufe vergleichbarere Ergebnisse
# randerr = randn(size(I_vorgabe))
# I_vorgabe+= rfac*randerr
s = I_vorgabe[:,:,auswahl_vorgabe(auslassen, n_samples)]
# Aufloesungsunabhaengig, ohne Vorgabe
# s[:,:,1] = zeros(m,n)
# s[:,:,2] = zeros(m,n)
# for i = 1:m
# for j = 1:n
# s[i,j,1] = ((j-1)*dx >= 0.3) & ((j-1)*dx<=0.6) & (((i-1)*dx) >=0.3 )& (((i-1)*dx)<=0.6)
# s[i,j,2] = ((j-1)*dx >= 0.5) & ((j-1)*dx<=0.8) & (((i-1)*dx) >=0.5 )& (((i-1)*dx)<=0.8)
# end
# end
velocities_at == "centers" && begin
u = 0* ones( m, n, T-1 )
v = 0* ones( m, n, T-1 )
end
velocities_at == "interfaces" && begin
u = 0* ones( m, n-1, T-1 )
v = 0* ones( m-1, n, T-1 )
end
include("verfahren.jl")
echo=_echolog
@everywhere rootdir = "../out/new/2quadrate_eng/time_reg_$(time_regularization)/$(m)_x_$(n)_$(n_samples)_$(n_zwischensamples)_$(alpha)_$(beta)/"
make_output_dir(rootdir)
@time I, u, v, p, L2_errs, H1_errs, Js, H1_J_ws, steps = verfahren_grad(s, u, v, 1, 1.0, grad_bound)
save_endergebnis(rootdir)
# # Differenz zur Vorgabe
vorgabe_fehler = diff_vorgabe(I_vorgabe, I, auslassen, zwischen_ausgelassen)
echo("L2( I-I_vorgabe )", L2norm(vorgabe_fehler))
echo("linf( I-I_vorgabe )", l_inf(vorgabe_fehler))
echo("PNSR( I-I_vorgabe )", psnr(vorgabe_fehler))
echo("Gradnorm", H1_J_ws[end])
demo_table("test", "test")
# save_demo_quadrat_schwer([(".png", 100), (".eps", 1200)])
# save_all()
_="fertig"
pygui(true)