SPEED
COMPUTE_ENERGY_ERROR.f90
Go to the documentation of this file.
1! Copyright (C) 2012 The SPEED FOUNDATION
2! Author: Ilario Mazzieri
3!
4! This file is part of SPEED.
5!
6! SPEED is free software; you can redistribute it and/or modify it
7! under the terms of the GNU Affero General Public License as
8! published by the Free Software Foundation, either version 3 of the
9! License, or (at your option) any later version.
10!
11! SPEED is distributed in the hope that it will be useful, but
12! WITHOUT ANY WARRANTY; without even the implied warranty of
13! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
14! Affero General Public License for more details.
15!
16! You should have received a copy of the GNU Affero General Public License
17! along with SPEED. If not, see <http://www.gnu.org/licenses/>.
18
78
80
81 subroutine compute_energy_error(nnod_loc, u1, v1, time, ne_loc, cs_loc, cs_nnz_loc,&
82 nm,prop_mat, sdeg_mat,tag_mat, &
83 loc_n_num, local_el_num, nelem_dg_local, &
84 alfa11,alfa12,alfa13,&
85 alfa21,alfa22,alfa23,&
86 alfa31,alfa32,alfa33,&
87 beta11,beta12,beta13,&
88 beta21,beta22,beta23,&
89 beta31,beta32,beta33,&
90 gamma1,gamma2,gamma3,&
91 delta1,delta2,delta3,&
92 xs_loc,ys_loc,zs_loc,&
93 mpi_id,el_new, &
94 nelem_dg_global, &
95 nsd_jump,node_sd_jump, &
96 nrv_jump,node_rv_jump, &
97 send_length_jump, recv_length_jump, &
98 mpi_np, mpi_comm, &
99 cs_nnz_dg, cs_dg)
100
101
102 use max_var
103 use dgjump
104 use speed_exit_codes
105
106 implicit none
107
108 include 'SPEED.MPI'
109
110 type(el4loop), dimension(ne_loc), intent(in) :: el_new
111
112 character*100000 :: input_line
113 character*4 :: keyword
114 character*70 :: file_mpi
115
116 integer*4 :: nnod_loc, nfunc, mpi_id, cs_nnz_loc, nelem_dg_global, nsd_jump, nrv_jump
117 integer*4 :: nm, unit_mpi, ne_loc, status, ic, id, ie_curr, jb,kb,ib,isb
118 integer*4 :: mpi_np, mpi_comm, mpi_ierr, trofa, posit, cs_nnz_dg
119 integer*4 :: ielem, iene, iene_curr, imne, jshift,mm, nelem_dg_local
120 integer*4 :: ie, im, nn, is,in,iaz,i,j,k,ileft,iright,nval
121 integer*4, dimension(0:cs_nnz_loc) :: cs_loc
122 integer*4, dimension(0:cs_nnz_dg) :: cs_dg
123 integer*4, dimension(nnod_loc) :: loc_n_num
124 integer*4, dimension(ne_loc) :: local_el_num
125 integer*4, dimension(nsd_jump):: node_sd_jump
126 integer*4, dimension(nrv_jump):: node_rv_jump
127 integer*4, dimension(nm) :: sdeg_mat
128 integer*4, dimension(nm) :: tag_mat
129 integer*4, dimension(SPEED_STATUS_SIZE) :: mpi_stat
130 integer*4, dimension(mpi_np) :: send_length_jump, recv_length_jump
131
132 real*8 :: time, term1,term2,term3, l2_err, en_err, term_ener,term_l2
133 real*8 :: linf_err, en_err_tot, l2_err_tot, linf_err_tot, xp, yp, zp, phi
134
135 real*8 :: l2_err_vel, h1_err, term_h1, h1_err_tot, term_l2_vel
136 real*8 :: linf_err_vel, l2_err_vel_tot, linf_err_vel_tot
137 real*8 :: e11_ex, e12_ex, e13_ex, e22_ex, e23_ex, e33_ex
138 real*8 :: sxx_ex, syy_ex, szz_ex, sxy_ex, szx_ex, syz_ex
139
140
141 real*8 :: a1,a2,a3,b1,b2,b3,c1,c2,c3,w1, rho,lambda,mu, cost1,cost2,cost3,pi,x,y,z, u1_ctm, u2_ctm, u3_ctm
142 real*8 :: dxdx,dxdy,dxdz,dydx,dydy,dydz,dzdx,dzdy,dzdz,det_j, u1_ex,u2_ex,u3_ex, v1_ex, v2_ex, v3_ex
143
144 real*8, dimension(4) :: linf_array, linf_array_vel
145 real*8, dimension(nm,4) :: prop_mat
146 real*8, dimension(ne_loc) :: alfa11,alfa12,alfa13,alfa21,alfa22,alfa23,alfa31,alfa32,alfa33
147 real*8, dimension(ne_loc) :: beta11,beta12,beta13,beta21,beta22,beta23,beta31,beta32,beta33
148 real*8, dimension(ne_loc) :: gamma1,gamma2,gamma3,delta1,delta2,delta3
149 real*8, dimension(3*nnod_loc) :: u1, v1
150 real*8, dimension(nnod_loc) :: xs_loc, ys_loc, zs_loc
151 real*8, dimension(3*nsd_jump) :: send_buffer_error
152 real*8, dimension(3*nrv_jump) :: recv_buffer_error, jump_minus_error
153
154 real*8, dimension(:), allocatable :: ct, ww, ctm, wwm, u_p, u_m, jump_all
155 real*8, dimension(:,:), allocatable :: dd, ddm
156 real*8, dimension(:), allocatable :: val_data, basis_functions
157 real*8, dimension(:,:,:), allocatable :: ux_el,uy_el,uz_el
158 real*8, dimension(:,:,:), allocatable :: duxdx_el,duydx_el,duzdx_el
159 real*8, dimension(:,:,:), allocatable :: duxdy_el,duydy_el,duzdy_el
160 real*8, dimension(:,:,:), allocatable :: duxdz_el,duydz_el,duzdz_el
161 real*8, dimension(:,:,:), allocatable :: sxx_el,syy_el,szz_el
162 real*8, dimension(:,:,:), allocatable :: syz_el,szx_el,sxy_el
163 real*8, dimension(:,:,:), allocatable :: lambda_el,mu_el
164
165 pi = 4.d0*datan(1.d0)
166
167
168
169 l2_err = 0.d0; en_err = 0.d0; linf_err = 0.d0;
170 l2_err_tot = 0.d0; en_err_tot = 0.d0; linf_err_tot = 0.d0;
171
172 l2_err_vel = 0.d0; h1_err = 0.d0; linf_err_vel = 0.d0;
173 l2_err_vel_tot = 0.d0; h1_err_tot = 0.d0; linf_err_vel_tot = 0.d0;
174
175 do ie = 1,ne_loc
176
177 im = cs_loc(cs_loc(ie -1))
178 nn = sdeg_mat(im) +1;
179 mm = nn ! 8 ! gll nodes for computing norms
180 allocate(ct(nn),ww(nn),dd(nn,nn)); allocate(ctm(mm),wwm(mm),ddm(mm,mm));
181
182 allocate(ux_el(nn,nn,nn),uy_el(nn,nn,nn),uz_el(nn,nn,nn))
183 allocate(duxdx_el(nn,nn,nn),duydx_el(nn,nn,nn),duzdx_el(nn,nn,nn))
184 allocate(duxdy_el(nn,nn,nn),duydy_el(nn,nn,nn),duzdy_el(nn,nn,nn))
185 allocate(duxdz_el(nn,nn,nn),duydz_el(nn,nn,nn),duzdz_el(nn,nn,nn))
186 allocate(sxx_el(nn,nn,nn),syy_el(nn,nn,nn),szz_el(nn,nn,nn))
187 allocate(syz_el(nn,nn,nn),szx_el(nn,nn,nn),sxy_el(nn,nn,nn))
188 allocate(lambda_el(nn,nn,nn),mu_el(nn,nn,nn))
189
190 call make_lgl_nw(nn,ct,ww,dd)
191 call make_lgl_nw(mm,ctm,wwm,ddm)
192
193
194 rho = prop_mat(im,1); lambda = prop_mat(im,2); mu = prop_mat(im,3)
195
196
197
198 do k = 1,nn
199 do j = 1,nn
200 do i = 1,nn
201 is = nn*nn*(k -1) +nn*(j -1) +i
202 in = cs_loc(cs_loc(ie -1) + is)
203
204
205 iaz = 3*(in -1) +1; ux_el(i,j,k) = u1(iaz)
206 iaz = 3*(in -1) +2; uy_el(i,j,k) = u1(iaz)
207 iaz = 3*(in -1) +3; uz_el(i,j,k) = u1(iaz)
208
209 enddo
210 enddo
211 enddo
212
213
214 call make_strain_tensor(nn,ct,ww,dd,&
215 alfa11(ie),alfa12(ie),alfa13(ie),&
216 alfa21(ie),alfa22(ie),alfa23(ie),&
217 alfa31(ie),alfa32(ie),alfa33(ie),&
218 beta11(ie),beta12(ie),beta13(ie),&
219 beta21(ie),beta22(ie),beta23(ie),&
220 beta31(ie),beta32(ie),beta33(ie),&
221 gamma1(ie),gamma2(ie),gamma3(ie),&
222 delta1(ie),delta2(ie),delta3(ie),&
223 ux_el,uy_el,uz_el,&
224 duxdx_el,duydx_el,duzdx_el,&
225 duxdy_el,duydy_el,duzdy_el,&
226 duxdz_el,duydz_el,duzdz_el)
227
228
229
230 lambda_el = lambda; mu_el = mu;
231
232
233 call make_stress_tensor(nn,lambda_el,mu_el,&
234 duxdx_el,duydx_el,duzdx_el,&
235 duxdy_el,duydy_el,duzdy_el,&
236 duxdz_el,duydz_el,duzdz_el,&
237 sxx_el,syy_el,szz_el,&
238 syz_el,szx_el,sxy_el)
239
240
241 do k = 1,mm
242 do j = 1,mm
243 do i = 1,mm
244 dxdx = alfa11(ie) +beta12(ie)*ctm(k) &
245 + beta13(ie)*ctm(j) +gamma1(ie)*ctm(j)*ctm(k)
246 dydx = alfa21(ie) +beta22(ie)*ctm(k) &
247 + beta23(ie)*ctm(j) +gamma2(ie)*ctm(j)*ctm(k)
248 dzdx = alfa31(ie) +beta32(ie)*ctm(k) &
249 + beta33(ie)*ctm(j) +gamma3(ie)*ctm(j)*ctm(k)
250
251 dxdy = alfa12(ie) +beta11(ie)*ctm(k) &
252 + beta13(ie)*ctm(i) +gamma1(ie)*ctm(k)*ctm(i)
253 dydy = alfa22(ie) +beta21(ie)*ctm(k) &
254 + beta23(ie)*ctm(i) +gamma2(ie)*ctm(k)*ctm(i)
255 dzdy = alfa32(ie) +beta31(ie)*ctm(k) &
256 + beta33(ie)*ctm(i) +gamma3(ie)*ctm(k)*ctm(i)
257
258 dxdz = alfa13(ie) +beta11(ie)*ctm(j) &
259 + beta12(ie)*ctm(i) +gamma1(ie)*ctm(i)*ctm(j)
260 dydz = alfa23(ie) +beta21(ie)*ctm(j) &
261 + beta22(ie)*ctm(i) +gamma2(ie)*ctm(i)*ctm(j)
262 dzdz = alfa33(ie) +beta31(ie)*ctm(j) &
263 + beta32(ie)*ctm(i) +gamma3(ie)*ctm(i)*ctm(j)
264
265 det_j = dxdz * (dydx*dzdy - dzdx*dydy) &
266 - dydz * (dxdx*dzdy - dzdx*dxdy) &
267 + dzdz * (dxdx*dydy - dydx*dxdy)
268
269 !Local to global node transformation
270 xp = alfa11(ie)*ctm(i) + alfa12(ie)*ctm(j) &
271 + alfa13(ie)*ctm(k) + beta11(ie)*ctm(j)*ctm(k) &
272 + beta12(ie)*ctm(i)*ctm(k) + beta13(ie)*ctm(i)*ctm(j) &
273 + gamma1(ie)*ctm(i)*ctm(j)*ctm(k) + delta1(ie)
274
275 yp = alfa21(ie)*ctm(i) + alfa22(ie)*ctm(j) &
276 + alfa23(ie)*ctm(k) + beta21(ie)*ctm(j)*ctm(k) &
277 + beta22(ie)*ctm(i)*ctm(k) + beta23(ie)*ctm(i)*ctm(j) &
278 + gamma2(ie)*ctm(i)*ctm(j)*ctm(k) + delta2(ie)
279
280 zp = alfa31(ie)*ctm(i) + alfa32(ie)*ctm(j) &
281 + alfa33(ie)*ctm(k) + beta31(ie)*ctm(j)*ctm(k) &
282 + beta32(ie)*ctm(i)*ctm(k) + beta33(ie)*ctm(i)*ctm(j) &
283 + gamma3(ie)*ctm(i)*ctm(j)*ctm(k) + delta3(ie)
284
285 !STATIONARY POLYNOMIAL SOLUTIONS
286 !u1_ex = dsin(3.d0*pi*time)*((xp*(xp-1.d0))*(yp*(yp-1.d0))*(zp*(zp-1.d0)));
287 !u2_ex = 0.d0;
288 !u3_ex = 0.d0;
289
290
291 is = nn*nn*(k -1) +nn*(j -1) +i
292 in = cs_loc(cs_loc(ie -1) +is)
293
294 !POPLYNOMIAL SOLUTION
295 !u1_ex = dsin(3.d0*pi*time)*((xs_loc(in)*(xs_loc(in)-1.d0))**2.d0 &
296 ! *(ys_loc(in)*(ys_loc(in)-1.d0))**2.d0 &
297 ! *(zs_loc(in)*(zs_loc(in)-1.d0))**2.d0);
298 !u2_ex = dsin(3.d0*pi*time)*((xs_loc(in)*(xs_loc(in)-1.d0))**2.d0 &
299 ! *(ys_loc(in)*(ys_loc(in)-1.d0))**2.d0 &
300 ! *(zs_loc(in)*(zs_loc(in)-1.d0))**2.d0);
301 !u3_ex = dsin(3.d0*pi*time)*((xs_loc(in)*(xs_loc(in)-1.d0))**2.d0 &
302 ! *(ys_loc(in)*(ys_loc(in)-1.d0))**2.d0 &
303 ! *(zs_loc(in)*(zs_loc(in)-1.d0))**2.d0);
304 !
305 !v1_ex = + 3.d0*pi*dcos(3.d0*pi*time)*((xs_loc(in)*(xs_loc(in)-1.d0))**2.d0 &
306 ! *(ys_loc(in)*(ys_loc(in)-1.d0))**2.d0 &
307 ! *(zs_loc(in)*(zs_loc(in)-1.d0))**2.d0);
308 !v2_ex = + 3.d0*pi*dcos(3.d0*pi*time)*((xs_loc(in)*(xs_loc(in)-1.d0))**2.d0 &
309 ! *(ys_loc(in)*(ys_loc(in)-1.d0))**2.d0 &
310 ! *(zs_loc(in)*(zs_loc(in)-1.d0))**2.d0);
311 !v3_ex = + 3.d0*pi*dcos(3.d0*pi*time)*((xs_loc(in)*(xs_loc(in)-1.d0))**2.d0 &
312 ! *(ys_loc(in)*(ys_loc(in)-1.d0))**2.d0 &
313 ! *(zs_loc(in)*(zs_loc(in)-1.d0))**2.d0);
314
315
316
317
318 !PAPER WITH BLANCA
319 u1_ex = -dsin(3.d0*pi*time)*(dsin(pi*xp))**2.d0*dsin(2.d0*pi*yp)*dsin(2.d0*pi*zp);
320 u2_ex = +dsin(3.d0*pi*time)*(dsin(pi*yp))**2.d0*dsin(2.d0*pi*xp)*dsin(2.d0*pi*zp);
321 u3_ex = +dsin(3.d0*pi*time)*(dsin(pi*zp))**2.d0*dsin(2.d0*pi*xp)*dsin(2.d0*pi*yp);
322
323 v1_ex = - 3.d0*pi*dcos(3.d0*pi*time)*(dsin(pi*xs_loc(in)))**2.d0*dsin(2.d0*pi*ys_loc(in))*dsin(2.d0*pi*zs_loc(in));
324 v2_ex = + 3.d0*pi*dcos(3.d0*pi*time)*(dsin(pi*ys_loc(in)))**2.d0*dsin(2.d0*pi*xs_loc(in))*dsin(2.d0*pi*zs_loc(in));
325 v3_ex = + 3.d0*pi*dcos(3.d0*pi*time)*(dsin(pi*zs_loc(in)))**2.d0*dsin(2.d0*pi*xs_loc(in))*dsin(2.d0*pi*ys_loc(in));
326
327 x = xs_loc(in); y = ys_loc(in); z = zs_loc(in);
328 e11_ex = -2.d0*pi*dcos(pi*x)*dsin(3.d0*pi*time)*dsin(pi*x)*dsin(2.d0*pi*y)*dsin(2.d0*pi*z)
329 e22_ex = 2.d0*pi*dcos(pi*y)*dsin(3.d0*pi*time)*dsin(2.d0*pi*x)*dsin(pi*y)*dsin(2.d0*pi*z)
330
331 e33_ex = 2.d0*pi*dcos(pi*z)*dsin(3.d0*pi*time)*dsin(2.d0*pi*x)*dsin(2.d0*pi*y)*dsin(pi*z)
332
333 e12_ex = pi*dcos(2.d0*pi*x)*dsin(3.d0*pi*time)*dsin(pi*y)**2.d0*dsin(2.d0*pi*z) &
334 - pi*dcos(2.d0*pi*y)*dsin(3.d0*pi*time)*dsin(pi*x)**2.d0*dsin(2.d0*pi*z)
335
336 e13_ex = pi*dcos(2.d0*pi*x)*dsin(3.d0*pi*time)*dsin(2.d0*pi*y)*dsin(pi*z)**2.d0 &
337 - pi*dcos(2.d0*pi*z)*dsin(3.d0*pi*time)*dsin(pi*x)**2.d0*dsin(2.d0*pi*y)
338
339 e23_ex = pi*dcos(2.d0*pi*z)*dsin(3.d0*pi*time)*dsin(2.d0*pi*x)*dsin(pi*y)**2.d0 &
340 + pi*dcos(2.d0*pi*y)*dsin(3.d0*pi*time)*dsin(2.d0*pi*x)*dsin(pi*z)**2.d0
341
342
343 sxx_ex = lambda*(2.d0*pi*dcos(pi*y)*dsin(3.d0*pi*time)*dsin(2.d0*pi*x)*dsin(pi*y)*dsin(2.d0*pi*z) &
344 - 2.d0*pi*dcos(pi*x)*dsin(3.d0*pi*time)*dsin(pi*x)*dsin(2.d0*pi*y)*dsin(2.d0*pi*z) &
345 + 2.d0*pi*dcos(pi*z)*dsin(3.d0*pi*time)*dsin(2.d0*pi*x)*dsin(2.d0*pi*y)*dsin(pi*z)) &
346 - 4.d0*pi*mu*dcos(pi*x)*dsin(3.d0*pi*time)*dsin(pi*x)*dsin(2.d0*pi*y)*dsin(2.d0*pi*z)
347
348 syy_ex = lambda*(2.d0*pi*dcos(pi*y)*dsin(3.d0*pi*time)*dsin(2.d0*pi*x)*dsin(pi*y)*dsin(2.d0*pi*z) &
349 - 2.d0*pi*dcos(pi*x)*dsin(3.d0*pi*time)*dsin(pi*x)*dsin(2.d0*pi*y)*dsin(2.d0*pi*z) &
350 + 2.d0*pi*dcos(pi*z)*dsin(3.d0*pi*time)*dsin(2.d0*pi*x)*dsin(2.d0*pi*y)*dsin(pi*z)) &
351 + 4.d0*pi*mu*dcos(pi*y)*dsin(3.d0*pi*time)*dsin(2.d0*pi*x)*dsin(pi*y)*dsin(2.d0*pi*z)
352
353 szz_ex = lambda*(2*pi*dcos(pi*y)*dsin(3.d0*pi*time)*dsin(2.d0*pi*x)*dsin(pi*y)*dsin(2.d0*pi*z) &
354 - 2.d0*pi*dcos(pi*x)*dsin(3.d0*pi*time)*dsin(pi*x)*dsin(2.d0*pi*y)*dsin(2.d0*pi*z) &
355 + 2.d0*pi*dcos(pi*z)*dsin(3.d0*pi*time)*dsin(2.d0*pi*x)*dsin(2.d0*pi*y)*dsin(pi*z)) &
356 + 4.d0*pi*mu*dcos(pi*z)*dsin(3.d0*pi*time)*dsin(2.d0*pi*x)*dsin(2.d0*pi*y)*dsin(pi*z)
357
358 sxy_ex = mu*(2.d0*pi*dcos(2.d0*pi*x)*dsin(3.d0*pi*time)*dsin(pi*y)**2.d0*dsin(2.d0*pi*z) &
359 - 2.d0*pi*dcos(2.d0*pi*y)*dsin(3.d0*pi*time)*dsin(pi*x)**2.d0*dsin(2.d0*pi*z))
360
361 szx_ex = mu*(2.d0*pi*dcos(2.d0*pi*x)*dsin(3.d0*pi*time)*dsin(2.d0*pi*y)*dsin(pi*z)**2.d0 &
362 - 2.d0*pi*dcos(2.d0*pi*z)*dsin(3.d0*pi*time)*dsin(pi*x)**2.d0*dsin(2.d0*pi*y))
363
364 syz_ex = mu*(2.d0*pi*dcos(2.d0*pi*z)*dsin(3.d0*pi*time)*dsin(2.d0*pi*x)*dsin(pi*y)**2.d0 &
365 + 2.d0*pi*dcos(2.d0*pi*y)*dsin(3.d0*pi*time)*dsin(2.d0*pi*x)*dsin(pi*z)**2.d0)
366
367
368 u1_ctm = 0.d0; u2_ctm = 0.d0; u3_ctm = 0.d0
369
370 allocate(basis_functions(mm**3)) ! basis function evaluation
371 do kb = 1, nn
372 do jb = 1, nn
373 do ib =1, nn
374 isb = nn*nn*(kb -1) + nn*(jb -1) +ib
375 in = cs_loc(cs_loc(ie -1) +isb)
376 iaz = 3*(in -1) + 1
377
378 basis_functions(isb) = phi(nn-1, ct(ib), ctm(i)) * &
379 phi(nn-1, ct(jb), ctm(j)) * &
380 phi(nn-1, ct(kb), ctm(k))
381
382 u1_ctm = u1_ctm + basis_functions(isb)*u1(iaz);
383 u2_ctm = u2_ctm + basis_functions(isb)*u1(iaz+1);
384 u3_ctm = u3_ctm + basis_functions(isb)*u1(iaz+2);
385
386 enddo
387 enddo
388 enddo
389 deallocate(basis_functions)
390
391
392
393
394
395 !L2 NORM OLD
396 is = nn*nn*(k -1) +nn*(j -1) +i
397 in = cs_loc(cs_loc(ie -1) +is)
398
399 iaz = 3*(in -1) + 1
400
401 term_l2_vel = det_j * ww(i) * ww(j) * ww(k) * &
402 ( (v1_ex - v1(iaz))**2.d0 + (v2_ex - v1(iaz+1))**2.d0 +(v3_ex - v1(iaz+2))**2.d0)
403
404
405 l2_err_vel = l2_err_vel + term_l2_vel
406
407
408
409 !LINF NORM - DISPL & VEL
410 linf_array = (/ linf_err,abs(u1_ex - u1_ctm),abs(u2_ex - u2_ctm),abs(u3_ex - u3_ctm) /)
411 linf_err = maxval(linf_array)
412
413
414 !L2 NORM NEW
415 term_l2 = det_j * wwm(i) * wwm(j) * wwm(k) * &
416 ( (u1_ex - u1_ctm)**2.d0 + (u2_ex - u2_ctm)**2.d0 +(u3_ex - u3_ctm)**2.d0)
417 l2_err = l2_err + term_l2
418
419
420
421 !ENERGY NORM
422 iaz = 3*(in -1) + 1
423 term_ener = det_j * ww(i) * ww(j) * ww(k) * ( &
424 rho * (v1_ex - v1(iaz+0))**2.d0 &
425 + rho * (v2_ex - v1(iaz+1))**2.d0 &
426 + rho * (v3_ex - v1(iaz+2))**2.d0 &
427 + (sxx_el(i,j,k) - sxx_ex) * (duxdx_el(i,j,k) - e11_ex) &
428 + (syy_el(i,j,k) - syy_ex) * (duydy_el(i,j,k) - e22_ex) &
429 + (szz_el(i,j,k) - szz_ex) * (duzdz_el(i,j,k) - e33_ex) &
430 + 2.d0*(sxy_el(i,j,k) - sxy_ex) * (0.5d0*duxdy_el(i,j,k) + 0.5d0*duydx_el(i,j,k) - e12_ex) &
431 + 2.d0*(szx_el(i,j,k) - szx_ex) * (0.5d0*duxdz_el(i,j,k) + 0.5d0*duzdx_el(i,j,k) - e13_ex) &
432 + 2.d0*(syz_el(i,j,k) - syz_ex) * (0.5d0*duydz_el(i,j,k) + 0.5d0*duzdy_el(i,j,k) - e23_ex) &
433 )
434
435
436 en_err = en_err + term_ener
437
438 term_h1 = det_j * ww(i) * ww(j) * ww(k) * ( &
439 + (sxx_el(i,j,k) - sxx_ex) * (duxdx_el(i,j,k) - e11_ex) &
440 + (syy_el(i,j,k) - syy_ex) * (duydy_el(i,j,k) - e22_ex) &
441 + (szz_el(i,j,k) - szz_ex) * (duzdz_el(i,j,k) - e33_ex) &
442 + 2.d0*(sxy_el(i,j,k) - sxy_ex) * (0.5d0*duxdy_el(i,j,k) + 0.5d0*duydx_el(i,j,k) - e12_ex) &
443 + 2.d0*(szx_el(i,j,k) - szx_ex) * (0.5d0*duxdz_el(i,j,k) + 0.5d0*duzdx_el(i,j,k) - e13_ex) &
444 + 2.d0*(syz_el(i,j,k) - syz_ex) * (0.5d0*duydz_el(i,j,k) + 0.5d0*duzdy_el(i,j,k) - e23_ex) &
445 )
446
447 h1_err = h1_err + term_h1
448
449 enddo
450 enddo
451 enddo
452
453
454 deallocate(ct,ww,dd, ctm, wwm, ddm)
455 deallocate(ux_el,uy_el,uz_el)
456 deallocate(duxdx_el,duydx_el,duzdx_el,duxdy_el,duydy_el,duzdy_el,duxdz_el,duydz_el,duzdz_el)
457 deallocate(sxx_el,syy_el,szz_el,syz_el,szx_el,sxy_el)
458 deallocate(lambda_el,mu_el)
459
460
461 enddo
462
463
464
465!----------------------------------------------------------------------------------------------------------
466! NOW COMPUTE NON CONFORMING TERM S_e int_e [u-uh][u-uh]
467!----------------------------------------------------------------------------------------------------------
468 if(nelem_dg_global .gt. 0) then
469
470 jump_minus_error = 0.d0
471
472 do i = 1,nsd_jump
473 in = node_sd_jump(i)
474 !PAPER WITH BLANCA
475 u1_ex = -dsin(3.d0*pi*time)*(dsin(pi*xs_loc(in)))**2.d0*dsin(2.d0*pi*ys_loc(in))*dsin(2.d0*pi*zs_loc(in));
476 u2_ex = +dsin(3.d0*pi*time)*(dsin(pi*ys_loc(in)))**2.d0*dsin(2.d0*pi*xs_loc(in))*dsin(2.d0*pi*zs_loc(in));
477 u3_ex = +dsin(3.d0*pi*time)*(dsin(pi*zs_loc(in)))**2.d0*dsin(2.d0*pi*xs_loc(in))*dsin(2.d0*pi*ys_loc(in));
478
479
480 send_buffer_error(3*(i -1) +1) = abs(u1(3*(in -1) +1) - u1_ex)
481 send_buffer_error(3*(i -1) +2) = abs(u1(3*(in -1) +2) - u2_ex)
482 send_buffer_error(3*(i -1) +3) = abs(u1(3*(in -1) +3) - u3_ex)
483
484 enddo
485
486 call exchange_double(3*nsd_jump, send_buffer_error, 3*nrv_jump, recv_buffer_error,&
487 mpi_np, send_length_jump, recv_length_jump,&
488 mpi_comm, mpi_stat, mpi_ierr,mpi_id)
489
490 do i = 1,nrv_jump
491 in = node_rv_jump(i)
492 jump_minus_error(3*(i -1) +1) = recv_buffer_error(3*(in -1) +1)
493 jump_minus_error(3*(i -1) +2) = recv_buffer_error(3*(in -1) +2)
494 jump_minus_error(3*(i -1) +3) = recv_buffer_error(3*(in -1) +3)
495 enddo
496
497 endif
498
499
500
501 do ie = 1, nelem_dg_local
502
503 ielem = el_new(ie)%ind; nn = el_new(ie)%deg
504 call get_indloc_from_indglo(local_el_num, ne_loc, ielem, ie_curr)
505
506 allocate(u_p(3*nn**3)); u_p = 0.d0;
507
508 do k = 1,nn
509 do j = 1,nn
510 do i = 1,nn
511 is = nn*nn*(k -1) +nn*(j -1) +i
512 in = cs_loc(cs_loc(ie_curr -1) + is)
513
514 !PAPER WITH BLANCA
515 u1_ex = -dsin(3.d0*pi*time)*(dsin(pi*xs_loc(in)))**2.d0*dsin(2.d0*pi*ys_loc(in))*dsin(2.d0*pi*zs_loc(in));
516 u2_ex = +dsin(3.d0*pi*time)*(dsin(pi*ys_loc(in)))**2.d0*dsin(2.d0*pi*xs_loc(in))*dsin(2.d0*pi*zs_loc(in));
517 u3_ex = +dsin(3.d0*pi*time)*(dsin(pi*zs_loc(in)))**2.d0*dsin(2.d0*pi*xs_loc(in))*dsin(2.d0*pi*ys_loc(in));
518
519
520 iaz = 3*(in -1) +1; u_p(is) = abs(u1(iaz) - u1_ex)
521 iaz = 3*(in -1) +2; u_p(nn**3+is) = abs(u1(iaz) - u2_ex)
522 iaz = 3*(in -1) +3; u_p(2*nn**3 + is) = abs(u1(iaz) - u3_ex)
523
524
525 enddo
526 enddo
527 enddo
528
529
530 allocate(u_m(el_new(ie)%nnz_col_only_uv)); u_m = 0.d0; jshift = 0
531
532 do ic = 1, el_new(ie)%num_of_ne
533
534 iene = el_new(ie)%el_conf(ic,1)
535 imne = el_new(ie)%el_conf(ic,0)
536
537 mm = 2
538 do i = 1, nm
539 if(tag_mat(i) .eq. imne ) then
540 mm = sdeg_mat(i) +1
541 endif
542 enddo
543
544 trofa = 0
545 if(mpi_np .gt. 1) call check_mpi(cs_nnz_dg, cs_dg, iene, trofa, posit)
546
547 if(trofa .eq. 1) then
548 do k = 1,mm
549 do j = 1,mm
550 do i = 1,mm
551 is = mm*mm*(k -1) +mm*(j -1) +i
552 id = cs_dg(cs_dg(posit-1)+is)
553 iaz = 3*(id -1) +1; u_m(is + jshift) = jump_minus_error(iaz)
554 iaz = 3*(id -1) +2; u_m(mm**3 + is + jshift) = jump_minus_error(iaz)
555 iaz = 3*(id -1) +3; u_m(2*(mm**3) + is + jshift) = jump_minus_error(iaz)
556 enddo
557 enddo
558 enddo
559 else
560
561 call get_indloc_from_indglo(local_el_num, ne_loc, iene, iene_curr)
562
563 do k = 1,mm
564 do j = 1,mm
565 do i = 1,mm
566 is = mm*mm*(k -1) +mm*(j -1) +i
567 id = cs_loc(cs_loc(iene_curr-1)+is)
568
569 !PAPER WITH BLANCA
570 u1_ex = -dsin(3.d0*pi*time)*(dsin(pi*xs_loc(id)))**2.d0*dsin(2.d0*pi*ys_loc(id))*dsin(2.d0*pi*zs_loc(id));
571 u2_ex = +dsin(3.d0*pi*time)*(dsin(pi*ys_loc(id)))**2.d0*dsin(2.d0*pi*xs_loc(id))*dsin(2.d0*pi*zs_loc(id));
572 u3_ex = +dsin(3.d0*pi*time)*(dsin(pi*zs_loc(id)))**2.d0*dsin(2.d0*pi*xs_loc(id))*dsin(2.d0*pi*ys_loc(id));
573
574
575 iaz = 3*(id -1) +1; u_m(is + jshift) = abs(u1(iaz) - u1_ex)
576 iaz = 3*(id -1) +2; u_m(mm**3 + is + jshift) = abs(u1(iaz) - u2_ex)
577 iaz = 3*(id -1) +3; u_m(2*(mm**3) + is + jshift) = abs(u1(iaz) - u3_ex)
578
579
580 enddo
581 enddo
582 enddo
583 endif
584
585 ! write(*,*) '--------------- fine +- ------------------------------'
586 ! read(*,*)
587
588
589 jshift = jshift + 3*mm**3
590
591 enddo
592
593 allocate(jump_all(3*nn**3))
594 jump_all = 0.d0
595
596
597 call matmul_sparse(el_new(ie)%matMin_only_uv, el_new(ie)%nnz_minus_only_uv, el_new(ie)%JMin_only_uv, &
598 el_new(ie)%IMin_only_uv, jump_all, 3*nn**3, &
599 u_m, el_new(ie)%nnz_col_only_uv,1)
600
601
602 en_err = en_err + dot_product(jump_all,u_p)
603 h1_err = h1_err + dot_product(jump_all,u_p)
604
605
606 deallocate(u_m, jump_all)
607 allocate(jump_all(3*nn**3)); jump_all = 0.d0
608
609 call matmul_sparse(el_new(ie)%matPlus_only_uv, el_new(ie)%nnz_plus_only_uv,el_new(ie)%JPlus_only_uv, &
610 el_new(ie)%IPlus_only_uv, jump_all, 3*nn**3, &
611 u_p, 3*nn**3,1)
612
613
614 en_err = en_err + dot_product(jump_all,u_p)
615 h1_err = h1_err + dot_product(jump_all,u_p)
616
617
618 deallocate(u_p, jump_all)
619
620 enddo
621
622
623!----------------------------------------------------------------------------------------------------------
624! NOW COMPUTE NON CONFORMING TERM S_e int_e [v-vh][v-vh] v=du/dt
625!----------------------------------------------------------------------------------------------------------
626 if(nelem_dg_global .gt. 0) then
627
628 jump_minus_error = 0.d0
629
630 do i = 1,nsd_jump
631 in = node_sd_jump(i)
632 !PAPER WITH BLANCA
633 v1_ex = - 3.d0*pi*dcos(3.d0*pi*time)*(dsin(pi*xs_loc(in)))**2.d0*dsin(2.d0*pi*ys_loc(in))*dsin(2.d0*pi*zs_loc(in));
634 v2_ex = + 3.d0*pi*dcos(3.d0*pi*time)*(dsin(pi*ys_loc(in)))**2.d0*dsin(2.d0*pi*xs_loc(in))*dsin(2.d0*pi*zs_loc(in));
635 v3_ex = + 3.d0*pi*dcos(3.d0*pi*time)*(dsin(pi*zs_loc(in)))**2.d0*dsin(2.d0*pi*xs_loc(in))*dsin(2.d0*pi*ys_loc(in));
636
637
638 send_buffer_error(3*(i -1) +1) = abs(v1(3*(in -1) +1) - v1_ex)
639 send_buffer_error(3*(i -1) +2) = abs(v1(3*(in -1) +2) - v2_ex)
640 send_buffer_error(3*(i -1) +3) = abs(v1(3*(in -1) +3) - v3_ex)
641
642 enddo
643
644 call exchange_double(3*nsd_jump, send_buffer_error, 3*nrv_jump, recv_buffer_error,&
645 mpi_np, send_length_jump, recv_length_jump,&
646 mpi_comm, mpi_stat, mpi_ierr,mpi_id)
647
648 do i = 1,nrv_jump
649 in = node_rv_jump(i)
650 jump_minus_error(3*(i -1) +1) = recv_buffer_error(3*(in -1) +1)
651 jump_minus_error(3*(i -1) +2) = recv_buffer_error(3*(in -1) +2)
652 jump_minus_error(3*(i -1) +3) = recv_buffer_error(3*(in -1) +3)
653 enddo
654
655 endif
656
657
658 do ie = 1, nelem_dg_local
659
660 ielem = el_new(ie)%ind; nn = el_new(ie)%deg
661 call get_indloc_from_indglo(local_el_num, ne_loc, ielem, ie_curr)
662
663 allocate(u_p(3*nn**3)); u_p = 0.d0;
664
665 do k = 1,nn
666 do j = 1,nn
667 do i = 1,nn
668 is = nn*nn*(k -1) +nn*(j -1) +i
669 in = cs_loc(cs_loc(ie_curr -1) + is)
670
671 !PAPER WITH BLANCA
672 v1_ex = - 3.d0*pi*dcos(3.d0*pi*time)*(dsin(pi*xs_loc(in)))**2.d0*dsin(2.d0*pi*ys_loc(in))*dsin(2.d0*pi*zs_loc(in));
673 v2_ex = + 3.d0*pi*dcos(3.d0*pi*time)*(dsin(pi*ys_loc(in)))**2.d0*dsin(2.d0*pi*xs_loc(in))*dsin(2.d0*pi*zs_loc(in));
674 v3_ex = + 3.d0*pi*dcos(3.d0*pi*time)*(dsin(pi*zs_loc(in)))**2.d0*dsin(2.d0*pi*xs_loc(in))*dsin(2.d0*pi*ys_loc(in));
675
676
677 iaz = 3*(in -1) +1; u_p(is) = abs(v1(iaz) - v1_ex)
678 iaz = 3*(in -1) +2; u_p(nn**3+is) = abs(v1(iaz) - v2_ex)
679 iaz = 3*(in -1) +3; u_p(2*nn**3 + is) = abs(v1(iaz) - v3_ex)
680
681
682 enddo
683 enddo
684 enddo
685
686
687 allocate(u_m(el_new(ie)%nnz_col_only_uv)); u_m = 0.d0; jshift = 0
688
689 do ic = 1, el_new(ie)%num_of_ne
690
691 iene = el_new(ie)%el_conf(ic,1)
692 imne = el_new(ie)%el_conf(ic,0)
693
694 mm = 2
695 do i = 1, nm
696 if(tag_mat(i) .eq. imne ) then
697 mm = sdeg_mat(i) +1
698 endif
699 enddo
700
701 trofa = 0
702 if(mpi_np .gt. 1) call check_mpi(cs_nnz_dg, cs_dg, iene, trofa, posit)
703
704 if(trofa .eq. 1) then
705 do k = 1,mm
706 do j = 1,mm
707 do i = 1,mm
708 is = mm*mm*(k -1) +mm*(j -1) +i
709 id = cs_dg(cs_dg(posit-1)+is)
710 iaz = 3*(id -1) +1; u_m(is + jshift) = jump_minus_error(iaz)
711 iaz = 3*(id -1) +2; u_m(mm**3 + is + jshift) = jump_minus_error(iaz)
712 iaz = 3*(id -1) +3; u_m(2*(mm**3) + is + jshift) = jump_minus_error(iaz)
713 enddo
714 enddo
715 enddo
716 else
717
718 call get_indloc_from_indglo(local_el_num, ne_loc, iene, iene_curr)
719
720 do k = 1,mm
721 do j = 1,mm
722 do i = 1,mm
723 is = mm*mm*(k -1) +mm*(j -1) +i
724 id = cs_loc(cs_loc(iene_curr-1)+is)
725
726 !PAPER WITH BLANCA
727 v1_ex = - 3.d0*pi*dcos(3.d0*pi*time)*(dsin(pi*xs_loc(id)))**2.d0*dsin(2.d0*pi*ys_loc(id))*dsin(2.d0*pi*zs_loc(id));
728 v2_ex = + 3.d0*pi*dcos(3.d0*pi*time)*(dsin(pi*ys_loc(id)))**2.d0*dsin(2.d0*pi*xs_loc(id))*dsin(2.d0*pi*zs_loc(id));
729 v3_ex = + 3.d0*pi*dcos(3.d0*pi*time)*(dsin(pi*zs_loc(id)))**2.d0*dsin(2.d0*pi*xs_loc(id))*dsin(2.d0*pi*ys_loc(id));
730
731
732
733 iaz = 3*(id -1) +1; u_m(is + jshift) = abs(v1(iaz) - v1_ex)
734 iaz = 3*(id -1) +2; u_m(mm**3 + is + jshift) = abs(v1(iaz) - v2_ex)
735 iaz = 3*(id -1) +3; u_m(2*(mm**3) + is + jshift) = abs(v1(iaz) - v3_ex)
736
737
738 enddo
739 enddo
740 enddo
741 endif
742
743 ! write(*,*) '--------------- fine +- ------------------------------'
744 ! read(*,*)
745
746
747 jshift = jshift + 3*mm**3
748
749 enddo
750
751 allocate(jump_all(3*nn**3))
752 jump_all = 0.d0
753
754
755 call matmul_sparse(el_new(ie)%matMin_only_uv, el_new(ie)%nnz_minus_only_uv, el_new(ie)%JMin_only_uv, &
756 el_new(ie)%IMin_only_uv, jump_all, 3*nn**3, &
757 u_m, el_new(ie)%nnz_col_only_uv,1)
758
759
760 en_err = en_err + dot_product(jump_all,u_p)
761 !write(*,*) dot_product(jump_all,u_p)
762 !H1_err = H1_err + dot_product(jump_all,u_p)
763
764
765 deallocate(u_m, jump_all)
766 allocate(jump_all(3*nn**3)); jump_all = 0.d0
767
768 call matmul_sparse(el_new(ie)%matPlus_only_uv, el_new(ie)%nnz_plus_only_uv,el_new(ie)%JPlus_only_uv, &
769 el_new(ie)%IPlus_only_uv, jump_all, 3*nn**3, &
770 u_p, 3*nn**3,1)
771
772
773 en_err = en_err + dot_product(jump_all,u_p)
774 !write(*,*) dot_product(jump_all,u_p)
775 !read(*,*)
776 !H1_err = H1_err + dot_product(jump_all,u_p)
777
778
779 deallocate(u_p, jump_all)
780
781 enddo
782
783
784
785
786 linf_err_tot = dabs(linf_err)
787 !Linf_err_vel_tot = dabs(Linf_err_vel)
788
789 l2_err_tot = dabs(l2_err)
790 l2_err_vel_tot = dabs(l2_err_vel)
791 h1_err_tot = dabs(h1_err)
792 en_err_tot = dabs(en_err)
793
794 if (en_err_tot .ge. 100.d0) call exit(exit_energy_error)
795
796 if(mpi_id .eq. 0) then
797 file_mpi = 'EN.ERR'; unit_mpi = 40;
798 open(unit_mpi,file=file_mpi, position='APPEND')
799 write(unit_mpi,*) sqrt(h1_err_tot), sqrt(en_err_tot)!, sqrt(L2_err_tot),sqrt(L2_err_vel_tot)
800
801 close(unit_mpi)
802 endif
803
804
805
806
807
808
809
810
811
812 end subroutine compute_energy_error
check_mpi
subroutine check_mpi(n, v, ind, tt, pos)
Checks if a an element is present on a vector and give its position.
Definition CHECK_MPI.f90:30
compute_energy_error
subroutine compute_energy_error(nnod_loc, u1, v1, time, ne_loc, cs_loc, cs_nnz_loc, nm, prop_mat, sdeg_mat, tag_mat, loc_n_num, local_el_num, nelem_dg_local, alfa11, alfa12, alfa13, alfa21, alfa22, alfa23, alfa31, alfa32, alfa33, beta11, beta12, beta13, beta21, beta22, beta23, beta31, beta32, beta33, gamma1, gamma2, gamma3, delta1, delta2, delta3, xs_loc, ys_loc, zs_loc, mpi_id, el_new, nelem_dg_global, nsd_jump, node_sd_jump, nrv_jump, node_rv_jump, send_length_jump, recv_length_jump, mpi_np, mpi_comm, cs_nnz_dg, cs_dg)
Computes norms of error in test mode case.
Definition COMPUTE_ENERGY_ERROR.f90:100
exchange_double
subroutine exchange_double(nsend, buff_send, nrecv, buff_recv, nproc, proc_send, proc_recv, comm, status, ierr, myid)
Exchanges double between MPI processes.
Definition EXCHANGE_DOUBLE.f90:39
get_indloc_from_indglo
subroutine get_indloc_from_indglo(local_el, nel_loc, ie, ic)
Returns local id from global id.
Definition GET_INDLOC_FROM_INDGLO.f90:30
make_lgl_nw
subroutine make_lgl_nw(nb_pnt, xq, wq, dd)
Makes Gauss-Legendre-Lobatto nodes, weigths and spectral derivatives.
Definition MAKE_LGL_NW.f90:29
make_strain_tensor
subroutine make_strain_tensor(nn, ct, ww, dd, alfa11, alfa12, alfa13, alfa21, alfa22, alfa23
Computes spatial derivatives of displacement.
Definition MAKE_STRAIN_TENSOR.f90:67
make_stress_tensor
subroutine make_stress_tensor(nn, lambda, mu, duxdx, duydx, duzdx, duxdy, duydy, duzdy, duxdz, duydz, duzdz, sxx, syy, szz, syz, szx, sxy)
Computes the stress tensor.
Definition MAKE_STRESS_TENSOR.f90:49
matmul_sparse
subroutine matmul_sparse(as, nnz, jsp, isp, vet_out, n, vet_in,
Matrix-vector multiplication for sparse matrices (RCS format).
Definition MATMUL_SPARSE.f90:34
phi
real *8 function phi(ng, csii, csij)
Evaluates Lengendre basis function on a specific point.
Definition PHI.f90:31
dgjump
Contains structure for jump matrices.
Definition MODULES.f90:155
max_var
Set maximal bounds.
Definition MODULES.f90:54
speed_exit_codes
SPEED exit codes.
Definition MODULES.f90:25
speed_exit_codes::exit_energy_error
integer, parameter exit_energy_error
Definition MODULES.f90:36
dgjump::el4loop
element structure for time loop (RAM saving)
Definition MODULES.f90:167