ETS  \$Id: Doxyfile 2162 2020-02-26 14:16:09Z g2dpc $
 All Classes Files Functions Variables Pages
solution7.f90
Go to the documentation of this file.
1 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
7 ! + + + + + + + + + NUMERICAL SOLUTION + + + + + + + + +
8 
9 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
10 ! + + + + + + + + + NUMERICAL SOLUTION + + + + + + + + +
11 
12  SUBROUTINE solution7 (SOLVER, ifail)
13 
14 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
15 ! +
16 ! This subroutine is prepared to solve single +
17 ! transport equation in standardised form +
18 ! adopted by the ETS. +
19 ! +
20 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
21 ! Source: 1-D transport code RITM +
22 ! Developers: M.Tokar, D.Kalupin +
23 ! Kontacts: M.Tokar@fz-juelich.de +
24 ! D.Kalupin@fz-juelich.de +
25 ! +
26 ! Comments: equation is solved in +
27 ! DIFFERENTIAL form +
28 ! +
29 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
30 
31 
32  USE type_solver
33 
34  IMPLICIT NONE
35 
36  INTEGER, PARAMETER :: dp = kind(1.0d0) ! Double precision
37 
38  INTEGER, INTENT (INOUT) :: ifail
39 
40 ! +++ Input/Output with ETS:
41  TYPE (numerics), INTENT (INOUT) :: solver !contains all I/O quantities to numerics part
42 
43 
44 ! +++ Internal input/output parameters:
45  INTEGER :: idim, ndim !equation index and total number of equations to be solved
46 
47  INTEGER :: irho, nrho !radius index, number of radial points
48 
49  INTEGER :: flag !flag for equation: 0 - interpretative (not solved), 1 - predictive (solved)
50 
51  REAL (DP) :: rho(solver%nrho) !radii
52 
53  REAL (DP) :: amix !fraction of new sollution mixed
54 
55  REAL (DP) :: y(solver%nrho), ym(solver%nrho) !function at the current amd previous time steps
56  REAL (DP) :: dy(solver%nrho) !derivative of function
57 
58  REAL (DP) :: a(solver%nrho), b(solver%nrho) !coefficients
59  REAL (DP) :: c(solver%nrho), d(solver%nrho) !coefficients
60  REAL (DP) :: e(solver%nrho), f(solver%nrho) !coefficients
61  REAL (DP) :: g(solver%nrho), h !coefficients
62  REAL (DP) :: dd(solver%nrho), de(solver%nrho) !derivatives of coefficients
63 
64  REAL (DP) :: v(2), u(2), w(2) !boundary conditions
65 
66 
67 ! +++ Coefficients used by internal numerical solver:
68  INTEGER :: i, n !radius index, number of radial points
69  INTEGER :: isp !spline flag
70 
71  REAL (DP) :: x(solver%nrho) !radii
72 
73  REAL (DP) :: sol(solver%nrho) !solution
74  REAL (DP) :: dsol(solver%nrho) !derivative of solution
75 
76  REAL (DP) :: as(solver%nrho), bs(solver%nrho) !coefficients
77  REAL (DP) :: cs(solver%nrho) !coefficients
78 
79  REAL (DP) :: vs(2), us(2), ws(2) !boundary conditions
80 
81 
82 
83 ! + + + + + + + + + INTERFACE PART + + + + + + + + + + +
84 
85 
86 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
87 ! +++ Set up local variables (input, obtained from ETS):
88 ! Control parameters:
89  ndim = solver%NDIM
90  nrho = solver%NRHO
91  amix = solver%AMIX
92 
93 
94 
95 ! +++ Solution of equations starting from 1 to NDIM:
96  equation_loop1: DO idim = 1, ndim
97 
98  flag = solver%EQ_FLAG(idim)
99 
100 
101  IF (flag.EQ.0) goto 20 !equation is not solved
102 
103 
104 
105 ! +++ Numerical coefficients obtained from physics part in form:
106 !
107 ! (A*Y-B*Y(t-1))/H + 1/C * (-D*Y' + E*Y) = F - G*Y
108 
109  rho_loop1: DO irho=1,nrho
110  rho(irho) = solver%RHO(irho)
111 
112  y(irho) = solver%Y(idim,irho)
113  dy(irho) = solver%DY(idim,irho)
114  ym(irho) = solver%YM(idim,irho)
115 
116  a(irho) = solver%A(idim,irho)
117  b(irho) = solver%B(idim,irho)
118  c(irho) = solver%C(idim,irho)
119  d(irho) = solver%D(idim,irho)
120  e(irho) = solver%E(idim,irho)
121  f(irho) = solver%F(idim,irho)
122  g(irho) = solver%G(idim,irho)
123 
124  END DO rho_loop1
125 
126  h = solver%H
127 
128  CALL derivn7(nrho,rho,d,dd) !Derivation of coefficient D
129  CALL derivn7(nrho,rho,e,de) !Derivation of coefficient E
130 
131 ! +++ Boundary conditions for numerical solver in form:
132 !
133 ! V*Y' + U*Y = W
134 
135 ! +++ On axis:
136 
137  v(1) = solver%V(idim,1)
138  u(1) = solver%U(idim,1)
139  w(1) = solver%W(idim,1)
140 
141 ! +++ At the edge:
142 
143  v(2) = solver%V(idim,2)
144  u(2) = solver%U(idim,2)
145  w(2) = solver%W(idim,2)
146 
147 
148 
149 
150 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
151 ! +++ Set up numerical coefficients used in the solver
152 ! (translation is done for differential form
153 ! of transport equation):
154 
155 ! indexies:
156  isp = 1
157  n = nrho
158 ! coefficients:
159  rho_loop2: DO i = 2,n
160  x(i) = rho(i)
161  as(i) = -e(i)/d(i) + dd(i)/d(i)
162  bs(i) = a(i)*c(i)/h/d(i) + c(i)*g(i)/d(i) + de(i)/d(i)
163  cs(i) = c(i)/d(i) * (f(i) + b(i)*ym(i)/h)
164  END DO rho_loop2
165  x(1) = rho(1)
166  CALL axis7(n, x, as)
167  CALL axis7(n, x, bs)
168  CALL axis7(n, x, cs)
169 
170 ! boundary conditions:
171  vs(1) = v(1)
172  us(1) = u(1)
173  ws(1) = w(1)
174  vs(2) = v(2)
175  us(2) = u(2)
176  ws(2) = w(2)
177 
178 
179 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
180 ! +++ Solution of transport equation:
181 
182  CALL solutn7(x,n, vs,us,ws, as,bs,cs, sol,dsol)
183 ! radii b.c. coeff. solution
184 
185 
186  rho_loop3: DO irho = 1,nrho
187 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
188 ! +++ New function and derivative:
189  y(irho) = y(irho)*(1.d0-amix) + sol(irho)*amix
190  dy(irho) = dy(irho)*(1.d0-amix) + dsol(irho)*amix
191 
192 
193 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
194 ! +++ Return solution to ETS:
195  solver%Y(idim,irho) = y(irho)
196  solver%DY(idim,irho) = dy(irho)
197 
198  END DO rho_loop3
199 
200 
201 
202 
203 
204  20 CONTINUE
205 
206 
207  END DO equation_loop1
208 
209 
210 
211 
212  RETURN
213 
214 
215 
216  END SUBROUTINE solution7
217 
218 
219 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
220 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
221 
222 
223 
224 
225 
226 
227 
228 
229 ! + + + + + + + + + NUMERICAL PART + + + + + + + + + + +
230 
231 
232 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
233 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
234  SUBROUTINE solutn7 (RHO,N, V,U,W, A,B,F, Y,DY)
235 ! radii b.c. coeff. solution
236 !-------------------------------------------------------!
237 ! !
238 ! This subroutine provides solution for the equation!
239 ! !
240 ! d(dy/dx)/dx + a(x)*dy/dx = b(x)*y - f(x) !
241 ! !
242 !-------------------------------------------------------!
243  USE itm_types
244 
245  IMPLICIT NONE
246 
247  INTEGER i,j,k,l,m,n
248 
249  REAL (R8) rho(n),h(n+1),a(n),b(n),f(n),y(n),dy(n),lam(10001,2)
250  REAL (R8) u(2),v(2),w(2)
251  REAL (R8) y0(n),y1(10001,2,2)
252  REAL (R8) et(2,2),mu(2),det
253  REAL (R8) c(n),d(n),g(2),x
254  REAL (R8) cdy(n,2)
255 
256  h = 0.d0
257 
258  DO i=1,n-1
259  h(i) = rho(i+1)-rho(i)
260  END DO
261  h(n) = h(n-1)
262 
263  DO i=2,n
264  y0(i) = f(i)/b(i)
265  det = a(i)**2/4.d0+b(i)
266 
267  DO k=1,2
268  lam(i,k) = -a(i)/2.d0-(-1)**k*dsqrt(det)
269 
270 
271  FORALL(l=1:2) &
272  y1(i,k,l) = dexp(dmin1(200.d0,(-1)**l*lam(i,k) &
273  *h(i-2+l)/2.d0))
274  END DO
275 
276  END DO
277 
278 
279  c(3) = v(1)/(v(1)-h(3)*u(1))
280  d(3) = h(3)*w(1)/(h(3)*u(1)-v(1))
281 
282 
283  DO i=3,n-1
284 
285  DO m=1,2
286  l = 3-m
287  j = i+(-1)**m
288 
289  DO k=1,2
290  et(k,m) = -(-1)**m*(lam(j,m)-lam(i,k) &
291  -(lam(j,l)-lam(i,k))*y1(j,l,l) &
292  /y1(j,m,l))/(lam(j,1)-lam(j,2))
293  END DO
294 
295  mu(m) = -(-1)**m*(y0(i)-y0(j))/y1(i,l,m) &
296  *(lam(j,m)-lam(j,l)*y1(j,l,l) &
297  /y1(j,m,l))/(lam(j,1)-lam(j,2))
298  END DO
299 
300  det = et(1,2)*et(2,1)-et(1,1)*et(2,2) &
301  *y1(i,l,1)/y1(i,2,1) &
302  *y1(i,2,2)/y1(i,l,2)
303 
304  dy(i) = 0.0d0
305 
306  DO m=1,2
307  l = 3-m
308  g(m) = y1(i+(-1)**m,l,l)/y1(i,l,m)* &
309  (et(m,l)-y1(i,l,l)/y1(i,m,l) &
310  *et(l,l))/det
311 ! derivative:
312  dy(i) = dy(i)+lam(i,m)* &
313  (mu(m)*y1(i,l,l)/y1(i,m,l)* &
314  et(l,l)-mu(l)*et(l,m))/det
315  cdy(i,m) = (et(m,l)*lam(i,l)- &
316  y1(i,l,l)/y1(i,m,l)*et(l,l)*lam(i,m))/det
317  END DO
318 
319  x = (mu(1)*(y1(i,2,2)/y1(i,1,2) &
320  *et(2,2)-et(1,2))+mu(2)*(y1(i,1,1) &
321  /y1(i,2,1)*et(1,1)-et(2,1)))/det- &
322  y0(i-1)*g(1)-y0(i+1)*g(2)+y0(i)
323  c(i+1) = g(2)/(1.d0-g(1)*c(i))
324  d(i+1) = (g(1)*d(i)+x)/(1.d0-g(1)*c(i))
325 
326  END DO
327 
328  y(n) = (h(n-1)*w(2)+v(2)*d(n)) &
329  /(h(n-1)*u(2)+v(2)-v(2)*c(n))
330 
331  DO i=n,3,-1
332  y(i-1) = c(i)*y(i)+d(i)
333  END DO
334 
335 
336 ! derivative:
337  DO i=3,n-1
338  DO m=1,2
339  l = 3-m
340  j = i+(-1)**m
341  dy(i) = dy(i)+(y(j)-y0(j))*y1(j,l,l)/y1(i,l,m)*cdy(i,m)
342  END DO
343  END DO
344 
345  dy(2)=(y(3)-y(3))/h(2)
346  dy(n)=(y(n)-y(n-1))/h(n-1)
347 
348 
349 ! CALL DERIVN7(N,RHO,y,dy)
350 
351  RETURN
352 
353 
354  END SUBROUTINE solutn7
355 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
356 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
357 
358 
359 
360 
361 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
362 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
363 ! These subroutines calculate first and second derivatives,
364 ! DY1 and DY2, of function Y respect to argument X
365  SUBROUTINE derivn7(N,X,Y,DY1)
366 
367  IMPLICIT NONE
368 
369  INTEGER, PARAMETER :: dp = kind(1.0d0) ! Double precision
370  INTEGER :: n ! number of radial points (input)
371  INTEGER :: i
372 
373  REAL (DP) :: x(n), & ! argument array (input)
374  y(n), & ! function array (input)
375  dy1(n) ! function derivative array (output)
376  REAL (DP) :: h(n),dy2(n)
377 
378  DO i=1,n-1
379  h(i)=x(i+1)-x(i)
380  END DO
381 
382  DO i=2,n-1
383  dy1(i)=((y(i+1)-y(i))*h(i-1)/h(i)+(y(i)-y(i-1))*h(i)/h(i-1)) &
384  /(h(i)+h(i-1))
385  dy2(i)=2.d0*((y(i-1)-y(i))/h(i-1)+(y(i+1)-y(i))/h(i)) &
386  /(h(i)+h(i-1))
387  END DO
388 
389  dy1(1)=dy1(2)-dy2(2)*h(1)
390  dy1(n)=dy1(n-1)+dy2(n-1)*h(n-1)
391 
392  RETURN
393  END SUBROUTINE derivn7
394 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
395 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
396 
397 
398 
399 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
400 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
401  SUBROUTINE axis7(n, r, f)
402 !-------------------------------------------------------!
403 ! !
404 ! This subroutine finds !
405 ! f(r_1=0) from f(r_2), f(r_3) and f(r_4) !
406 ! !
407 !-------------------------------------------------------!
408 
409  IMPLICIT NONE
410 
411  INTEGER n, i
412  REAL *8 h(n), r(n), f(n)
413 
414 
415  DO i=1,3
416  h(i)=r(i+1)-r(i)
417  END DO
418 
419  f(1) = ((f(2)*r(4)/h(2)+f(4)*r(2)/h(3))*r(3) &
420  -f(3)*(r(2)/h(2)+r(4)/h(3))*r(2)*r(4)/r(3)) &
421  /(r(4)-r(2))
422 
423 
424  RETURN
425 
426 
427  END SUBROUTINE axis7
428 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
429 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
solver
Definition: solver.f90:1
type_solver
Definition: type_solver.f90:1
solutn7
subroutine solutn7(RHO, N, V, U, W, A, B, F, Y, DY)
Definition: solution7.f90:234
type_solver::numerics
Definition: type_solver.f90:22
axis7
subroutine axis7(n, r, f)
Definition: solution7.f90:401
solution7
subroutine solution7(SOLVER, ifail)
This subroutine is prepared to solve single transport equation in standardised form adopted by the ET...
Definition: solution7.f90:12
derivn7
subroutine derivn7(N, X, Y, DY1)
Definition: solution7.f90:365