-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathMain_glacial_quadrTH_new.m
248 lines (218 loc) · 7.96 KB
/
Main_glacial_quadrTH_new.m
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
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
% Modified Taylor-Hood stable discretization
%
%%%%%%%%%%%%%%%%%% 2D elasticity with advection terms %%%%%%%%%%%%%%%%%%%%%%%%%%
% Works on any given coarse mesh, described as {Node, Edge, Face}
% 1. The coarse mesh is refined 'levels' times
% 2. The mesh obtained in 1. is the pressure mesh. It is then refined
% once more to get the mesh for the displacements. At this level
% parent-child relation is kept to utilize the assembly of the
% stiffness matrices
% 3. The viscoelastic solver is run It includes the assembly of the
% corresponding STIFFNESS and MASS matrices
% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
% Block-factorized preconditioner, two-by-two block structure introduced
% Parent-child relation in the mesh refinement introduced
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
close all
%% ------------------------------- global definitions
global average_inner how_many outer_solv inner_solv inner_solver_type
global restart restart_tol
global avcgit_d avcgit_u
global sigma
global h
global solver_type
global actionILU actionFEM
global L11d U11d L11u U11u % lu(A11)/luinc(A11)
global kkk lll
global eps_inner
global theta
% Debug variable
% 1 - print debug texts
% 0 - no debug
global debug;
debug = 0;
%% Variables
% 0 - No additional output(figure and text)
% 1 - Only figures
% 2 - everything
verbose = 1;
np = 4; % number of points in the finite element
nip = 4; % number of integration points
dim = 2; % number of degree of freedom (elasticity part)
ndof = 8; % ndof = nip*dim
% actionILU =[999,1e-1,1e-2,1e-3];
% actionILUt=['exact solve(A11)','cholinc(0.0100) ','cholinc(0.0010) ','cholinc(0.0001) '];
% actionFEM =[0,1,2,3,4,5,1e6,999];
% actionFEMt=['iter=0 ','iter=1 ','iter=2 ','iter=3 ','iter=4 ',...
% 'iter=5 ','iter=50','Z12 '];
% lll = 2;
% kkk = 1;
%% Preparing parameters and mesh
tic
%test_problem = 0; % pure elasticity in unit square, benchmark
test_problem = 1; % visco-elasticity
if test_problem==0,
wh = 'g0';
no_domains = 1;
Emagn = 1;
[L,H,l_ice,h_ice,rho_ice,rho_earth,...
Disco,Discoef,grav,load_surf,...
L_char, S_char, U_char, N_char, T_char, ...
T_LGM, T_EOG, T, delta_t_char] = Elasto_parameters(no_domains,wh,Emagn);
H0=-max(abs(H));
Nx = 3;
Ny = 3;
[xc,yc,hx,hy,Nx,Ny] = Glace_coord_vectors_TH(L,H0,Nx,Ny);
else
wh = 'gs';
no_domains = 2;
Emagn = 1; % can be 1, 10, 100 (jump in E between the two subdomains)
% -- delta_t_char corresponds to one scaled year --
[L,H,l_ice,h_ice,rho_ice,rho_earth,...
Disco,Discoef,grav,load_surf,...
L_char, S_char, U_char, N_char, T_char, ...
T_LGM, T_EOG, T, delta_t_char] = Visco_parameters_new4(no_domains,wh,Emagn);
% T_LGM, T_EOG, T] = Visco_parameters(no_domains,wh,Emagn);
H0=-max(abs(H));
Nx = L/(2*l_ice)+1; % ensure mesh aligned with the ice after one refinement
Ny = abs(H0)/(2*l_ice)+1;
% [xc,yc,hx,hy,Nx,Ny] = Glace_coord_vectors_TH(L,H0,Nx,Ny);
[xc,yc,hx,hy,Nx,Ny] = Glace_coord_vectors_Wu_coarse(L,H,l_ice);
end
[Node,Edge,Face,...
Node_flagx,Node_flagy,...
Edge_flagx,Edge_flagy,...
Face_flag,Face_thick] = Rectan_glace_vect(L,H,xc,yc,Nx,Ny,...
no_domains,Disco,wh);
% Visualise the mesh
if(verbose ~= 0)
figure(1),clf,Bvisual_mesh(Node,Edge,Face,1,1,1,0,16)
end
disp(['Time to create initial mesh: ' num2str(toc)])
% -------------------- Input parameters ---------
levels = input('How many times to refine: ');
h = min(abs(hx),abs(hy))/(2^levels);
sigma = h^2;
disp('sigma = h^2')
solver_type = 1;
tic
for lvl=1:levels,
[Node,Edge,Face,...
Node_flagx,Node_flagy,...
Edge_flagx,Edge_flagy,...
Face_flag,Face_thick] = my_Refine_quadr(Node,Edge,Face,...
Node_flagx,Node_flagy,...
Edge_flagx,Edge_flagy,...
Face_flag,Face_thick);
% figure(1),Bvisual_mesh(Node,Edge,Face,1,1,1,0,16)
end
nnodeP = size(Node,2);
nedgeP = size(Edge,2);
nfaceP = size(Face,2); % number of subdomains (for the pressure)
nallP = nnodeP; % number of pressure variables
disp(['Total number of quadrilaterals for the pressure: ' int2str(nfaceP)])
%figure(1),clf,Bvisual_mesh(Node,Edge,Face,1,1,1,3,16)
disp(['Total number of coarse quadrilaterals: ' int2str(nfaceP)])
%
% ---------------- Create mesh-related matrices -----
%Edge_Node = spalloc(nedgeP,nnodeP,2*nedgeP);
%Face_Edge = spalloc(nfaceP,nedgeP,np*nfaceP);
%for iedge=1:nedgeP,
% Edge_Node(iedge,Edge(1,iedge))=1;
% Edge_Node(iedge,Edge(2,iedge))=1;
%end
%for iface=1:nfaceP,
% Face_Edge(iface,Face(1,iface))=1;
% Face_Edge(iface,Face(2,iface))=1;
% Face_Edge(iface,Face(3,iface))=1;
% Face_Edge(iface,Face(4,iface))=1;
%end
%
%Face_Node_S = Face_Edge*Edge_Node; % Face-to-Node-to-Face relation
%for iface=1:nfaceP,
% rr = find(Face_Node_S(iface,:));
% Face_Node(1:4,iface,1)=rr';
%nd
%clear Face_Node_S
%
%Node_Face_no = sum(Face_Node)./2; % to how many faces each node belongs
% Refine once to obtain the mesh for the displacements
lvl_total = levels + 1;
lvl_coars = max(1,levels-1);
nface_lvl = zeros(lvl_total,1);
nface_lvl(1) = size(Face,2);
nnode_lvl(1) = size(Node,2);
for lvl=1:1, % only once
[Node,Edge,Face,...
Node_flagx,Node_flagy,...
Edge_flagx,Edge_flagy,...
Face_flag,Face_thick,...
Face_Node,Face_Parent,...
Face_eorder11,Face_eorder22,...
nface_lvl,nnode_lvl] = my_Refine_quadr_hier(Node,Edge,Face,...
Node_flagx,Node_flagy,Edge_flagx,Edge_flagy,...
Face_flag,Face_thick,nface_lvl,nnode_lvl,lvl);
% Bvisual_mesh(Node,Edge,Face,1,1,1,1,16)
end
nnode = size(Node,2);
nedge = size(Edge,2);
nface = size(Face,2); % number of subdomains (for the displacements)
% ---------------- Create mesh-related matrices (again) -----
Edge_Node = spalloc(nedge,nnode,2*nedge);
Face_Edge = spalloc(nface,nedge,np*nface);
for iedge=1:nedge,
Edge_Node(iedge,Edge(1,iedge))=1;
Edge_Node(iedge,Edge(2,iedge))=1;
end
for iface=1:nface,
Face_Edge(iface,Face(1,iface))=1;
Face_Edge(iface,Face(2,iface))=1;
Face_Edge(iface,Face(3,iface))=1;
Face_Edge(iface,Face(4,iface))=1;
end
% ------------ detect boundary edges
bndry_edge = zeros(nedge,1);
bndry_edge = sum(Face_Edge,1); % The boundary edges are with sum '1'
[noi,Bndry_Edges]=find(bndry_edge==1); % Bndry_Edges is a list of boundary edges ONLY!
clear bndry_edge
% ------------ detect boundary edges under the ice
% if wh=='gs' % ice for y=0, 0<=x<=length_ice
Surface_Nodes = find(Node(2,:)==0); % all surface nodes
[noi,noj] = find(Node(1,Surface_Nodes)<=l_ice);
Ice_Nodes = Surface_Nodes(noj);
% end
%[noi,noj]=find(Edge_Node(:,Surface_Nodes));
[noi,noj]=find(Edge_Node(:,Ice_Nodes));
noi=unique(noi);
clear Bndry_Ice, lb=0;
for i=1:length(noi),
if (Node(2,Edge(:,noi(i)))==[0 0])&(prod(Node(1,Edge(:,noi(i))))<=l_ice^2),
lb = lb + 1; Bndry_Ice(lb)=noi(i);
end
end
% Node_Ice = [];
% Bndry_Ice = [];
% if wh=='gs' % ice for y=0, 0<=x<=length_ice
% noj = find(Node(1,Surface_Nodes)<=l_ice);
% Node_Ice = Surface_Nodes(noj); % Node number of the nodes under the ice
% Bndry_Ice = 0*Node_Ice;
% for i=1:length(Node_Ice)
% temp = find(Edge(1,:) == Node_Ice(i));
% temp2 = find(Node(2,Edge(2,temp)) == 0);
% Bndry_Ice(i) = temp(temp2);
% end
% end
Node_Face_noP = sum(Face_Node(:,:,1))./2; %to how many faces each P node belongs
Node_Face_noD = sum(Face_Node(:,:,2))./2; %to how many faces each D node belongs
%depth = -2e-4;
%[Surface_Nodes, Depth_Nodes] = Level_Nodes(Node,depth);
%Surface_NodesP = Surface_Nodes(1:Nx+(Nx-1)*(2^levels-1));
disp(['End refinement. Elapsed: ' num2str(toc)])
if size(Node,2)<625, figure(1),clf,Bvisual_mesh(Node,Edge,Face,0,0,0,3,12), end
disp(' Ready with the mesh.')
disp(' Start the visco-elastic solver.')
Tmax = T;%input(' Run till time: ');
eps_inner = 5e-1;
Visco_elastic_solverTH_new
return