ETS  \$Id: Doxyfile 2162 2020-02-26 14:16:09Z g2dpc $
 All Classes Files Functions Variables Pages
solution4.f90
Go to the documentation of this file.
1 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
7 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
8 SUBROUTINE solution4 (SOLVER,ifail)
9 
10 !-------------------------------------------------------!
11 ! This subroutine solves transport equations !
12 ! using matrix PROGONKA method !
13 !-------------------------------------------------------!
14 ! Source: -- !
15 ! Developers: R.Stankiewicz !
16 ! Contacts: romsta@ifpilm.waw.pl !
17 ! !
18 ! Comments: -- !
19 ! !
20 !-------------------------------------------------------!
21 
22 
23 ! +++ Declaration of variables:
24  use itm_types
25  USE type_solver
26 
27  IMPLICIT NONE
28 
29  TYPE(numerics) :: solver
30 
31  INTEGER :: ifail
32  INTEGER :: ndim,idim1,idim2,iupwind
33  INTEGER :: nrho,irho1,irho2,irhop,irhom
34 
35  REAL (R8) :: tau,h12
36  REAL (R8) :: rho(solver%nrho)
37 
38  REAL (R8) :: rho_derivn(solver%nrho) !AF 4.Oct.2011
39 
40  REAL (R8) :: y(solver%ndim,solver%nrho)
41  REAL (R8) :: ym(solver%ndim,solver%nrho)
42  REAL (R8) :: dy(solver%ndim,solver%nrho)
43  REAL (R8) :: a(solver%ndim,solver%nrho)
44  REAL (R8) :: b(solver%ndim,solver%nrho)
45  REAL (R8) :: c(solver%ndim,solver%nrho)
46  REAL (R8) :: d(solver%ndim,solver%nrho)
47  REAL (R8) :: e(solver%ndim,solver%nrho)
48  REAL (R8) :: f(solver%ndim,solver%nrho)
49  REAL (R8) :: g(solver%ndim,solver%nrho)
50  REAL (R8) :: cm1(solver%ndim,solver%ndim,solver%nrho)
51  REAL (R8) :: u(solver%ndim,2)
52  REAL (R8) :: v(solver%ndim,2)
53  REAL (R8) :: w(solver%ndim,2)
54 
55  REAL (R8) :: aa(solver%nrho,solver%ndim,solver%ndim)
56  REAL (R8) :: bb(solver%nrho,solver%ndim,solver%ndim)
57  REAL (R8) :: cc(solver%nrho,solver%ndim,solver%ndim)
58  REAL (R8) :: ff(solver%nrho,solver%ndim)
59  REAL (R8) :: h(solver%nrho)
60  REAL (R8) :: uu(2,solver%ndim,solver%ndim)
61  REAL (R8) :: vv(2,solver%ndim,solver%ndim)
62  REAL (R8) :: ww(2,solver%ndim)
63 
64  REAL (R8) :: yy(solver%nrho,solver%ndim)
65  REAL (R8) :: amix
66 
67 
68 
69 ! +++ Initial coefficients are set to zero:
70  iupwind = 0
71 
72  aa(:,:,:) = 0.0_r8
73  bb(:,:,:) = 0.0_r8
74  cc(:,:,:) = 0.0_r8
75  cm1(:,:,:) = 0.0_r8
76  ff(:,:) = 0.0_r8
77 
78  uu(:,:,:) = 0.0_r8
79  vv(:,:,:) = 0.0_r8
80  ww(:,:) = 0.0_r8
81 
82  yy(:,:) = 0.0_r8
83 
84  ndim = solver%NDIM
85  nrho = solver%NRHO
86  amix = solver%AMIX
87 
88 
89 ! +++ Input coefficients from the work flow:
90  rho_loop1: DO irho1 = 1,nrho
91  rho(irho1) = solver%RHO(irho1)
92 
93  dim_loop1: DO idim1 = 1,ndim
94  y(idim1,irho1) = solver%Y(idim1,irho1)
95  ym(idim1,irho1) = solver%YM(idim1,irho1)
96  dy(idim1,irho1) = solver%DY(idim1,irho1)
97 
98  a(idim1,irho1) = solver%A(idim1,irho1)
99  b(idim1,irho1) = solver%B(idim1,irho1)
100  c(idim1,irho1) = solver%C(idim1,irho1)
101  d(idim1,irho1) = solver%D(idim1,irho1)
102  e(idim1,irho1) = solver%E(idim1,irho1)
103  f(idim1,irho1) = solver%F(idim1,irho1)
104  g(idim1,irho1) = solver%G(idim1,irho1)
105 
106  dim_loop2: DO idim2 = 1,ndim
107  cm1(idim1,idim2,irho1) = solver%CM1(idim1,idim2,irho1)
108  END DO dim_loop2
109 
110  u(idim1,1) = solver%U(idim1,1)
111  v(idim1,1) = solver%V(idim1,1)
112  w(idim1,1) = solver%W(idim1,1)
113  u(idim1,2) = solver%U(idim1,2)
114  v(idim1,2) = solver%V(idim1,2)
115  w(idim1,2) = solver%W(idim1,2)
116 
117  END DO dim_loop1
118 
119  END DO rho_loop1
120  tau = solver%H
121 
122 
123 ! +++ Set up coefficients for matrix solver:
124  DO idim2 = 1,ndim
125  DO idim1 = 1,2
126  ww(idim1,idim2) = w(idim2,idim1)
127  ENDDO
128  ENDDO
129 
130 
131  DO irho1 = 2,nrho
132  irhom = irho1-1
133  h(irhom) = rho(irho1)-rho(irhom)
134  ENDDO
135  DO irho1=1,nrho
136  DO idim1=1,ndim
137  DO idim2=1,ndim
138  f(idim2,irho1)=f(idim2,irho1)!-CM1(IDIM1,IDIM2,IRHO1)*YM(IDIM1,IRHO1)
139  enddo
140  enddo
141  ENDDO
142  DO idim1 = 1,ndim
143  uu(1,idim1,idim1) = - 0.5*v(idim1,1)/h(1)+0.5*u(idim1,1)
144  vv(1,idim1,idim1) = 0.5*v(idim1,1)/h(1)+0.5*u(idim1,1)
145  uu(2,idim1,idim1) = 0.5* v(idim1,2)/h(nrho-1)+0.5*u(idim1,2)
146  vv(2,idim1,idim1) = - 0.5*v(idim1,2)/h(nrho-1)+0.5*u(idim1,2)
147  ENDDO
148  irho1=2
149  irhop = irho1+1
150  irhom = irho1-1
151  h12 = 0.5*(h(irho1)+2.*h(irhom))
152 
153  doloop2idim1 : DO idim1 = 1,ndim
154  ff(irho1,idim1) = b(idim1,irho1)*ym(idim1,irho1)+f(idim1,irho1)*tau
155  aa(irho1,idim1,idim1) = -0.5/c(idim1,irho1)/h(irho1)/ h12 &
156  * (d(idim1,irhop)+d(idim1,irho1))
157  cc(irho1,idim1,idim1) = -0.5/c(idim1,irho1)/h(irhom)/ h12 &
158  * d(idim1,irhom)
159  bb(irho1,idim1,idim1) = -aa(irho1,idim1,idim1) -cc(irho1,idim1,idim1)
160 
161  IF (iupwind.NE.1) THEN
162  aa(irho1,idim1,idim1) = aa(irho1,idim1,idim1) &
163  + 0.25 /c(idim1,irho1)/h12 &
164  * (e(idim1,irhop)+e(idim1,irho1))
165  cc(irho1,idim1,idim1) = cc(irho1,idim1,idim1) &
166  - 0.5 /c(idim1,irho1)/h12 &
167  * e(idim1,irhom)
168  bb(irho1,idim1,idim1) = bb(irho1,idim1,idim1) &
169  + 0.25 /c(idim1,irho1)/h12 &
170 ! * (E(IDIM1,IRHO1)+E(IDIM1,IRHOM)-2.*E(IDIM1,IRHOM)) !AF 2.Dec.2011 - Roman found this bug that affects the axis
171  * (e(idim1,irho1)+e(idim1,irhop)-2.*e(idim1,irhom)) !AF 2.Dec.2011 - Roman found this bug that affects the axis
172  ELSE
173  IF (e(idim1,irho1).GE.0.) THEN
174  cc(irho1,idim1,idim1) = cc(irho1,idim1,idim1) &
175  - e(idim1,irhom)/c(idim1,irho1)/h12
176  bb(irho1,idim1,idim1) = bb(irho1,idim1,idim1) &
177  + e(idim1,irho1) / c(idim1,irho1)/h12
178  ELSE
179  aa(irho1,idim1,idim1) = aa(irho1,idim1,idim1) &
180  + e(idim1,irhop)/c(idim1,irhop)/h12
181  bb(irho1,idim1,idim1) = bb(irho1,idim1,idim1) &
182  + e(idim1,irho1) / c(idim1,irho1)/h12
183  ENDIF
184 
185  ENDIF
186 
187  bb(irho1,idim1,idim1) = bb(irho1,idim1,idim1)+g(idim1,irho1)
188 
189  DO idim2 = 1,ndim
190  bb(irho1,idim1,idim2) = bb(irho1,idim1,idim2)!- CM1(IDIM1,IDIM2,IRHO1)
191  ENDDO
192 
193 
194  ENDDO doloop2idim1
195  irho1=nrho-1
196 
197 
198  irhop = irho1+1
199  irhom = irho1-1
200  h12 = 0.5*(2.*h(irho1)+h(irhom))
201 
202  doloop2idim2 : DO idim1 = 1,ndim
203  ff(irho1,idim1) = b(idim1,irho1)*ym(idim1,irho1)+f(idim1,irho1)*tau
204  aa(irho1,idim1,idim1) = -0.5/c(idim1,irho1)/h(irho1)/ h12 &
205  * d(idim1,irhop)
206  cc(irho1,idim1,idim1) = -0.5/c(idim1,irho1)/h(irhom)/ h12 &
207  * (d(idim1,irho1)+ d(idim1,irhom))
208  bb(irho1,idim1,idim1) = -aa(irho1,idim1,idim1) -cc(irho1,idim1,idim1)
209 
210  IF (iupwind.NE.1) THEN
211  aa(irho1,idim1,idim1) = aa(irho1,idim1,idim1) &
212  + 0.5 /c(idim1,irho1)/h12 &
213  * e(idim1,irhop)
214  cc(irho1,idim1,idim1) = cc(irho1,idim1,idim1) &
215  - 0.25 /c(idim1,irho1)/h12 &
216  * (e(idim1,irho1)+e(idim1,irhom))
217  bb(irho1,idim1,idim1) = bb(irho1,idim1,idim1) &
218  + 0.25 /c(idim1,irho1)/h12 &
219  * (2.*e(idim1,irhop)-e(idim1,irhom)-e(idim1,irho1))
220  ELSE
221  IF (e(idim1,irho1).GE.0.) THEN
222  cc(irho1,idim1,idim1) = cc(irho1,idim1,idim1) &
223  - e(idim1,irhom)/c(idim1,irho1)/h12
224  bb(irho1,idim1,idim1) = bb(irho1,idim1,idim1) &
225  + e(idim1,irho1) / c(idim1,irho1)/h12
226  ELSE
227  aa(irho1,idim1,idim1) = aa(irho1,idim1,idim1) &
228  + e(idim1,irhop)/c(idim1,irhop)/h12
229  bb(irho1,idim1,idim1) = bb(irho1,idim1,idim1) &
230  + e(idim1,irho1) / c(idim1,irho1)/h12
231  ENDIF
232 
233  ENDIF
234 
235  bb(irho1,idim1,idim1) = bb(irho1,idim1,idim1)+g(idim1,irho1)
236 
237  DO idim2 = 1,ndim
238  bb(irho1,idim1,idim2) = bb(irho1,idim1,idim2)!- CM1(IDIM1,IDIM2,IRHO1)
239  ENDDO
240 
241 
242  ENDDO doloop2idim2
243 
244 
245  doloop1irho1: DO irho1 = 3,nrho-2
246  irhop = irho1+1
247  irhom = irho1-1
248  h12 = 0.5*(h(irho1)+h(irhom))
249 
250  doloop2idim3 : DO idim1 = 1,ndim
251  ff(irho1,idim1) = b(idim1,irho1)*ym(idim1,irho1)+f(idim1,irho1)*tau
252  aa(irho1,idim1,idim1) = -0.5/c(idim1,irho1)/h(irho1)/ h12 &
253  * (d(idim1,irho1)+ d(idim1,irhop))
254  cc(irho1,idim1,idim1) = -0.5/c(idim1,irho1)/h(irhom)/ h12 &
255  * (d(idim1,irho1)+ d(idim1,irhom))
256  bb(irho1,idim1,idim1) = -aa(irho1,idim1,idim1) -cc(irho1,idim1,idim1)
257 
258  IF (iupwind.NE.1) THEN
259  aa(irho1,idim1,idim1) = aa(irho1,idim1,idim1) &
260  + 0.25 /c(idim1,irho1)/h12 &
261  * (e(idim1,irho1)+e(idim1,irhop))
262  cc(irho1,idim1,idim1) = cc(irho1,idim1,idim1) &
263  - 0.25 /c(idim1,irho1)/h12 &
264  * (e(idim1,irho1)+e(idim1,irhom))
265  bb(irho1,idim1,idim1) = bb(irho1,idim1,idim1) &
266  + 0.25 /c(idim1,irho1)/h12 &
267  * (e(idim1,irhop)-e(idim1,irhom))
268  ELSE
269  IF (e(idim1,irho1).GE.0.) THEN
270  cc(irho1,idim1,idim1) = cc(irho1,idim1,idim1) &
271  - e(idim1,irhom)/c(idim1,irho1)/h12
272  bb(irho1,idim1,idim1) = bb(irho1,idim1,idim1) &
273  + e(idim1,irho1) / c(idim1,irho1)/h12
274  ELSE
275  aa(irho1,idim1,idim1) = aa(irho1,idim1,idim1) &
276  + e(idim1,irhop)/c(idim1,irhop)/h12
277  bb(irho1,idim1,idim1) = bb(irho1,idim1,idim1) &
278  + e(idim1,irho1) / c(idim1,irho1)/h12
279  ENDIF
280 
281  ENDIF
282 
283  bb(irho1,idim1,idim1) = bb(irho1,idim1,idim1)+g(idim1,irho1)
284 
285  DO idim2 = 1,ndim
286  bb(irho1,idim1,idim2) = bb(irho1,idim1,idim2)!- CM1(IDIM1,IDIM2,IRHO1)
287  ENDDO
288 
289 
290  ENDDO doloop2idim3
291 
292 
293  ENDDO doloop1irho1
294 
295 
296 
297 
298 !------------------------------COUPLING TERMS ----------------------------
299 ! DO IRHO1 = 2,NRHO-1
300 ! IRHOP = IRHO1+1
301 ! IRHOM = IRHO1-1
302 ! DO IDIM1 = 1,NDIM
303 ! H12 = 0.5*(H(IRHO1)+H(IRHOM))
304 ! if(IUPWIND.NE.1) THEN
305 ! AA(IRHO1,NDIM,IDIM1) = AA(IRHO1,NDIM,IDIM1)+ &
306 ! 0.25 *S2(IRHO1)/H12*(S3(IDIM1,IRHO1)+S3(IDIM1,IRHOP))
307 ! CC(IRHO1,NDIM,IDIM1) = CC(IRHO1,NDIM,IDIM1)- &
308 ! -0.25 *S2(IRHO1)/H12*(S3(IDIM1,IRHO1)+S3(IDIM1,IRHOM))
309 ! BB(IRHO1,NDIM,IDIM1) = BB(IRHO1,NDIM,IDIM1) +&
310 ! 0.25 /S2(IRHO1)/H12*(S3(IDIM1,IRHOP)-S3(IDIM1,IRHOM))
311 ! ELSE
312 ! IF(E(IDIM1,IRHO1).GE.0.) THEN
313 ! CC(IRHO1,NDIM,IDIM1) = CC(IRHO1,NDIM,IDIM1) - &
314 ! S2(IRHOM)*S3(IDIM1,IRHOM)/H12
315 ! BB(IRHO1,NDIM,IDIM1) = BB(IRHO1,NDIM,IDIM1)+ &
316 ! S2(IRHO1) * S3(IDIM1,IRHO1)/H12
317 ! ELSE
318 ! AA(IRHO1,NDIM,IDIM1) = AA(IRHO1,NDIM,IDIM1)+ &
319 ! S2(IRHOP)*S3(IDIM1,IRHOP)/H12
320 ! BB(IRHO1,NDIM,IDIM1) = BB(IRHO1,NDIM,IDIM1)+ &
321 ! S2(IRHO1) * S3(IDIM1,IRHO1)/H12
322 ! ENDIF
323 ! ENDIF
324 
325 
326 
327  bb = bb*tau
328  aa = aa*tau
329  cc = cc*tau
330 
331 
332  DO irho1 = 1,nrho
333  DO idim1 = 1,ndim
334  bb(irho1,idim1,idim1) = bb(irho1,idim1,idim1)+ a(idim1,irho1)
335  ENDDO
336  ENDDO
337 
338 
339 
340 ! +++ Call the matrix progonka:
341  CALL mprogonka_4(nrho,ndim,aa,bb,cc,ff,uu,vv,ww,yy)
342 
343 
344 ! +++ New solution:
345  DO idim1 = 1,ndim
346  DO irho1 = 1,nrho
347  y(idim1,irho1) = yy(irho1,idim1)*amix+y(idim1,irho1)*(1.e0_r8-amix)
348  ENDDO
349  ENDDO
350 
351  !AF 4.Oct.2011 - at this stage IRHO.eq.1 and IRHO.eq.NRHO are still the ghost points, so DERIVN needs the corresponding (regular) RHO grid,
352  !otherwise DY at points 2 and N-1 will be wrong
353  rho_derivn = rho
354  rho_derivn(1) = rho(2) - ( rho(3) - rho(2) )
355  rho_derivn(nrho) = rho(nrho-1) + ( rho(nrho-1) - rho(nrho-2) )
356  !AF - End
357 
358 ! CALL DERIVN3_4(NRHO,NDIM,RHO,Y,DY) !AF 4.Oct.2011
359  CALL derivn3_4(nrho,ndim,rho_derivn,y,dy) !AF 4.Oct.2011
360  DO idim1=1,ndim
361 ! DY(IDIM1,1)=0.5*(Y(IDIM1,2)-Y(IDIM1,1))/H(1) !AF 6.Oct.2011
362  dy(idim1,1)=0.5*(dy(idim1,2)+dy(idim1,1)) !AF 6.Oct.2011 - this uses DY from DERIVN3_4 - seems to make no difference
363  y(idim1,1)=0.5*(y(idim1,2)+y(idim1,1))
364  !DY(IDIM1,NRHO)=0.5*(Y(IDIM1,NRHO-1)-Y(IDIM1,NRHO))/H(NRHO-1) !AF 23.Sep.2011
365 ! DY(IDIM1,NRHO)=0.5*(Y(IDIM1,NRHO)-Y(IDIM1,NRHO-1))/H(NRHO-1) !AF 23.Sep.2011, 6.Oct.2011
366  dy(idim1,nrho)=0.5*(dy(idim1,nrho)+dy(idim1,nrho-1)) !AF 23.Sep.2011, 6.Oct.2011 - this uses DY from DERIVN3_4 - seems to make no difference
367  y(idim1,nrho)=0.5*(y(idim1,nrho)+y(idim1,nrho-1))
368  ENDDO
369 ! +++ Return solution to the work flow:
370  DO idim1 = 1,ndim
371  DO irho1 = 1,nrho
372  solver%Y(idim1,irho1) = y(idim1,irho1)
373  solver%DY(idim1,irho1) = dy(idim1,irho1)
374  ENDDO
375  ENDDO
376 
377 
378 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
379 ! +++ Return solution to ETS:
380 
381 
382  RETURN
383 
384 
385 
386 END SUBROUTINE solution4 !AF, 14.Sep.2011
387 !END SUBROUTINE SOLUTION3 !AF, 14.Sep.2011
388 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
389 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
390 
391 
392 
393 
394 
395 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
407 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
408 SUBROUTINE mprogonka_4 (NP,NDIM,A,B,C,F,U,V,W,Y) !AF, 14.Sep.2011
409 !SUBROUTINE MPROGONKA (NP,NDIM,A,B,C,F,U,V,W,Y) !AF, 14.Sep.2011
410 ! Finds Solution of discrete equation of the form
411 ! A_i * Y_i+1 + B_i * Y_i +C_i * Y_i-1 = F_i
412 ! where:
413 ! Y_i - unknown wektor at point i (Y_i(NDIM))
414 ! F_i - Rhs at point i (F_i(NDIM))
415 ! A_i,B_i,C_i - equation coefficients (matrixes (NDIM,NDIM))
416 ! U,V,W - matrixes with boundary conditions
417 
418 
419  use itm_types
420  IMPLICIT NONE
421 !! Number of points
422  INTEGER, INTENT(IN) :: np
423 !! Rank of the block matrix
424  INTEGER, INTENT(IN) :: ndim
425 !! Equation coefficients: A, B, C
426  REAL(R8), DIMENSION(NP,NDIM,NDIM), INTENT(IN) :: a, b, c
427 !! Right hand side: F
428  REAL(R8), DIMENSION(NP,NDIM), INTENT(IN) :: f
429 !! Boundary conditions arrays U,V
430  REAL(R8), DIMENSION(2,NDIM,NDIM), INTENT(IN) :: u,v
431 !! Boundary condition array W
432  REAL(R8), DIMENSION(2,NDIM), INTENT(IN) :: w
433 !! Solution: Y
434  REAL(R8), DIMENSION(NP,NDIM), INTENT(OUT) :: y
435 
436 ! working variables
437 !! Progonka coefficients: alpha, beta
438  REAL(R8), DIMENSION(NP,NDIM,NDIM) :: alf
439  REAL(R8), DIMENSION(NP,NDIM) :: bet
440  REAL(R8), DIMENSION(NDIM,NDIM) :: temp1,temp2,temp3
441  REAL(R8), DIMENSION(NP,NDIM) :: wekt1
442  REAL(R8), DIMENSION(NDIM) :: wekt2
443  INTEGER :: i, j, j2, k, n, nn, n2
444  REAL(R8) :: stemp1,stemp2,det
445 
446 
447 
448 
449  IF(ndim .EQ. 1) THEN
450 ! BOUNDARY CONDITION
451  alf(2,1,1) =-v(1,1,1)/u(1,1,1)
452  bet(2,1)= w(1,1)/u(1,1,1)
453 
454  DO j=2,np-1
455  j2=j+1
456  det=b(j,1,1)+alf(j,1,1)*c(j,1,1)
457  IF(abs(det).GT.1.e-30_r8) THEN
458  alf(j2,1,1)=-a(j,1,1)/det
459  bet(j2,1)=(f(j,1)-bet(j,1)*c(j,1,1))/det
460  ELSE
461  alf(j2,1,1)=0.
462  bet(j2,1)=0.
463  ENDIF
464  ENDDO
465 
466 ! BOUNDARY CONDITION
467  y(np,1)=0.
468  det=v(2,1,1)*alf(np,1,1)+u(2,1,1)
469  IF(abs(det).GT.1.e-30_r8) y(np,1)=(w(2,1)-bet(np,1)*v(2,1,1))/det
470 
471  DO j=np-1,1,-1
472  j2=j+1
473  y(j,1)=alf(j2,1,1)*y(j2,1)+bet(j2,1)
474  ENDDO
475 
476 
477  ELSE
478 ! BOUNDARY CONDITION
479  DO i=1,ndim
480  DO j=1,ndim
481  temp1(i,j)=u(1,i,j)
482  ENDDO
483  ENDDO
484  CALL invmatrix_4(ndim, temp1, temp2)
485 
486  DO i=1,ndim
487  bet(2,i)=0.0
488  DO j=1,ndim
489  alf(2,i,j)=0.0
490  DO k=1,ndim
491  alf(2,i,j)=alf(2,i,j)-temp2(i,k)*v(1,k,j)
492  ENDDO
493  bet(2,i)=bet(2,i)+temp2(i,j)*w(1,j)
494  ENDDO
495  ENDDO
496 
497  DO n=2,np-1
498  DO i=1,ndim
499  wekt1(n,i)=f(n,i)
500  DO j=1,ndim
501  temp1(i,j)=b(n,i,j)
502  DO k=1,ndim
503  temp1(i,j)=temp1(i,j)+c(n,i,k)*alf(n,k,j)
504  ENDDO
505  wekt1(n,i)=wekt1(n,i)-c(n,i,j)*bet(n,j)
506  ENDDO
507  ENDDO
508  CALL invmatrix_4(ndim, temp1, temp2)
509 
510  n2=n+1
511  DO i=1,ndim
512  bet(n2,i)=0.
513  DO j=1,ndim
514  alf(n2,i,j)=0.
515  DO k=1,ndim
516  alf(n2,i,j)=alf(n2,i,j)-temp2(i,k)*a(n,k,j)
517  ENDDO
518  bet(n2,i)=bet(n2,i)+temp2(i,j)*wekt1(n,j)
519  ENDDO
520  ENDDO
521  ENDDO
522 
523 
524 1001 FORMAT(1x,3i4,1p,8e12.4)
525 
526 
527 ! BOUNDARY CONDITION
528  DO i=1,ndim
529  wekt2(i)=w(2,i)
530  DO j=1,ndim
531  temp1(i,j)=u(2,i,j)
532  DO k=1,ndim
533  temp1(i,j)=temp1(i,j)+v(2,i,k)*alf(np,k,j)
534  ENDDO
535  wekt2(i)=wekt2(i)-v(2,i,j)*bet(np,j)
536  ENDDO
537  ENDDO
538  CALL invmatrix_4(ndim, temp1, temp2)
539 
540  DO i=1,ndim
541  y(np,i)=0.
542  DO j=1,ndim
543  y(np,i)=y(np,i)+temp2(i,j)*wekt2(j)
544  ENDDO
545  ENDDO
546 
547  DO n=np-1,1,-1
548  n2=n+1
549  DO i=1,ndim
550  y(n,i)=bet(n2,i)
551  DO j=1,ndim
552  y(n,i)=y(n,i)+alf(n2,i,j)*y(n2,j)
553  ENDDO
554  ENDDO
555  ENDDO
556 
557 
558  ENDIF
559 
560 
561  RETURN
562 
563 END SUBROUTINE mprogonka_4 !AF, 14.Sep.2011
564 !END SUBROUTINE MPROGONKA !AF, 14.Sep.2011
565 
566 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
571 ! Outputs
577 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
578 
579 SUBROUTINE invmatrix_4(n, AA, Ainv) !AF, 14.Sep.2011
580 !SUBROUTINE INVMATRIX(n, AA, Ainv) !AF, 14.Sep.2011
581 
582 ! Computes inverse of matrix
583 ! Input
584 !-------------------------------------------------
585 ! AA - Matrix A (n by n)
586 ! N - Dimension of matrix AA
587 !-------------------------------------------------
588 ! Outputs
589 ! Ainv - Inverse of matrix AA (n by n)
590 !-------------------------------------------------
591 
592  use itm_types
593 
594  IMPLICIT NONE
595 
596 !! order of matrix
597  INTEGER, INTENT(IN) :: n
598 !! Input Matrix
599  REAL(R8), DIMENSION(n,n), INTENT(IN) :: aa
600 
601 !! Output Inverted Matrix
602  REAL(R8), DIMENSION(n,n), INTENT(OUT) :: ainv
603 
604 
605  REAL(R8), DIMENSION(n) :: scale ! Scale factor
606  REAL(R8), DIMENSION(n,n) :: b ! Work array
607  REAL(R8), DIMENSION(n,n) :: a ! Working copy of input matrix
608  INTEGER, DIMENSION(n) :: index ! Index Vector
609  INTEGER :: i, j, k, jpivot, indexj
610  REAL(R8) :: scalemax, ratio, ratiomax, coeff, determ, sum
611 
612 
613 ! Copy matrix A so as not to modify original
614  do i=1,n
615  do j=1,n
616  a(i,j) = aa(i,j)
617  enddo
618  enddo
619 
620 !* Matrix b is initialized to the identity matrix
621  do i=1,n
622  do j=1,n
623  if( i .eq. j ) then
624  b(i,j) = 1.0
625  else
626  b(i,j) = 0.0
627  endif
628  enddo
629  enddo
630 
631 !* Set scale factor, scale(i) = max( |a(i,j)| ), for each row
632  do i=1,n
633  index(i) = i ! Initialize row index list
634  scalemax = 0.0
635  do j=1,n
636  if( abs(a(i,j)) .gt. scalemax ) then
637  scalemax = abs(a(i,j))
638  endif
639  enddo
640  scale(i) = scalemax
641  enddo
642 
643 !* Loop over rows k = 1, ..., (N-1)
644  do k=1,n-1
645 !* Select pivot row from max( |a(j,k)/s(j)| )
646  ratiomax = 0.0
647  jpivot = k
648  do i=k,n
649  ratio = abs(a(index(i),k))/scale(index(i))
650  if( ratio .gt. ratiomax ) then
651  jpivot=i
652  ratiomax = ratio
653  endif
654  enddo
655 !* Perform pivoting using row index list
656  indexj = index(k)
657  if( jpivot .ne. k ) then ! Pivot
658  indexj = index(jpivot)
659  index(jpivot) = index(k) ! Swap index jPivot and k
660  index(k) = indexj
661  endif
662 !* Perform forward elimination
663  do i=k+1,n
664  coeff = a(index(i),k)/a(indexj,k)
665  do j=k+1,n
666  a(index(i),j) = a(index(i),j) - coeff*a(indexj,j)
667  enddo
668  a(index(i),k) = coeff
669  do j=1,n
670  b(index(i),j) = b(index(i),j) - a(index(i),k)*b(indexj,j)
671  enddo
672  enddo
673  enddo
674 
675 
676 !* Perform backsubstitution
677  do k=1,n
678  ainv(n,k) = b(index(n),k)/a(index(n),n)
679  do i=n-1,1,-1
680  sum = b(index(i),k)
681  do j=i+1,n
682  sum = sum - a(index(i),j)*ainv(j,k)
683  enddo
684  ainv(i,k) = sum/a(index(i),i)
685  enddo
686  enddo
687 
688  return
689 
690 END SUBROUTINE invmatrix_4 !AF, 14.Sep.2011
691 !END SUBROUTINE INVMATRIX !AF, 14.Sep.2011
692 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
693 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
694 
695 
696 
697 
698 
699 
700 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
707 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
708 SUBROUTINE derivn3_4(N,NFUN,X,FUN,DFUN) !AF, 14.Sep.2011
709 !SUBROUTINE DERIVN3(N,NFUN,X,FUN,DFUN) !AF, 14.Sep.2011
710 
711  use itm_types
712 
713  IMPLICIT NONE
714 
715  INTEGER :: n ! number of radial points (input)
716  INTEGER :: nfun ! number of functions
717  INTEGER :: i, ifun
718 
719  REAL (R8) :: x(n), & ! argument array (input)
720  fun(nfun,n), & ! function array (input)
721  dfun(nfun,n), & ! function derivative array (output)
722  y(n), & ! function array 1-D (local)
723  dy1(n) ! function derivative array 1-D (local)
724  REAL (R8) :: h(n),dy2(n)
725 
726  REAL (R8) :: ddy !AF 6.Oct.2011
727 
728 
729 
730  DO i=1,n-1
731  h(i) = x(i+1)-x(i)
732  END DO
733 
734  DO ifun = 1,nfun
735 
736  DO i = 1,n
737  y(i) = fun(ifun,i)
738  END DO
739 
740  DO i=2,n-1
741  dy1(i) = ((y(i+1)-y(i))*h(i-1)/h(i)+(y(i)-y(i-1))*h(i)/h(i-1)) &
742  /(h(i)+h(i-1))
743 ! DY2(I) = 2.e0_R8*((Y(I-1)-Y(I))/H(I-1)+(Y(I+1)-Y(I))/H(I)) & !AF 6.Oct.2011
744 ! /(H(I)+H(I-1))
745  END DO
746 
747 ! DY1(1) = DY1(2)-DY2(2)*H(1) !AF 6.Oct.2011
748 ! DY1(N) = DY1(N-1)+DY2(N-1)*H(N-1) !AF 6.Oct.2011
749 
750  ddy = 2.e0_r8*((y(1)-y(2))/h(1)+(y(3)-y(2))/h(2))/(h(2)+h(1))
751  dy1(1) = dy1(2)-ddy*h(1)
752  ddy = 2.e0_r8*((y(n-2)-y(n-1))/h(n-2)+(y(n)-y(n-1))/h(n-1))/(h(n-1)+h(n-2))
753  dy1(n) = dy1(n-1)+ddy*h(n-1)
754 
755  DO i = 1,n
756  dfun(ifun,i) = dy1(i)
757  END DO
758 
759 
760  END DO
761 
762  RETURN
763 
764 
765 END SUBROUTINE derivn3_4 !AF, 14.Sep.2011
766 !END SUBROUTINE DERIVN3 !AF, 14.Sep.2011
767 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
768 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
solver
Definition: solver.f90:1
type_solver
Definition: type_solver.f90:1
fun
subroutine fun(X, F)
Definition: Ev2.f:10
derivn3_4
subroutine derivn3_4(N, NFUN, X, FUN, DFUN)
These subroutines calculate first and second derivatives, DY1 and DY2, of function Y respect to argum...
Definition: solution4.f90:708
invmatrix_4
subroutine invmatrix_4(n, AA, Ainv)
Computes inverse of matrix Input AA - Matrix A (n by n) N - Dimension of matrix AA.
Definition: solution4.f90:579
mprogonka_4
subroutine mprogonka_4(NP, NDIM, A, B, C, F, U, V, W, Y)
Finds Solution of discrete equation of the form A_i * Y_i+1 + B_i * Y_i +C_i * Y_i-1 = F_i where: Y_i...
Definition: solution4.f90:408
solution4
subroutine solution4(SOLVER, ifail)
This subroutine solves transport equations using matrix PROGONKA method.
Definition: solution4.f90:8
phys_functions::p
real(r8) function p(a, x, xr, xs, yr, ys, psi, psir, F_dia)
Definition: phys_functions.f90:383
type_solver::numerics
Definition: type_solver.f90:22