ETS  \$Id: Doxyfile 2162 2020-02-26 14:16:09Z g2dpc $
 All Classes Files Functions Variables Pages
solution2.f90
Go to the documentation of this file.
1 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
11 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
12 
13 SUBROUTINE solution2 (SOLVER, ifail)
14 
15 !-------------------------------------------------------!
16 ! !
17 ! This subroutine is prepared to solve single !
18 ! transport equation in standardised form adopted !
19 ! by the ETS. !
20 ! !
21 !-------------------------------------------------------!
22 ! !
23 ! Source: provided by M.Tokar on 17.02.2009 !
24 ! Developers: M.Tokar, D.Kalupin !
25 ! Kontacts: M.Tokar@fz-juelich.de !
26 ! D.Kalupin@fz-juelich.de !
27 ! !
28 ! Comments: equation is solved in !
29 ! INTEGRAL form !
30 ! !
31 !-------------------------------------------------------!
32 
33 
34 
35 
36  USE type_solver
37 
38  USE itm_types
39 
40  IMPLICIT NONE
41 
42  INTEGER, INTENT (INOUT) :: ifail
43 
44 ! +++ Input/Output with ETS:
45  TYPE (numerics), INTENT (INOUT) :: solver !contains all I/O quantities to numerics part
46 
47 
48 ! +++ Internal input/output parameters:
49  INTEGER :: idim, ndim !equation index and total number of equations to be solved
50 
51  INTEGER :: irho, nrho !radius index, number of radial points
52 
53  INTEGER :: flag !flag for equation: 0 - interpretative (not solved), 1 - predictive (solved)
54 
55  REAL (R8) :: rho(solver%nrho) !radii
56 
57  REAL (R8) :: amix !fraction of new sollution mixed
58 
59  REAL (R8) :: y(solver%nrho), ym(solver%nrho) !function at the current amd previous time steps
60  REAL (R8) :: dy(solver%nrho) !derivative of function
61 
62  REAL (R8) :: a(solver%nrho), b(solver%nrho) !coefficients
63  REAL (R8) :: c(solver%nrho), d(solver%nrho) !coefficients
64  REAL (R8) :: e(solver%nrho), f(solver%nrho) !coefficients
65  REAL (R8) :: g(solver%nrho), h !coefficients
66  REAL (R8) :: dd(solver%nrho), de(solver%nrho) !derivatives of coefficients
67 
68  REAL (R8) :: v(2), u(2), w(2) !boundary conditions
69 
70 
71 
72 
73 
74 
75 ! + + + + + + + + + INTERFACE PART + + + + + + + + + + +
76 
77 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
78 ! +++ Set up local variables (input, obtained from ETS):
79 ! Control parameters:
80  ndim = solver%NDIM
81  nrho = solver%NRHO
82  amix = solver%AMIX
83 
84 
85 
86 ! +++ Solution of equations starting from 1 to NDIM:
87  equation_loop1: DO idim = 1, ndim
88 
89  flag = solver%EQ_FLAG(idim)
90 
91  IF (flag.EQ.0) goto 20 !equation is not solved
92 
93 
94 
95 ! +++ Numerical coefficients obtained from physics part in form:
96 !
97 ! (A*Y-B*Y(t-1))/H + 1/C * (-D*Y' + E*Y) = F - G*Y
98 
99  h = solver%H
100 
101  rho_loop1: DO irho=1,nrho
102  rho(irho) = solver%RHO(irho)
103 
104  y(irho) = solver%Y(idim,irho)
105  dy(irho) = solver%DY(idim,irho)
106  ym(irho) = solver%YM(idim,irho)
107 
108  a(irho) = solver%A(idim,irho)
109  b(irho) = solver%B(idim,irho)
110  c(irho) = solver%C(idim,irho)
111  d(irho) = solver%D(idim,irho)
112  e(irho) = solver%E(idim,irho)
113  f(irho) = solver%F(idim,irho)
114  g(irho) = solver%G(idim,irho)
115  END DO rho_loop1
116 
117 
118 ! +++ Boundary conditions for numerical solver in form:
119 !
120 ! V*Y' + U*Y = W
121 
122 ! +++ On axis:
123 
124  v(1) = solver%V(idim,1)
125  u(1) = solver%U(idim,1)
126  w(1) = solver%W(idim,1)
127 
128 ! +++ At the edge:
129 
130  v(2) = solver%V(idim,2)
131  u(2) = solver%U(idim,2)
132  w(2) = solver%W(idim,2)
133 
134 
135 
136 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
137 ! +++ Solution of transport equation:
138 
139  CALL solver_ritm_integral &
140  (nrho,rho,amix, y,ym,dy, a,b,c,d,e,f,g,h, u,v,w)
141 ! radii b.c. coeff. solution spline
142 
143 
144  rho_loop3: DO irho = 1,nrho
145 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
146 ! +++ Return solution to ETS:
147  solver%Y(idim,irho) = y(irho)
148  solver%DY(idim,irho) = dy(irho)
149 
150  END DO rho_loop3
151 
152 
153 
154 20 CONTINUE
155 
156 
157 
158  END DO equation_loop1
159 
160 
161 
162  RETURN
163 
164 
165 
166 END SUBROUTINE solution2
167 
168 
169 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
170 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
171 
172 
173 
174 
175 
176 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
177 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
178 
179 ! + + + + + + + + + NUMERICAL PART + + + + + + + + + + +
180  SUBROUTINE solver_ritm_integral (n,r,amix,z,zm,dz, &
181  aets,bets,cets,dets,eets,fets,gets,tau,uets,vets,wets)
182 
183  USE itm_types
184 
185  IMPLICIT NONE
186 
187  INTEGER i,n
188 
189  REAL (R8) tau,xi,p,amix,error,errornew
190 !DPC 2009-02-28 Changed h(n) to h(n+1) to match definition in solut3
191  REAL (R8) h(n+1),r(n),z(n),zm(n)
192  REAL (R8) dz(n),dzm(n)
193  REAL (R8) aets(n),bets(n),cets(n)
194  REAL (R8) dets(n),eets(n),gets(n),fets(n),gamets(n)
195  REAL (R8) uets(2),vets(2),wets(2)
196  REAL (R8) cs(n),csn(n)
197  REAL (R8) al(n),be(n),so(n),nu(n),dif(n),vel(n),dfl(n)
198  REAL (R8) om(n),dom(n),s(n),ints(n)
199  REAL (R8) a(n),b(n),f(n),u(2),v(2),w(2)
200  REAL (R8) y(n),d(n)
201  REAL (R8) dd(n)
202 
203 ! +++ ETS - standard equation in cylindrical geometry
204  DO i=2,n
205  al(i) = aets(i)*cets(i)/r(i)
206  be(i) = bets(i)*cets(i)/r(i)
207  so(i) = fets(i)*cets(i)/r(i)
208  nu(i) = gets(i)*cets(i)/r(i)
209  dif(i) = dets(i)/r(i)
210  vel(i) = eets(i)/r(i)
211  gamets(i) = 0.d0
212  dfl(i) = gamets(i)/r(i)
213  h(i-1) = r(i)-r(i-1)
214  END DO
215  h(n)=0.0_r8 !DPC 2009-02-28
216 
217  FORALL(i=1:n) cs(i)=0.d0
218 
219  CALL derivn2(n,r,zm,dzm)
220 
221 ! +++ parameters for integral representation
222  10 xi=1.d-1
223 
224  DO i=2,n
225  p = (1.d0-xi)*al(i)+tau*nu(i)
226  om(i) = xi*al(i)+(p+dabs(p))/2.d0
227  s(i) = so(i)+((p-dabs(p))/2.d0*z(i)-(om(i)-be(i))*zm(i))/tau
228  END DO
229 
230  CALL axis(n,h,r,om)
231 
232  CALL axis(n,h,r,s)
233 
234 
235 ! +++ coefficients for differential equation
236  CALL derivn2(n,r,om,dom)
237  CALL integral(n,h,r,s,ints)
238 
239  DO i=2,n
240  a(i) = 3.d0/r(i)-dom(i)/om(i)-vel(i)/dif(i)
241  b(i) = om(i)/tau/dif(i)+2.d0/r(i)*(3.d0/r(i)-a(i))
242  f(i) = om(i)/r(i)/dif(i)* &
243  (ints(i)/r(i)+dif(i)*dzm(i)-vel(i)*zm(i)-dfl(i))
244  END DO
245 
246 
247 ! +++ boundary conditions
248  IF(wets(1).NE.0.d0) THEN
249  write(6,*) 'wrong boundary condition: w_1 is not 0'
250  stop
251  END IF
252 
253  IF(uets(1).EQ.0.d0) THEN !dz/dr(0)=0
254  u(1)=b(2)
255  v(1)=-1.d0/r(2)-a(2)
256  w(1)=f(2)
257  END IF
258 
259  IF(vets(1).EQ.0.d0) THEN !z(0)=0
260  u(1) = b(2)+2.d0/r(2)**2
261  v(1) = -a(2)-2.d0/r(2)
262  w(1) = f(2)
263  END IF
264 
265  w(2) = vets(2)*f(n)+ &
266  om(n)/r(n)*(wets(2)-uets(2)*zm(n)-vets(2)*dzm(n))
267  v(2) = uets(2)+vets(2)*vel(n)/dif(n)
268  u(2) = vets(2)*om(n)/tau/dif(n)+2.d0/r(n)*v(2)
269 
270 
271 ! +++ solution: cs=z-zm
272  CALL solut3(n,h,a,b,f,u,v,w,y,d)
273 
274  CALL axis(n,h,r,y)
275 
276  IF(vets(1).EQ.0.d0) y(1)=0.d0
277 
278 ! CALL DERIVN2 (n,r,y,d)
279 
280  DO i=1,n
281  csn(i) = (r(i)*d(i)+2.d0*y(i))/om(i)
282  cs(i) = cs(i)*(1.d0-amix)+csn(i)*amix
283  z(i) = zm(i)+cs(i)
284  END DO
285 
286  CALL derivn2(n,r,z,dz)
287 
288 ! +++ error in solution
289  error=0.d0
290 
291  DO i=1,n
292  errornew =dsqrt(dabs(1.d0-csn(i)/cs(i)))
293  IF(cs(i).NE.0.d0.AND.error.LT.errornew) error=errornew
294  END DO
295 
296  IF (error.GT.1.d-5*amix) goto 10
297 
298 
299  RETURN
300 
301 
302  END SUBROUTINE solver_ritm_integral
303 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
304 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
305 
306 
307 
308 
309 
310 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
311 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
312  SUBROUTINE solut3 (n,h,a,b,f,u,v,w,y,dy)
313 !-------------------------------------------------------!
314 ! !
315 ! This subroutine provides solution for the equation!
316 ! !
317 ! d(dy/dx)/dx + a(x)*dy/dx = b(x)*y - f(x) !
318 ! !
319 !-------------------------------------------------------!
320  USE itm_types
321 
322  IMPLICIT NONE
323 
324  INTEGER i,j,k,l,m,n
325 
326  REAL (R8) h(n+1),a(n),b(n),f(n),y(n),dy(n),lam(10001,2)
327  REAL (R8) u(2),v(2),w(2)
328  REAL (R8) y0(n),y1(10001,2,2)
329  REAL (R8) et(2,2),mu(2),det
330  REAL (R8) c(n),d(n),g(2),x
331  REAL (R8) cdy(n,2)
332 
333  DO i=2,n
334  y0(i) = f(i)/b(i)
335  det = a(i)**2/4.d0+b(i)
336 
337  DO k=1,2
338  lam(i,k) = -a(i)/2.d0-(-1)**k*dsqrt(det)
339 
340 
341  FORALL(l=1:2) &
342  y1(i,k,l) = dexp(dmin1(200.d0,(-1)**l*lam(i,k) &
343  *h(i-2+l)/2.d0))
344  END DO
345 
346  END DO
347 
348 
349  c(3) = v(1)/(v(1)-h(3)*u(1))
350  d(3) = h(3)*w(1)/(h(3)*u(1)-v(1))
351 
352 
353  DO i=3,n-1
354 
355  DO m=1,2
356  l = 3-m
357  j = i+(-1)**m
358 
359  DO k=1,2
360  et(k,m) = -(-1)**m*(lam(j,m)-lam(i,k) &
361  -(lam(j,l)-lam(i,k))*y1(j,l,l) &
362  /y1(j,m,l))/(lam(j,1)-lam(j,2))
363  END DO
364 
365  mu(m) = -(-1)**m*(y0(i)-y0(j))/y1(i,l,m) &
366  *(lam(j,m)-lam(j,l)*y1(j,l,l) &
367  /y1(j,m,l))/(lam(j,1)-lam(j,2))
368  END DO
369 
370  det = et(1,2)*et(2,1)-et(1,1)*et(2,2) &
371  *y1(i,l,1)/y1(i,2,1) &
372  *y1(i,2,2)/y1(i,l,2)
373 
374  dy(i) = 0.0d0
375 
376  DO m=1,2
377  l = 3-m
378  g(m) = y1(i+(-1)**m,l,l)/y1(i,l,m)* &
379  (et(m,l)-y1(i,l,l)/y1(i,m,l) &
380  *et(l,l))/det
381 ! derivative:
382  dy(i) = dy(i)+lam(i,m)* &
383  (mu(m)*y1(i,l,l)/y1(i,m,l)* &
384  et(l,l)-mu(l)*et(l,m))/det
385  cdy(i,m) = (et(m,l)*lam(i,l)- &
386  y1(i,l,l)/y1(i,m,l)*et(l,l)*lam(i,m))/det
387  END DO
388 
389  x = (mu(1)*(y1(i,2,2)/y1(i,1,2) &
390  *et(2,2)-et(1,2))+mu(2)*(y1(i,1,1) &
391  /y1(i,2,1)*et(1,1)-et(2,1)))/det- &
392  y0(i-1)*g(1)-y0(i+1)*g(2)+y0(i)
393  c(i+1) = g(2)/(1.d0-g(1)*c(i))
394  d(i+1) = (g(1)*d(i)+x)/(1.d0-g(1)*c(i))
395 
396  END DO
397 
398  y(n) = (h(n-1)*w(2)+v(2)*d(n)) &
399  /(h(n-1)*u(2)+v(2)-v(2)*c(n))
400 
401  DO i=n,3,-1
402  y(i-1) = c(i)*y(i)+d(i)
403  END DO
404 
405 
406 ! derivative:
407  DO i=3,n-1
408  DO m=1,2
409  l = 3-m
410  j = i+(-1)**m
411  dy(i) = dy(i)+(y(j)-y0(j))*y1(j,l,l)/y1(i,l,m)*cdy(i,m)
412  END DO
413  END DO
414 
415  dy(2)=(y(3)-y(3))/h(2)
416  dy(n)=(y(n)-y(n-1))/h(n-1)
417 
418 
419  RETURN
420 
421 
422  END SUBROUTINE solut3
423 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
424 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
425 
426 
427 
428 
429 
430 
431 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
432 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
433  SUBROUTINE derivn2(N,X,Y,DY1)
434 !-------------------------------------------------------!
435 ! !
436 ! This subroutine calculates first and !
437 ! second derivatives, DY1 and DY2, !
438 ! of function Y respect to argument X !
439 ! !
440 !-------------------------------------------------------!
441  USE itm_types
442 
443  IMPLICIT NONE
444 
445  INTEGER n, i
446  REAL (R8) x(n), y(n), dy1(n)
447  REAL (R8) h(n), dy2(n)
448 
449  SAVE
450 
451 
452 
453 
454  DO i=2,n-1
455  dy1(i) = (y(i+1)-y(i-1))/(x(i+1)-x(i-1))
456  END DO
457 
458 
459  dy1(1) = (y(2)-y(1))/(x(2)-x(1))
460  dy1(n) = (y(n)-y(n-1))/(x(n)-x(n-1))
461 
462 
463  RETURN
464 
465 
466  END SUBROUTINE derivn2
467 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
468 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
469 
470 
471 
472 
473 
474 
475 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
476 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
477  SUBROUTINE derivn21(N,X,Y,DY1)
478 !-------------------------------------------------------!
479 ! !
480 ! This subroutine calculates first and !
481 ! second derivatives, DY1 and DY2, !
482 ! of function Y respect to argument X !
483 ! !
484 !-------------------------------------------------------!
485  USE itm_types
486 
487  IMPLICIT NONE
488 
489  INTEGER n, i
490  REAL (R8) x(n), y(n), dy1(n)
491  REAL (R8) h(n), dy2(n)
492 
493  SAVE
494 
495 
496  DO i=1,n-1
497  h(i) = x(i+1)-x(i)
498  END DO
499 
500 
501  DO i=2,n-1
502  dy1(i) = ((y(i+1)-y(i))*h(i-1)/h(i) &
503  +(y(i)-y(i-1))*h(i)/h(i-1))/(h(i)+h(i-1))
504  dy2(i) = 2.d0*((y(i-1)-y(i))/h(i-1) &
505  +(y(i+1)-y(i))/h(i))/(h(i)+h(i-1))
506  END DO
507 
508 
509  dy1(1) = dy1(2)-dy2(2)*h(1)
510  dy1(n) = dy1(n-1)+dy2(n-1)*h(n-1)
511 
512 
513  RETURN
514 
515 
516  END SUBROUTINE derivn21
517 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
518 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
519 
520 
521 
522 
523 
524 
525 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
526 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
527  SUBROUTINE integral (n, h, r, f, int)
528 !-------------------------------------------------------!
529 ! !
530 ! This subroutine calculates integral !
531 ! from 0 to r_n of f*r byinterpolating f !
532 ! with a local quadratic spline !
533 ! !
534 !-------------------------------------------------------!
535 
536  USE itm_types
537 
538  IMPLICIT NONE
539 
540  INTEGER i, n
541 
542  REAL (R8) h(n), r(n), f(n), int(n)
543  REAL (R8) den
544  REAL (R8) a(n), b(n)
545 
546 
547  DO i=2,n-1
548  den = h(i-1)*(h(i-1)+h(i))*h(i)
549  a(i) = ((f(i+1)-f(i))*h(i-1)**2-(f(i-1)-f(i))*h(i)**2)/den
550  b(i) = ((f(i+1)-f(i))*h(i-1)+(f(i-1)-f(i))*h(i))/den
551  END DO
552 
553 
554  int(1) = 0.d0
555 
556  int(2) = f(2)*r(2)**2/2.d0-a(2)*r(2)**3/6.d0+b(2)*r(2)**4/12.d0
557 
558 
559  DO i=2,n-2
560  int(i+1) = int(i)+(r(i)*f(i)+r(i+1)*f(i+1))*h(i)/2.d0+ &
561  (f(i)+r(i)*a(i)-f(i+1)-r(i+1)*a(i+1))*h(i)**2/8.d0+ &
562  (a(i)+r(i)*b(i)+a(i+1)+r(i+1)*b(i+1))*h(i)**3/24.d0+ &
563  (b(i)-b(i+1))*h(i)**4/64.d0
564  END DO
565 
566 
567  int(n) = int(n-1)+(f(n-1)-a(n-1)*r(n-1)+b(n-1)*r(n-1)**2)* &
568  (r(n)**2-r(n-1)**2)/2.d0+(a(n-1)-2.d0*b(n-1)*r(n-1))* &
569  (r(n)**3-r(n-1)**3)/3.d0+b(n-1)*(r(n)**4-r(n-1)**4)/4.d0
570 
571 
572  RETURN
573 
574 
575  END SUBROUTINE integral
576 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
577 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
578 
579 
580 
581 
582 
583 
584 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
585 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
586  SUBROUTINE axis(n, h, r, f)
587 !-------------------------------------------------------!
588 ! !
589 ! This subroutine finds !
590 ! f(r_1=0) from f(r_2), f(r_3) and f(r_4) !
591 ! !
592 !-------------------------------------------------------!
593 
594  USE itm_types
595 
596  IMPLICIT NONE
597 
598  INTEGER n
599  REAL (R8) h(n), r(n), f(n)
600 
601  f(1) = ((f(2)*r(4)/h(2)+f(4)*r(2)/h(3))*r(3) &
602  -f(3)*(r(2)/h(2)+r(4)/h(3))*r(2)*r(4)/r(3)) &
603  /(r(4)-r(2))
604 
605 
606  RETURN
607 
608 
609  END SUBROUTINE axis
610 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
611 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
612 
613 
derivn21
subroutine derivn21(N, X, Y, DY1)
Definition: solution2.f90:477
solver
Definition: solver.f90:1
type_solver
Definition: type_solver.f90:1
solut3
subroutine solut3(n, h, a, b, f, u, v, w, y, dy)
Definition: solution2.f90:312
axis
subroutine axis(n, h, r, f)
Definition: solution2.f90:586
phys_functions::p
real(r8) function p(a, x, xr, xs, yr, ys, psi, psir, F_dia)
Definition: phys_functions.f90:383
solver_ritm_integral
subroutine solver_ritm_integral(n, r, amix, z, zm, dz, aets, bets, cets, dets, eets, fets, gets, tau, uets, vets, wets)
Definition: solution2.f90:180
type_solver::numerics
Definition: type_solver.f90:22
derivn2
subroutine derivn2(N, X, Y, DY1)
Definition: solution2.f90:433
integral
subroutine integral(n, h, r, f, int)
Definition: solution2.f90:527
solution2
subroutine solution2(SOLVER, ifail)
NUMERICAL SOLUTION.
Definition: solution2.f90:13