ETS  \$Id: Doxyfile 2162 2020-02-26 14:16:09Z g2dpc $
 All Classes Files Functions Variables Pages
solution6.f90
Go to the documentation of this file.
1 !-----------------------------------------------------------------------
2 ! GTS transport solver (1D stiff convection-diffusion-reaction solver)
3 !
4 ! - using cubic Hermite finite elements
5 ! - hyper diffusivity as stabilisation
6 ! - Crank-Nicholson or Gears scheme as time evolution scheme
7 !
8 ! Author : Guido Huysmans (Association Euratom-CEA)
9 ! version : 0.95
10 ! Date : 21-5-2009
11 !-----------------------------------------------------------------------
12 module gts_parameters
13  implicit none
14  integer :: neq, norder, nvar, gts_status ! neq is the number of equations
15  parameter(norder=3) ! the order of the finite elements
16  parameter(nvar=2) ! the number of variables per node
17 endmodule
18 
19 module gts_matrix
20  implicit none
21  real*8, ALLOCATABLE :: AA(:,:) ! fem matrix
22  real*8, ALLOCATABLE :: BB(:) ! right hand side
23  integer :: LDA, N, KL, KU ! LDA : the leading dimension of matrix A
24 endmodule ! N : the dimension of A and B
25  ! KL,KU : the number of non-zero diagonals below,above the main diagonal
26 module gts_grid
27  implicit none
28  real*8, ALLOCATABLE :: x_grid(:,:),xg(:) ! the positions of the gridpoints and gaussian points
29  integer :: n_grid ! the number of grid points
30 endmodule
31 
32 module gts_gauss ! contains the Gaussain quadrature points and weights
33  implicit none
34  integer,parameter :: Ngauss = 4 ! the order of the integration
35  real*8,parameter :: Xgauss(Ngauss) = (/ -.861136311594053D0, -.339981043584856D0, .339981043584856D0, .861136311594053D0 /)
36  real*8,parameter :: Wgauss(Ngauss) = (/ .347854845137454D0, .652145154862546D0, .652145154862546D0, .347854845137454D0 /)
37 endmodule
38 
39 module gts_values
40  implicit none
41  real*8, allocatable :: values(:), deltas(:) ! the state vector
42 endmodule
43 
44 module gts_schema ! parameters for the time integration scheme
45  implicit none
46  real*8 :: theta, zeta, gamma, hyper
47 endmodule
48 
49 module gts_coefficients ! the coefficients of the equations
50  implicit none
51  real*8, allocatable :: mass(:,:,:), diff(:,:,:), conv(:,:,:), reac(:,:,:), source(:,:), hyper2(:,:,:)
52 endmodule
53 
54 
55 subroutine solution6(solver,ifail)
56 !-----------------------------------------------------------------------
57 ! main GTS routine advancing the variables by n_steps time steps
58 ! - using single step Gears scheme (requires values from previous time step)
59 ! - or two step step TRBDF2 scheme
60 !-----------------------------------------------------------------------
61 use type_solver
62 
63 use gts_parameters
64 use gts_gauss
65 use gts_grid
66 use gts_matrix
67 use gts_values
68 use gts_schema
69 use gts_coefficients
70 
71 implicit none
72 
73 type (numerics) :: solver
74 integer :: ifail, n_steps, n_in, istep, i, j, ia
75 real*8 :: timestep, time, timestep1, timestep2
76 real*8 :: dxds, bnd_err1, bnd_err2
77 
78 !write(*,*) '*********************************'
79 !write(*,*) '* GTS solver *'
80 !write(*,*) '*********************************'
81 
82 neq = solver%ndim
83 n_in = solver%nrho
84 
85 time = 0.0d0 ! DPC addition
86 timestep = 1.d0
87 n_steps = 1
88 n_grid = n_in
89 hyper = 0.d-5 * solver%h
90 gamma = 0.d0
91 theta = 0.5d0
92 zeta = 0.d0
93 
94 if (.not. allocated(values)) allocate(values(n_grid * nvar * neq)) ! allocate state vector
95 
96 if (.not. allocated(deltas)) then
97  allocate(deltas(n_grid * nvar * neq)) ! allocate deltas previous step
98  deltas = 0.d0
99 endif
100 
101 !do i=1,solver%nrho
102 ! write(*,*) solver%y(1,i),solver%ym(1,i),solver%y(1,i)-solver%ym(1,i)
103 !enddo
104 
105 
106 call gts_construct_grid(solver%nrho,solver%rho) ! initialise the grid
107 
108 call gts_initialise_profiles(solver%nrho, solver%rho, solver%y, solver%ym) ! copy input to state vector
109 
110 call gts_initialise_coefficients(solver%nrho, solver%rho, solver%a, solver%b, solver%c, solver%d, &
111  solver%e, solver%f, solver%g, solver%h)
112 
113 
114 do istep=1, n_steps ! time stepping
115 
116  call gts_construct_matrix(solver%u, solver%v, solver%w, time, timestep) ! construct the FE matrix
117 
118  call gts_solve_matrix ! solve the system of equations
119 
120  values = values + bb ! add solution to state vector
121  deltas = bb ! keep change in values
122 
123  time = time + timestep
124 
125 enddo
126 
127 do i=1,n_grid
128  dxds = x_grid(2,i)
129  do j=1,neq ! copy state vector to output
130  ia = 2*neq*(i-1)+j
131  solver%ym(j,i) = solver%y(j,i) ! copy to previous timestep
132  solver%y(j,i) = values(ia)
133  solver%dy(j,i) = values(ia+neq) / dxds
134  enddo
135 enddo
136 
137 if (allocated(values)) deallocate(values) ! deallocate state vector
138 if (allocated(deltas)) deallocate(deltas) ! deallocate state vector
139 if (allocated(mass)) deallocate(mass)
140 if (allocated(diff)) deallocate(diff)
141 if (allocated(reac)) deallocate(reac)
142 if (allocated(conv)) deallocate(conv)
143 if (allocated(source)) deallocate(source)
144 if (allocated(hyper2)) deallocate(hyper2)
145 if (allocated(aa)) deallocate(aa)
146 if (allocated(bb)) deallocate(bb)
147 if (allocated(x_grid)) deallocate(x_grid)
148 if (allocated(xg)) deallocate(xg)
149 
150 return
151 
152 contains
153 
154 !************************************************************************
155  REAL*8 FUNCTION spwert(N,XWERT,A,B,C,D,X,ABLTG)
156 !-----------------------------------------------------------------------
157 ! INPUT:
158 !
159 ! N NUMBER OF GRID POINTS
160 ! XWERT STELLE AN DER FUNKTIONSWERTE BERECHNET WERDEN
161 ! A, B, C, D ARRAYS DER SPLINEKOEFFIZIENTEN (AUS SPLINE)
162 ! X ARRAY DER KNOTENPUNKTE
163 !
164 ! OUTPUT:
165 !
166 ! SPWERT FUNKTIONSWERT AN DER STELLE XWERT
167 ! ABLTG(I) I=1 : FIRST DERIVATIVE, ETC.
168 !-----------------------------------------------------------------------
169 
170  IMPLICIT NONE
171  INTEGER n
172  REAL*8 xwert, a(n), b(n), c(n), d(n), x(n), abltg(3), xx
173  INTEGER i, k, m
174 
175 ! SUCHE PASSENDES INTERVALL (BINAERE SUCHE)
176 
177  i = 1
178  k = n
179 
180 10 m = (i+k) / 2
181 
182  IF(m.NE.i) THEN
183  IF(xwert.LT.x(m)) THEN
184  k = m
185  ELSE
186  i = m
187  ENDIF
188  goto 10
189  ENDIF
190 
191  xx = xwert - x(i)
192 
193  abltg(1) = (3.0d0 * d(i) * xx + 2.0d0 * c(i)) * xx + b(i)
194  abltg(2) = 6.0d0 * d(i) * xx + 2.0d0 * c(i)
195  abltg(3) = 6.0d0 * d(i)
196 
197  spwert = ((d(i)*xx + c(i))*xx + b(i))*xx + a(i)
198 
199  RETURN
200  END FUNCTION spwert
201 
202 subroutine gts_initialise_coefficients(n_in, x_in, gts_a, gts_b, gts_c, gts_d, gts_e, gts_f, gts_g, gts_h)
203 !-----------------------------------------------------------------------------
204 ! subroutine to initialise the coeffients of the equations on the
205 ! Gaussian points
206 !-----------------------------------------------------------------------------
207 use gts_parameters
208 use gts_grid
209 use gts_coefficients
210 use gts_gauss
211 use gts_schema
212 implicit none
213 real*8 :: gts_a(neq,n_in), gts_b(neq,n_in), gts_c(neq,n_in), gts_d(neq,n_in)
214 real*8 :: gts_e(neq,n_in), gts_f(neq,n_in), gts_g(neq,n_in), gts_h
215 real*8 :: x_in(n_in), c_in(n_in), spline_a(n_in), spline_b(n_in), spline_c(n_in), spline_d(n_in)
216 real*8 :: derivs(3), scale
217 integer :: n_in, index, in, j, m, npoints, i
218 !!!real*8 :: spwert
219 
220 npoints = ngauss * (n_grid-1)
221 
222 if (.not. allocated(mass)) allocate(mass(npoints,neq,neq))
223 if (.not. allocated(diff)) allocate(diff(npoints,neq,neq))
224 if (.not. allocated(reac)) allocate(reac(npoints,neq,neq))
225 if (.not. allocated(conv)) allocate(conv(npoints,neq,neq))
226 if (.not. allocated(source)) allocate(source(npoints,neq))
227 if (.not. allocated(hyper2)) allocate(hyper2(npoints,neq,neq))
228 
229 mass = 0.d0
230 diff = 0.d0
231 conv = 0.d0
232 reac = 0.d0
233 source = 0.d0
234 hyper2 = 0.d0
235 
236 do in=1,neq
237 
238  c_in(1:n_in) = gts_b(in,1:n_in) * gts_c(in,1:n_in) ! mass coefficient
239 
240  scale = maxval(c_in)
241 
242  c_in = c_in / scale
243 
244  call spline(n_in,x_in,c_in,0.0d0,0.0d0,2,spline_a,spline_b,spline_c,spline_d)
245 
246  do j=1, n_grid-1
247  do m=1, ngauss
248  index = ngauss*(j-1) + m
249  mass(index,in,in) = spwert(n_in,xg(index),spline_a,spline_b,spline_c,spline_d,x_in,derivs)
250  enddo
251  enddo
252 
253  c_in(1:n_in) = gts_d(in,1:n_in) * gts_h ! diffusion coefficient
254 
255  c_in = c_in / scale
256 
257  call spline(n_in,x_in,c_in,0.0d0,0.0d0,2,spline_a,spline_b,spline_c,spline_d)
258 
259  do j=1, n_grid-1
260  do m=1, ngauss
261  index = ngauss*(j-1) + m
262  diff(index,in,in) = spwert(n_in,xg(index),spline_a,spline_b,spline_c,spline_d,x_in,derivs)
263  hyper2(index,in,in) = 0.001d0 * gts_h * derivs(2)**2
264  enddo
265  enddo
266 
267  c_in(1:n_in) = gts_e(in,1:n_in) * gts_h ! convection coefficient
268 
269  c_in = c_in / scale
270 
271  call spline(n_in,x_in,c_in,0.0d0,0.0d0,2,spline_a,spline_b,spline_c,spline_d)
272 
273  do j=1, n_grid-1
274  do m=1, ngauss
275  index = ngauss*(j-1) + m
276  conv(index,in,in) = spwert(n_in,xg(index),spline_a,spline_b,spline_c,spline_d,x_in,derivs)
277  enddo
278  enddo
279 
280  c_in(1:n_in) = gts_c(in,1:n_in) * (gts_b(in,1:n_in) - gts_a(in,1:n_in) - gts_g(in,1:n_in) * gts_h ) ! reaction coefficient
281 
282  c_in = c_in / scale
283 
284  call spline(n_in,x_in,c_in,0.0d0,0.0d0,2,spline_a,spline_b,spline_c,spline_d)
285 
286  do j=1, n_grid-1
287  do m=1, ngauss
288  index = ngauss*(j-1) + m
289  reac(index,in,in) = spwert(n_in,xg(index),spline_a,spline_b,spline_c,spline_d,x_in,derivs)
290  enddo
291  enddo
292 
293 
294  c_in(1:n_in) = gts_c(in,1:n_in) * gts_f(in,1:n_in) * gts_h ! source coefficient
295 
296  c_in = c_in / scale
297 
298  call spline(n_in,x_in,c_in,0.0d0,0.0d0,2,spline_a,spline_b,spline_c,spline_d)
299 
300  do j=1, n_grid-1
301  do m=1, ngauss
302  index = ngauss*(j-1) + m
303  source(index, in) = spwert(n_in,xg(index),spline_a,spline_b,spline_c,spline_d,x_in,derivs)
304  enddo
305  enddo
306 
307 enddo
308 
309 return
310 end subroutine gts_initialise_coefficients
311 
312 subroutine gts_initialise_profiles(n_in, x_in, y_in, ym_in)
313 !------------------------------------------------------------------------------
314 ! subroutine to initialise the state vector
315 !------------------------------------------------------------------------------
316 use gts_parameters
317 use gts_grid
318 use gts_values
319 implicit none
320 integer :: n_in, i, j, ia
321 real*8 :: x_in(n_in), y_in(neq,n_in), ym_in(neq,n_in), delta(n_in)
322 real*8 :: y_a(n_in), y_b(n_in), y_c(n_in), y_d(n_in)
323 real*8 :: d_a(n_in), d_b(n_in), d_c(n_in), d_d(n_in)
324 real*8 :: derivs(3), dxds
325 !!!real*8 :: spwert
326 
327 !------------------------------- convert to finite element grid and fill in derivatives
328 do j=1,neq
329 
330  call spline(n_in,x_in,y_in(j,1:n_in),0.0d0,0.0d0,2,y_a,y_b,y_c,y_d)
331 
332  delta = y_in(j,:) - ym_in(j,:)
333  call spline(n_in,x_in,delta,0.0d0,0.0d0,2,d_a,d_b,d_c,d_d)
334 
335  do i=1,n_grid
336 
337  dxds = x_grid(2,i)
338 
339  ia = 2*neq*(i-1)+j
340 
341  values(ia) = spwert(n_in,x_grid(1,i),y_a,y_b,y_c,y_d,x_in,derivs)
342  values(ia+neq) = derivs(1) * dxds
343 
344  deltas(ia) = spwert(n_in,x_grid(1,i),d_a,d_b,d_c,d_d,x_in,derivs)
345  deltas(ia+neq) = derivs(1) * dxds
346 
347  enddo
348 
349 enddo
350 
351 return
352 end subroutine gts_initialise_profiles
353 
354 subroutine gts_element_matrix(ife, time, timestep, ELM, RHS )
355 !----------------------------------------------------------
356 ! subroutine calculates the element matrix for one element
357 ! ife : the element number to be calculated (input)
358 ! norder : the order of the finite elements (input)
359 ! ELM : the element matrix (output)
360 ! RHS : the right hand side (output)
361 !----------------------------------------------------------
362 use gts_parameters
363 use gts_gauss
364 use gts_grid
365 use gts_values
366 use gts_coefficients
367 use gts_schema
368 implicit none
369 real*8 :: elm(2*nvar*neq,2*nvar*neq),rhs(2*nvar*neq)
370 real*8 :: fe(2*nvar), dfe(2*nvar), ddfe(2*nvar)
371 real*8 :: time, timestep, eq0(neq), deq0(neq), ddeq0(neq), delta_eq0(neq)
372 real*8 :: xjac, dxds, xright, xleft
373 real*8 :: sm, wm, xm
374 real*8 :: hyper_n(neq,neq)
375 real*8 :: v, dv, ddv
376 integer :: mgauss, ia, ja, i, j, ife, m, k, l, in, jn, index
377 
378 elm = 0.d0
379 rhs = 0.d0
380 
381 hyper_n = 0.d0
382 do in=1,neq
383  hyper_n(in,in) = hyper
384 enddo
385 
386 mgauss = ngauss ! the order of the Gaussian integration
387 
388 do m=1,mgauss ! loop over order of integration
389 
390  sm = xgauss(m) ! the gaussian integration points ( -1.< sm <1. )
391  wm = wgauss(m) ! the weights
392 
393  call gts_finite_elements(sm,fe,dfe,ddfe) ! the finite elements and its derivative
394 
395  index = (ife-1)*ngauss + m
396 
397  xm = x_grid(1,ife) * fe(1) + x_grid(2,ife) * fe(2) & ! the value of the coordinate x at the gaussian point
398  + x_grid(1,ife+1) * fe(3) + x_grid(2,ife+1) * fe(4)
399  dxds = x_grid(1,ife) * dfe(1) + x_grid(2,ife) * dfe(2) & ! scaling factor for derivatives in the master element
400  + x_grid(1,ife+1) * dfe(3) + x_grid(2,ife+1) * dfe(4)
401 
402  xjac = dxds ! the jacobian
403 
404 
405  eq0 = 0.d0 ! the value in between the nodes is the sum of contributions from two nodes
406  deq0 = 0.d0 ! the derivative
407  ddeq0 = 0.d0 ! second derivative
408  delta_eq0 =0.d0 ! value of delta previous timestep
409 
410  do in=1, neq ! reconstructed values of all variables
411  do k=1,2*nvar
412  ia = (ife-1) * nvar * neq + (k-1)*neq + in
413  eq0(in) = eq0(in) + values(ia) * fe(k)
414  deq0(in) = deq0(in) + values(ia) * dfe(k) / dxds
415  ddeq0(in) = ddeq0(in) + values(ia) * ddfe(k) / dxds**2 ! careful only for linear X(S)
416  delta_eq0(in) = delta_eq0(in) + deltas(ia) * fe(k)
417  enddo
418  enddo
419 
420  do in=1,neq
421 
422  do i=1,2*nvar
423 
424  ia = (i-1)*neq + in
425 
426  v = fe(i)
427  dv = dfe(i) / dxds
428  ddv = ddfe(i) / dxds**2 ! careful only for linear X(S)
429 
430  do jn=1,neq
431 
432  do j=1,2*nvar
433 
434  ja = (j-1)*neq + jn
435 
436  elm(ia,ja) = elm(ia,ja) + wm * xjac * ( mass(index,in,jn) * v * fe(j) * (1.d0 + zeta) &
437  + diff(index,in,jn) * theta * timestep * dv * dfe(j) / dxds &
438  - conv(index,in,jn) * theta * timestep * dv * fe(j) &
439  - reac(index,in,jn) * theta * timestep * v * fe(j) ) &
440 
441  + wm * xm * xjac * theta * timestep * hyper_n(in,jn) * ddv * ddfe(j)/dxds**2
442 
443 ! + wm * xm * Xjac * theta * timestep * hyper2(index,in,jn) * DDV * DDFE(j)/dXds**2
444  enddo
445 
446  enddo
447 
448  rhs(ia) = rhs(ia) + wm * xjac * timestep * ( &
449  - diff(index,in,in) * deq0(in) * dv &
450  + conv(index,in,in) * eq0(in) * dv &
451  + reac(index,in,in) * eq0(in) * v &
452  + source(index,in) * v ) &
453  - wm * xm * xjac * timestep * hyper_n(in,in) * ddeq0(in) * ddv &
454 
455 ! - wm * xm * Xjac * timestep * hyper2(index,in,in) * ddeq0(in) * DDV &
456 
457  + wm * xjac * v * zeta * mass(index,in,in) * delta_eq0(in)
458 
459  enddo
460 
461  enddo
462 
463 enddo
464 
465 return
466 end subroutine gts_element_matrix
467 
468 
469 subroutine gts_construct_grid(n_in,x_in)
470 !--------------------------------------------------
471 ! subroutine to define the grid
472 !--------------------------------------------------
473 use gts_parameters
474 use gts_grid
475 use gts_gauss
476 
477 implicit none
478 real*8 :: x_in(n_in)
479 integer :: n_in, i, m, index,k, l
480 real*8 :: s_in(n_in), spline_a(n_in), spline_b(n_in), spline_c(n_in), spline_d(n_in)
481 real*8 :: fe(2*nvar), dfe(2*nvar), ddfe(2*nvar)
482 real*8 :: s_out, derivs(3),x_tmp,dxds
483 !!! real*8 :: spwert
484 
485 !n_grid = n_in
486 
487 if (.not. allocated(x_grid)) allocate(x_grid(2,n_grid))
488 if (.not. allocated(xg)) allocate(xg(ngauss*(n_grid-1)))
489 
490 do i=1,n_in ! the global equidistant coordinate
491  s_in(i) = float(i-1)/float(n_in-1)
492 enddo
493 
494 call spline(n_in,s_in,x_in,0.0d0,0.0d0,2,spline_a,spline_b,spline_c,spline_d)
495 
496 x_grid(1,1:n_grid) = x_in(1:n_in)
497 
498 do i=1, n_grid
499  s_out = float(i-1)/float(n_grid-1)
500  x_tmp = spwert(n_in,s_out,spline_a,spline_b,spline_c,spline_d,s_in,derivs)
501  x_grid(2,i) = derivs(1) / (2.d0 * float(n_grid-1))
502 enddo
503 
504 xg = 0.d0 ! grid at the gaussian points
505 
506 do i=1, n_grid-1 ! definition of Gaussian points
507 
508  do m=1,ngauss
509 
510  call gts_finite_elements(xgauss(m),fe,dfe,ddfe) ! the finite elements and its derivative
511 
512  index = ngauss*(i-1) + m
513 
514  xg(index) = x_grid(1,i) * fe(1) + x_grid(2,i) * fe(2) + x_grid(1,i+1) * fe(3) + x_grid(2,i+1) * fe(4)
515 
516  enddo
517 
518 enddo
519 
520 return
521 end subroutine gts_construct_grid
522 
523 
524 subroutine gts_construct_matrix(u,v,w,time,timestep)
525 !------------------------------------------------------------------------
526 ! subroutine to assemble the fem matrices and apply boundary conditions
527 ! u(i,1)*dy(i,1) + v(i,1)*y(i,1) = w(i,1) (i=1,neq)
528 ! u(i,2)*dy(i,2) + v(i,2)*y(i,2) = w(i,2) (i=1,neq)
529 !------------------------------------------------------------------------
530 use gts_parameters
531 use gts_gauss
532 use gts_grid
533 use gts_matrix
534 use gts_values
535 implicit none
536 real*8 :: elm(2*nvar*neq,2*nvar*neq),rhs(2*nvar*neq) ! the element matrix and right hand side
537 real*8 :: u(neq,2), v(neq,2), w(neq,2), time, timestep, zbig
538 integer :: noff, nfe, ife, i, ia, j, ja, in
539 
540 
541 ku = 1 + 2*(neq*nvar-1) ! the number of diagonals above the main diagonal
542 kl = ku ! the number of diagonals below the main diagonal
543 lda = 2*kl+ku+1 ! the leading dimension
544 noff = kl+ku+1 ! the offset for band storage
545 
546 n = n_grid * nvar * neq ! the number of variables
547 
548 if (.not. allocated(aa)) ALLOCATE(aa(lda,n)) ! allocate a band matrix
549 if (.not. allocated(bb)) ALLOCATE(bb(n)) ! allocate the righthandside vector
550 
551 aa = 0.d0 ! initialise to zero
552 bb = 0.d0
553 
554 nfe = n_grid - 1 ! the number of finite elements
555 
556 do ife=1,nfe ! collect contributions from all elements
557 
558  call gts_element_matrix(ife,time,timestep,elm,rhs)! calculate the element matrix
559 
560  do i=1,2*nvar*neq ! add element matrix to main matrix
561 
562  ia = (ife-1) * nvar * neq + i
563 
564  do j=1,2*nvar*neq
565 
566  ja = (ife-1) * nvar * neq + j
567 
568  aa(noff+ia-ja,ja) = aa(noff+ia-ja,ja) + elm(i,j)
569 
570  enddo
571 
572  bb(ia) = bb(ia) + rhs(i)
573 
574  enddo
575 
576 enddo ! end of loop over elements
577 
578 !---------------------- boundary conditions, penalty 'method' : v*dy + u*y = w
579 zbig=1.d20
580 do in=1,neq
581  ia = neq + in
582  ja = neq + in
583  aa(noff+ia-ja,ja) = v(in,1) * zbig / x_grid(2,1)
584  bb(ia) = w(in,1) * zbig - v(in,1) * values(ja) * zbig / x_grid(2,1)
585 
586  ja = in
587  aa(noff+ia-ja,ja) = u(in,1) * zbig
588  bb(ia) = bb(ia) - u(in,1) * values(ja) * zbig
589 
590 enddo
591 
592 do in=1,neq
593  ia = (n_grid-1) * nvar * neq + in
594  ja = (n_grid-1) * nvar * neq + in
595 
596  aa(noff+ia-ja,ja) = u(in,2) * zbig
597  bb(ia) = w(in,2) * zbig - u(in,2) * values(ja) * zbig
598 
599  ja = (n_grid-1) * nvar * neq + in + neq
600  aa(noff+ia-ja,ja) = v(in,2) * zbig / x_grid(2,n_grid)
601  bb(ia) = bb(ia) - v(in,2) * values(ja) * zbig / x_grid(2,n_grid)
602 
603 enddo
604 
605 return
606 end subroutine gts_construct_matrix
607 
608 
609 subroutine gts_finite_elements(xm,FE,DFE,DDFE)
610 !--------------------------------------------------
611 ! the finite elements in the master element
612 ! defined between -1 and +1
613 !--------------------------------------------------
614 implicit none
615 real*8 :: fe(*),dfe(*),ddfe(*), xm
616 
617  fe(1) = -0.25d0*(xm - 1.d0)**2 * (-xm - 2.d0); dfe(1) = -0.5d0*(xm - 1.d0) * (-xm - 2.d0) + 0.25d0*(xm - 1.d0)**2
618  fe(3) = -0.25d0*(xm + 1.d0)**2 * ( xm - 2.d0); dfe(3) = -0.5d0*(xm + 1.d0) * ( xm - 2.d0) - 0.25d0*(xm + 1.d0)**2
619  fe(2) = -0.25d0*(xm - 1.d0)**2 * (-xm - 1.d0); dfe(2) = -0.5d0*(xm - 1.d0) * (-xm - 1.d0) + 0.25d0*(xm - 1.d0)**2
620  fe(4) = +0.25d0*(xm + 1.d0)**2 * ( xm - 1.d0); dfe(4) = +0.5d0*(xm + 1.d0) * ( xm - 1.d0) + 0.25d0*(xm + 1.d0)**2
621 
622  ddfe(1) = -0.5d0*(-xm - 2.d0) + 0.5d0*(xm - 1.d0) + 0.5d0*(xm - 1.d0)
623  ddfe(3) = -0.5d0*( xm - 2.d0) - 0.5d0*(xm + 1.d0) - 0.5d0*(xm + 1.d0)
624  ddfe(2) = -0.5d0*(-xm - 1.d0) + 0.5d0*(xm - 1.d0) + 0.5d0*(xm - 1.d0)
625  ddfe(4) = +0.5d0*( xm - 1.d0) + 0.5d0*(xm + 1.d0) + 0.5d0*(xm + 1.d0)
626 
627 return
628 end subroutine gts_finite_elements
629 
630 subroutine gts_solve_matrix
631 !--------------------------------------------------
632 ! Solves the banded system of equations
633 ! (for general up-down asymmetric matrices)
634 ! using Lapack routines dgbtrf/dgbtrs
635 !--------------------------------------------------
636 use gts_matrix
637 implicit none
638 integer, allocatable :: ipiv(:),iwork(:)
639 real*8, allocatable :: work(:)
640 real*8 :: info
641 
642 allocate(ipiv(n), work(3*n), iwork(n))
643 
644 call dgbtrf(n,n,kl,ku,aa,lda,ipiv,info) ! lapack routine for LU decomposition
645 
646 call dgbtrs('N',n,kl,ku,1,aa,lda,ipiv,bb,n,info) ! lapack solve equations (back substitution)
647 
648 deallocate(ipiv, work, iwork)
649 
650 return
651 end subroutine gts_solve_matrix
652 
653 
654 !***********************************************************************
655  SUBROUTINE spline(N,X,Y,ALFA,BETA,TYP,A,B,C,D)
656 !-----------------------------------------------------------------------
657 ! INPUT:
658 !
659 ! N NUMBER OF POINTS
660 ! X ARRAY X VECTOR
661 ! Y ARRAY Y VECTOR
662 ! ALFA BOUNDARY CONDITION IN X(1)
663 ! BETA " IN X(N)
664 ! TYP = 0 NOT-A-KNOT SPLINE
665 ! 1 ALFA, BETA 1. ABLEITUNGEN VORGEGEBEN
666 ! 2 " " 2. " "
667 ! 3 " " 3. " "
668 !
669 ! BEMERKUNG: MIT TYP = 2 UND ALFA = BETA = 0 ERHAELT MAN
670 ! EINEN NATUERLICHEN SPLINE
671 !
672 ! OUTPUT:
673 !
674 ! A, B, C, D ARRAYS OF SPLINE COEFFICIENTS
675 ! S = A(I) + B(I)*(X-X(I)) + C(I)*(X-X(I))**2+ D(I)*(X-X(I))**3
676 !
677 ! BEI ANWENDUNGSFEHLERN WIRD DAS PROGRAMM MIT ENTSPRECHENDER
678 ! FEHLERMELDUNG ABGEBROCHEN
679 !-----------------------------------------------------------------------
680 !
681 !
682  IMPLICIT NONE
683 
684  INTEGER n, typ
685  REAL*8 x(n), y(n), alfa, beta, a(n), b(n), c(n), d(n)
686  INTEGER i, ierr
687  REAL*8 h(n)
688 
689  IF((typ.LT.0).OR.(typ.GT.3)) THEN
690  WRITE(*,*) 'ERROR IN ROUTINE SPLINE: FALSE TYP'
691  stop
692  ENDIF
693 
694  IF (n.LT.3) THEN
695  WRITE(*,*) 'ERROR IN ROUTINE SPLINE: N < 3'
696  stop
697  ENDIF
698 
699 
700 ! BERECHNE DIFFERENZ AUFEINENDERFOLGENDER X-WERTE UND
701 ! UNTERSUCHE MONOTONIE
702 !
703  DO i = 1, n-1
704  h(i) = x(i+1)- x(i)
705  IF(h(i).LE.0.0d0) THEN
706  WRITE(*,*) 'NON MONOTONIC ABCISSAE IN SPLINE: X(I-1)>=X(I)'
707  stop
708  ENDIF
709  ENDDO
710 !
711 ! AUFSTELLEN DES GLEICHUNGSSYSTEMS
712 !
713  DO 20 i = 1, n-2
714  a(i) = 3.0d0 * ((y(i+2)-y(i+1)) / h(i+1) - (y(i+1)-y(i)) / h(i))
715  b(i) = h(i)
716  c(i) = h(i+1)
717  d(i) = 2.0d0 * (h(i) + h(i+1))
718  20 CONTINUE
719 !
720 ! BERUECKSICHTIGEN DER RANDBEDINGUNGEN
721 
722 ! NOT-A-KNOT
723 
724  IF(typ.EQ.0) THEN
725  a(1) = a(1) * h(2) / (h(1) + h(2))
726  a(n-2) = a(n-2) * h(n-2) / (h(n-1) + h(n-2))
727  d(1) = d(1) - h(1)
728  d(n-2) = d(n-2) - h(n-1)
729  c(1) = c(1) - h(1)
730  b(n-2) = b(n-2) - h(n-1)
731  ENDIF
732 
733 ! 1. ABLEITUNG VORGEGEBEN
734 
735  IF(typ.EQ.1) THEN
736  a(1) = a(1) - 1.5d0 * ((y(2)-y(1)) / h(1) - alfa)
737  a(n-2) = a(n-2) - 1.5d0 * (beta - (y(n)-y(n-1)) / h(n-1))
738  d(1) = d(1) - 0.5d0 * h(1)
739  d(n-2) = d(n-2) - 0.5d0 * h(n-1)
740  ENDIF
741 !
742 ! 2. ABLEITUNG VORGEGEBEN
743 !
744  IF(typ.EQ.2) THEN
745  a(1) = a(1) - 0.5d0 * alfa * h(1)
746  a(n-2) = a(n-2) - 0.5d0 * beta * h(n-1)
747  ENDIF
748 
749 ! 3. ABLEITUNG VORGEGEBEN
750 !
751  IF(typ.EQ.3 ) THEN
752  a(1) = a(1) + 0.5d0 * alfa * h(1) * h(1)
753  a(n-2) = a(n-2) - 0.5d0 * beta * h(n-1)* h(n-1)
754  d(1) = d(1) + h(1)
755  d(n-2) = d(n-2) + h(n-1)
756  ENDIF
757 
758 ! BERECHNUNG DER KOEFFIZIENTEN
759 !
760  CALL dgtsl(n-2,b,d,c,a,ierr)
761  IF(ierr.NE.0) THEN
762  WRITE(*,21)
763  stop
764  ENDIF
765 
766 ! UEBERSCHREIBEN DES LOESUNGSVEKTORS
767 
768  CALL dcopy(n-2,a,1,c(2),1)
769 !
770 ! IN ABHAENGIGKEIT VON DEN RANDBEDINGUNGEN WIRD DER 1. UND
771 ! DER LETZTE WERT VON C KORRIGIERT
772 !
773  IF(typ.EQ.0) THEN
774  c(1) = c(2) + h(1) * (c(2)-c(3)) / h(2)
775  c(n) = c(n-1) + h(n-1) * (c(n-1)-c(n-2)) / h(n-2)
776  ENDIF
777 
778  IF(typ.EQ.1) THEN
779  c(1) = 1.5d0*((y(2)-y(1)) / h(1) - alfa) / h(1) - 0.5d0 * c(2)
780  c(n) = -1.5d0*((y(n)-y(n-1)) / h(n-1)- beta) / h(n-1)- 0.5d0 * c(n-1)
781  ENDIF
782 
783  IF(typ.EQ.2) THEN
784  c(1) = 0.5d0 * alfa
785  c(n) = 0.5d0 * beta
786  ENDIF
787 
788  IF(typ.EQ.3) THEN
789  c(1) = c(2) - 0.5d0 * alfa * h(1)
790  c(n) = c(n-1) + 0.5d0 * beta * h(n-1)
791  ENDIF
792 
793  CALL dcopy(n,y,1,a,1)
794 
795  DO i = 1, n-1
796  b(i) = (a(i+1)-a(i)) / h(i) - h(i) * (c(i+1)+2.0d0 * c(i)) / 3.0d0
797  d(i) = (c(i+1)-c(i)) / (3.0d0 * h(i))
798  END DO
799 
800  b(n) = (3.0d0 * d(n-1) * h(n-1) + 2.0d0 * c(n-1)) * h(n-1) + b(n-1)
801 
802  RETURN
803 
804  21 FORMAT(1x,'ERROR IN SGTSL: MATRIX SINGULAR')
805  END SUBROUTINE spline
806 
807 
808 !DECK DGTSL CAS02750
809 !** FROM NETLIB, TUE AUG 28 08:28:34 EDT 1990 ***
810 !** COPIED FROM SGTSL AND RENAMED ***
811 !
812  SUBROUTINE dgtsl(N,C,D,E,B,INFO)
813  IMPLICIT NONE
814  INTEGER n,info
815  REAL*8 c(*),d(*),e(*),b(*)
816 
817 ! SGTSL GIVEN A GENERAL TRIDIAGONAL MATRIX AND A RIGHT HAND
818 ! SIDE WILL FIND THE SOLUTION.
819 !
820 ! ON ENTRY
821 !
822 ! N INTEGER
823 ! IS THE ORDER OF THE TRIDIAGONAL MATRIX.
824 !
825 !
826 ! C REAL(N)
827 ! IS THE SUBDIAGONAL OF THE TRIDIAGONAL MATRIX.
828 ! C(2) THROUGH C(N) SHOULD CONTAIN THE SUBDIAGONAL.
829 ! ON OUTPUT C IS DESTROYED.
830 !
831 ! D REAL(N)
832 ! IS THE DIAGONAL OF THE TRIDIAGONAL MATRIX.
833 ! ON OUTPUT D IS DESTROYED.
834 !
835 ! E REAL(N)
836 ! IS THE SUPERDIAGONAL OF THE TRIDIAGONAL MATRIX.
837 ! E(1) THROUGH E(N-1) SHOULD CONTAIN THE SUPERDIAGONAL.
838 ! ON OUTPUT E IS DESTROYED.
839 !
840 ! B REAL(N)
841 ! IS THE RIGHT HAND SIDE VECTOR.
842 !
843 ! ON RETURN
844 !
845 ! B IS THE SOLUTION VECTOR.
846 !
847 ! INFO INTEGER
848 ! = 0 NORMAL VALUE.
849 ! = K IF THE K-TH ELEMENT OF THE DIAGONAL BECOMES
850 ! EXACTLY ZERO. THE SUBROUTINE RETURNS WHEN
851 ! THIS IS DETECTED.
852 !
853 ! LINPACK. THIS VERSION DATED 08/14/78 .
854 ! JACK DONGARRA, ARGONNE NATIONAL LABORATORY.
855 !
856 ! NO EXTERNALS
857 ! FORTRAN ABS
858 !
859 ! INTERNAL VARIABLES
860 !
861  INTEGER k,kb,kp1,nm1,nm2
862  REAL*8 t
863 ! BEGIN BLOCK PERMITTING ...EXITS TO 100
864 
865  info = 0
866  c(1) = d(1)
867  nm1 = n - 1
868  IF (nm1 .LT. 1) go to 40
869  d(1) = e(1)
870  e(1) = 0.0d0
871  e(n) = 0.0d0
872 
873  DO 30 k = 1, nm1
874  kp1 = k + 1
875 
876 ! FIND THE LARGEST OF THE TWO ROWS
877 
878  IF (abs(c(kp1)) .LT. abs(c(k))) go to 10
879 
880 ! INTERCHANGE ROW
881 
882  t = c(kp1)
883  c(kp1) = c(k)
884  c(k) = t
885  t = d(kp1)
886  d(kp1) = d(k)
887  d(k) = t
888  t = e(kp1)
889  e(kp1) = e(k)
890  e(k) = t
891  t = b(kp1)
892  b(kp1) = b(k)
893  b(k) = t
894  10 CONTINUE
895 
896 ! ZERO ELEMENTS
897 
898  IF (c(k) .NE. 0.0d0) go to 20
899  info = k
900 ! ............EXIT
901  go to 100
902  20 CONTINUE
903  t = -c(kp1)/c(k)
904  c(kp1) = d(kp1) + t*d(k)
905  d(kp1) = e(kp1) + t*e(k)
906  e(kp1) = 0.0d0
907  b(kp1) = b(kp1) + t*b(k)
908  30 CONTINUE
909  40 CONTINUE
910  IF (c(n) .NE. 0.0d0) go to 50
911  info = n
912  go to 90
913  50 CONTINUE
914 
915 ! BACK SOLVE
916 
917  nm2 = n - 2
918  b(n) = b(n)/c(n)
919  IF (n .EQ. 1) go to 80
920  b(nm1) = (b(nm1) - d(nm1)*b(n))/c(nm1)
921  IF (nm2 .LT. 1) go to 70
922  DO 60 kb = 1, nm2
923  k = nm2 - kb + 1
924  b(k) = (b(k) - d(k)*b(k+1) - e(k)*b(k+2))/c(k)
925  60 CONTINUE
926  70 CONTINUE
927  80 CONTINUE
928  90 CONTINUE
929  100 CONTINUE
930 
931  RETURN
932  END SUBROUTINE dgtsl
933 
934 end subroutine solution6
gts_matrix
Definition: solution6.f90:19
solution6
subroutine solution6(solver, ifail)
Definition: solution6.f90:55
solver
Definition: solver.f90:1
type_solver
Definition: type_solver.f90:1
gts_finite_elements
subroutine gts_finite_elements(xm, FE, DFE, DDFE)
Definition: solution6.f90:609
gts_coefficients
Definition: solution6.f90:49
gts_initialise_profiles
subroutine gts_initialise_profiles(n_in, x_in, y_in, ym_in)
Definition: solution6.f90:312
gts_parameters
Definition: solution6.f90:12
spline
subroutine spline(N, X, Y, ALFA, BETA, TYP, A, B, C, D)
Definition: solution6.f90:655
spwert
REAL *8 function spwert(N, XWERT, A, B, C, D, X, ABLTG)
Definition: solution6.f90:155
gts_solve_matrix
subroutine gts_solve_matrix
Definition: solution6.f90:630
gts_initialise_coefficients
subroutine gts_initialise_coefficients(n_in, x_in, gts_a, gts_b, gts_c, gts_d, gts_e, gts_f, gts_g, gts_h)
Definition: solution6.f90:202
gts_values
Definition: solution6.f90:39
gts_element_matrix
subroutine gts_element_matrix(ife, time, timestep, ELM, RHS)
Definition: solution6.f90:354
gts_gauss
Definition: solution6.f90:32
gts_construct_grid
subroutine gts_construct_grid(n_in, x_in)
Definition: solution6.f90:469
dgtsl
subroutine dgtsl(N, C, D, E, B, INFO)
Definition: solution6.f90:812
type_solver::numerics
Definition: type_solver.f90:22
gts_grid
Definition: solution6.f90:26
gts_schema
Definition: solution6.f90:44
gts_construct_matrix
subroutine gts_construct_matrix(u, v, w, time, timestep)
Definition: solution6.f90:524