ETS  \$Id: Doxyfile 2162 2020-02-26 14:16:09Z g2dpc $
 All Classes Files Functions Variables Pages
solution3.f90
Go to the documentation of this file.
1 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
7 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
8 SUBROUTINE solution3 (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  REAL (R8) :: y(solver%ndim,solver%nrho)
38  REAL (R8) :: ym(solver%ndim,solver%nrho)
39  REAL (R8) :: dy(solver%ndim,solver%nrho)
40  REAL (R8) :: a(solver%ndim,solver%nrho)
41  REAL (R8) :: b(solver%ndim,solver%nrho)
42  REAL (R8) :: c(solver%ndim,solver%nrho)
43  REAL (R8) :: d(solver%ndim,solver%nrho)
44  REAL (R8) :: e(solver%ndim,solver%nrho)
45  REAL (R8) :: f(solver%ndim,solver%nrho)
46  REAL (R8) :: g(solver%ndim,solver%nrho)
47  REAL (R8) :: cm1(solver%ndim,solver%ndim,solver%nrho)
48  REAL (R8) :: u(solver%ndim,2)
49  REAL (R8) :: v(solver%ndim,2)
50  REAL (R8) :: w(solver%ndim,2)
51 
52  REAL (R8) :: aa(solver%nrho,solver%ndim,solver%ndim)
53  REAL (R8) :: bb(solver%nrho,solver%ndim,solver%ndim)
54  REAL (R8) :: cc(solver%nrho,solver%ndim,solver%ndim)
55  REAL (R8) :: ff(solver%nrho,solver%ndim)
56  REAL (R8) :: h(solver%nrho)
57  REAL (R8) :: uu(2,solver%ndim,solver%ndim)
58  REAL (R8) :: vv(2,solver%ndim,solver%ndim)
59  REAL (R8) :: ww(2,solver%ndim)
60 
61  REAL (R8) :: yy(solver%nrho,solver%ndim)
62  REAL (R8) :: amix
63 
64 
65 
66 ! +++ Initial coefficients are set to zero:
67  iupwind = 0
68 
69  aa(:,:,:) = 0.0_r8
70  bb(:,:,:) = 0.0_r8
71  cc(:,:,:) = 0.0_r8
72  cm1(:,:,:) = 0.0_r8
73  ff(:,:) = 0.0_r8
74 
75  uu(:,:,:) = 0.0_r8
76  vv(:,:,:) = 0.0_r8
77  ww(:,:) = 0.0_r8
78 
79  yy(:,:) = 0.0_r8
80 
81  ndim = solver%NDIM
82  nrho = solver%NRHO
83  amix = solver%AMIX
84 
85 
86 ! +++ Input coefficients from the work flow:
87  rho_loop1: DO irho1 = 1,nrho
88  rho(irho1) = solver%RHO(irho1)
89 
90  dim_loop1: DO idim1 = 1,ndim
91  y(idim1,irho1) = solver%Y(idim1,irho1)
92  ym(idim1,irho1) = solver%YM(idim1,irho1)
93  dy(idim1,irho1) = solver%DY(idim1,irho1)
94 
95  a(idim1,irho1) = solver%A(idim1,irho1)
96  b(idim1,irho1) = solver%B(idim1,irho1)
97  c(idim1,irho1) = solver%C(idim1,irho1)
98  d(idim1,irho1) = solver%D(idim1,irho1)
99  e(idim1,irho1) = solver%E(idim1,irho1)
100  f(idim1,irho1) = solver%F(idim1,irho1)
101  g(idim1,irho1) = solver%G(idim1,irho1)
102 
103  dim_loop2: DO idim2 = 1,ndim
104  cm1(idim1,idim2,irho1) = solver%CM1(idim1,idim2,irho1)
105  END DO dim_loop2
106 
107  u(idim1,1) = solver%U(idim1,1)
108  v(idim1,1) = solver%V(idim1,1)
109  w(idim1,1) = solver%W(idim1,1)
110  u(idim1,2) = solver%U(idim1,2)
111  v(idim1,2) = solver%V(idim1,2)
112  w(idim1,2) = solver%W(idim1,2)
113 
114  END DO dim_loop1
115 
116  END DO rho_loop1
117  tau = solver%H
118 
119 
120 ! +++ Set up coefficients for matrix solver:
121  DO idim2 = 1,ndim
122  DO idim1 = 1,2
123  ww(idim1,idim2) = w(idim2,idim1)
124  ENDDO
125  ENDDO
126 
127 
128  DO irho1 = 2,nrho
129  irhom = irho1-1
130  h(irhom) = rho(irho1)-rho(irhom)
131  ENDDO
132 
133 
134  DO idim1 = 1,ndim
135  uu(1,idim1,idim1) = - v(idim1,1)/h(1)+u(idim1,1)
136  vv(1,idim1,idim1) = v(idim1,1)/h(1)
137  uu(2,idim1,idim1) = v(idim1,2)/h(nrho-1)+u(idim1,2)
138  vv(2,idim1,idim1) = - v(idim1,2)/h(nrho-1)
139  ENDDO
140 
141  doloop1irho1: DO irho1 = 2,nrho-1
142  irhop = irho1+1
143  irhom = irho1-1
144  h12 = 0.5*(h(irho1)+h(irhom))
145 
146  doloop2idim1 : DO idim1 = 1,ndim
147  ff(irho1,idim1) = b(idim1,irho1)*ym(idim1,irho1)+f(idim1,irho1)*tau
148  aa(irho1,idim1,idim1) = -0.5/c(idim1,irho1)/h(irho1)/ h12 &
149  * (d(idim1,irho1)+ d(idim1,irhop))
150  cc(irho1,idim1,idim1) = -0.5/c(idim1,irho1)/h(irhom)/ h12 &
151  * (d(idim1,irho1)+ d(idim1,irhom))
152  bb(irho1,idim1,idim1) = -aa(irho1,idim1,idim1) -cc(irho1,idim1,idim1)
153 
154  IF (iupwind.NE.1) THEN
155  aa(irho1,idim1,idim1) = aa(irho1,idim1,idim1) &
156  + 0.25 /c(idim1,irho1)/h12 &
157  * (e(idim1,irho1)+e(idim1,irhop))
158  cc(irho1,idim1,idim1) = cc(irho1,idim1,idim1) &
159  - 0.25 /c(idim1,irho1)/h12 &
160  * (e(idim1,irho1)+e(idim1,irhom))
161  bb(irho1,idim1,idim1) = bb(irho1,idim1,idim1) &
162  + 0.25 /c(idim1,irho1)/h12 &
163  * (e(idim1,irhop)-e(idim1,irhom))
164  ELSE
165  IF (e(idim1,irho1).GE.0.) THEN
166  cc(irho1,idim1,idim1) = cc(irho1,idim1,idim1) &
167  - e(idim1,irhom)/c(idim1,irho1)/h12
168  bb(irho1,idim1,idim1) = bb(irho1,idim1,idim1) &
169  + e(idim1,irho1) / c(idim1,irho1)/h12
170  ELSE
171  aa(irho1,idim1,idim1) = aa(irho1,idim1,idim1) &
172  + e(idim1,irhop)/c(idim1,irhop)/h12
173  bb(irho1,idim1,idim1) = bb(irho1,idim1,idim1) &
174  + e(idim1,irho1) / c(idim1,irho1)/h12
175  ENDIF
176 
177  ENDIF
178 
179  bb(irho1,idim1,idim1) = bb(irho1,idim1,idim1)+g(idim1,irho1)
180 
181  DO idim2 = 1,ndim
182  bb(irho1,idim1,idim2) = bb(irho1,idim1,idim2)!- CM1(IDIM1,IDIM2,IRHO1)
183  ENDDO
184 
185 
186  ENDDO doloop2idim1
187 
188 
189  ENDDO doloop1irho1
190 
191 
192 
193 !------------------------------COUPLING TERMS ----------------------------
194 ! DO IRHO1 = 2,NRHO-1
195 ! IRHOP = IRHO1+1
196 ! IRHOM = IRHO1-1
197 ! DO IDIM1 = 1,NDIM
198 ! H12 = 0.5*(H(IRHO1)+H(IRHOM))
199 ! if(IUPWIND.NE.1) THEN
200 ! AA(IRHO1,NDIM,IDIM1) = AA(IRHO1,NDIM,IDIM1)+ &
201 ! 0.25 *S2(IRHO1)/H12*(S3(IDIM1,IRHO1)+S3(IDIM1,IRHOP))
202 ! CC(IRHO1,NDIM,IDIM1) = CC(IRHO1,NDIM,IDIM1)- &
203 ! -0.25 *S2(IRHO1)/H12*(S3(IDIM1,IRHO1)+S3(IDIM1,IRHOM))
204 ! BB(IRHO1,NDIM,IDIM1) = BB(IRHO1,NDIM,IDIM1) +&
205 ! 0.25 /S2(IRHO1)/H12*(S3(IDIM1,IRHOP)-S3(IDIM1,IRHOM))
206 ! ELSE
207 ! IF(E(IDIM1,IRHO1).GE.0.) THEN
208 ! CC(IRHO1,NDIM,IDIM1) = CC(IRHO1,NDIM,IDIM1) - &
209 ! S2(IRHOM)*S3(IDIM1,IRHOM)/H12
210 ! BB(IRHO1,NDIM,IDIM1) = BB(IRHO1,NDIM,IDIM1)+ &
211 ! S2(IRHO1) * S3(IDIM1,IRHO1)/H12
212 ! ELSE
213 ! AA(IRHO1,NDIM,IDIM1) = AA(IRHO1,NDIM,IDIM1)+ &
214 ! S2(IRHOP)*S3(IDIM1,IRHOP)/H12
215 ! BB(IRHO1,NDIM,IDIM1) = BB(IRHO1,NDIM,IDIM1)+ &
216 ! S2(IRHO1) * S3(IDIM1,IRHO1)/H12
217 ! ENDIF
218 ! ENDIF
219 
220 
221 
222  bb = bb*tau
223  aa = aa*tau
224  cc = cc*tau
225 
226 
227  DO irho1 = 1,nrho
228  DO idim1 = 1,ndim
229  bb(irho1,idim1,idim1) = bb(irho1,idim1,idim1)+ a(idim1,irho1)
230  ENDDO
231  ENDDO
232 
233 
234 
235 ! +++ Call the matrix progonka:
236  CALL mprogonka(nrho,ndim,aa,bb,cc,ff,uu,vv,ww,yy)
237 
238 
239 ! +++ New solution:
240  DO idim1 = 1,ndim
241  DO irho1 = 1,nrho
242  y(idim1,irho1) = yy(irho1,idim1)*amix+y(idim1,irho1)*(1.e0_r8-amix)
243  ENDDO
244  ENDDO
245 
246  CALL derivn3(nrho,ndim,rho,y,dy)
247 
248 
249 ! +++ Return solution to the work flow:
250  DO idim1 = 1,ndim
251  DO irho1 = 1,nrho
252  solver%Y(idim1,irho1) = y(idim1,irho1)
253  solver%DY(idim1,irho1) = dy(idim1,irho1)
254  ENDDO
255  ENDDO
256 
257 
258 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
259 ! +++ Return solution to ETS:
260 
261 
262  RETURN
263 
264 
265 
266 END SUBROUTINE solution3
267 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
268 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
269 
270 
271 
272 
273 
274 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
286 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
287 SUBROUTINE mprogonka (NP,NDIM,A,B,C,F,U,V,W,Y)
288 ! Finds Solution of discrete equation of the form
289 ! A_i * Y_i+1 + B_i * Y_i +C_i * Y_i-1 = F_i
290 ! where:
291 ! Y_i - unknown wektor at point i (Y_i(NDIM))
292 ! F_i - Rhs at point i (F_i(NDIM))
293 ! A_i,B_i,C_i - equation coefficients (matrixes (NDIM,NDIM))
294 ! U,V,W - matrixes with boundary conditions
295 
296 
297  use itm_types
298  IMPLICIT NONE
299 !! Number of points
300  INTEGER, INTENT(IN) :: np
301 !! Rank of the block matrix
302  INTEGER, INTENT(IN) :: ndim
303 !! Equation coefficients: A, B, C
304  REAL(R8), DIMENSION(NP,NDIM,NDIM), INTENT(IN) :: a, b, c
305 !! Right hand side: F
306  REAL(R8), DIMENSION(NP,NDIM), INTENT(IN) :: f
307 !! Boundary conditions arrays U,V
308  REAL(R8), DIMENSION(2,NDIM,NDIM), INTENT(IN) :: u,v
309 !! Boundary condition array W
310  REAL(R8), DIMENSION(2,NDIM), INTENT(IN) :: w
311 !! Solution: Y
312  REAL(R8), DIMENSION(NP,NDIM), INTENT(OUT) :: y
313 
314 ! working variables
315 !! Progonka coefficients: alpha, beta
316  REAL(R8), DIMENSION(NP,NDIM,NDIM) :: alf
317  REAL(R8), DIMENSION(NP,NDIM) :: bet
318  REAL(R8), DIMENSION(NDIM,NDIM) :: temp1,temp2,temp3
319  REAL(R8), DIMENSION(NP,NDIM) :: wekt1
320  REAL(R8), DIMENSION(NDIM) :: wekt2
321  INTEGER :: i, j, j2, k, n, nn, n2
322  REAL(R8) :: stemp1,stemp2,det
323 
324 
325 
326 
327  IF(ndim .EQ. 1) THEN
328 ! BOUNDARY CONDITION
329  alf(2,1,1) =-v(1,1,1)/u(1,1,1)
330  bet(2,1)= w(1,1)/u(1,1,1)
331 
332  DO j=2,np-1
333  j2=j+1
334  det=b(j,1,1)+alf(j,1,1)*c(j,1,1)
335  IF(abs(det).GT.1.e-30_r8) THEN
336  alf(j2,1,1)=-a(j,1,1)/det
337  bet(j2,1)=(f(j,1)-bet(j,1)*c(j,1,1))/det
338  ELSE
339  alf(j2,1,1)=0.
340  bet(j2,1)=0.
341  ENDIF
342  ENDDO
343 
344 ! BOUNDARY CONDITION
345  y(np,1)=0.
346  det=v(2,1,1)*alf(np,1,1)+u(2,1,1)
347  IF(abs(det).GT.1.e-30_r8) y(np,1)=(w(2,1)-bet(np,1)*v(2,1,1))/det
348 
349  DO j=np-1,1,-1
350  j2=j+1
351  y(j,1)=alf(j2,1,1)*y(j2,1)+bet(j2,1)
352  ENDDO
353 
354 
355  ELSE
356 ! BOUNDARY CONDITION
357  DO i=1,ndim
358  DO j=1,ndim
359  temp1(i,j)=u(1,i,j)
360  ENDDO
361  ENDDO
362  CALL invmatrix(ndim, temp1, temp2)
363 
364  DO i=1,ndim
365  bet(2,i)=0.0
366  DO j=1,ndim
367  alf(2,i,j)=0.0
368  DO k=1,ndim
369  alf(2,i,j)=alf(2,i,j)-temp2(i,k)*v(1,k,j)
370  ENDDO
371  bet(2,i)=bet(2,i)+temp2(i,j)*w(1,j)
372  ENDDO
373  ENDDO
374 
375  DO n=2,np-1
376  DO i=1,ndim
377  wekt1(n,i)=f(n,i)
378  DO j=1,ndim
379  temp1(i,j)=b(n,i,j)
380  DO k=1,ndim
381  temp1(i,j)=temp1(i,j)+c(n,i,k)*alf(n,k,j)
382  ENDDO
383  wekt1(n,i)=wekt1(n,i)-c(n,i,j)*bet(n,j)
384  ENDDO
385  ENDDO
386  CALL invmatrix(ndim, temp1, temp2)
387 
388  n2=n+1
389  DO i=1,ndim
390  bet(n2,i)=0.
391  DO j=1,ndim
392  alf(n2,i,j)=0.
393  DO k=1,ndim
394  alf(n2,i,j)=alf(n2,i,j)-temp2(i,k)*a(n,k,j)
395  ENDDO
396  bet(n2,i)=bet(n2,i)+temp2(i,j)*wekt1(n,j)
397  ENDDO
398  ENDDO
399  ENDDO
400 
401 
402 1001 FORMAT(1x,3i4,1p,8e12.4)
403 
404 
405 ! BOUNDARY CONDITION
406  DO i=1,ndim
407  wekt2(i)=w(2,i)
408  DO j=1,ndim
409  temp1(i,j)=u(2,i,j)
410  DO k=1,ndim
411  temp1(i,j)=temp1(i,j)+v(2,i,k)*alf(np,k,j)
412  ENDDO
413  wekt2(i)=wekt2(i)-v(2,i,j)*bet(np,j)
414  ENDDO
415  ENDDO
416  CALL invmatrix(ndim, temp1, temp2)
417 
418  DO i=1,ndim
419  y(np,i)=0.
420  DO j=1,ndim
421  y(np,i)=y(np,i)+temp2(i,j)*wekt2(j)
422  ENDDO
423  ENDDO
424 
425  DO n=np-1,1,-1
426  n2=n+1
427  DO i=1,ndim
428  y(n,i)=bet(n2,i)
429  DO j=1,ndim
430  y(n,i)=y(n,i)+alf(n2,i,j)*y(n2,j)
431  ENDDO
432  ENDDO
433  ENDDO
434 
435 
436  ENDIF
437 
438 
439  RETURN
440 
441 END SUBROUTINE mprogonka
442 
443 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
448 ! Outputs
454 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
455 
456 SUBROUTINE invmatrix(n, AA, Ainv)
457 
458 ! Computes inverse of matrix
459 ! Input
460 !-------------------------------------------------
461 ! AA - Matrix A (n by n)
462 ! N - Dimension of matrix AA
463 !-------------------------------------------------
464 ! Outputs
465 ! Ainv - Inverse of matrix AA (n by n)
466 !-------------------------------------------------
467 
468  use itm_types
469 
470  IMPLICIT NONE
471 
472 !! order of matrix
473  INTEGER, INTENT(IN) :: n
474 !! Input Matrix
475  REAL(R8), DIMENSION(n,n), INTENT(IN) :: aa
476 
477 !! Output Inverted Matrix
478  REAL(R8), DIMENSION(n,n), INTENT(OUT) :: ainv
479 
480 
481  REAL(R8), DIMENSION(n) :: scale ! Scale factor
482  REAL(R8), DIMENSION(n,n) :: b ! Work array
483  REAL(R8), DIMENSION(n,n) :: a ! Working copy of input matrix
484  INTEGER, DIMENSION(n) :: index ! Index Vector
485  INTEGER :: i, j, k, jpivot, indexj
486  REAL(R8) :: scalemax, ratio, ratiomax, coeff, determ, sum
487 
488 
489 ! Copy matrix A so as not to modify original
490  do i=1,n
491  do j=1,n
492  a(i,j) = aa(i,j)
493  enddo
494  enddo
495 
496 !* Matrix b is initialized to the identity matrix
497  do i=1,n
498  do j=1,n
499  if( i .eq. j ) then
500  b(i,j) = 1.0
501  else
502  b(i,j) = 0.0
503  endif
504  enddo
505  enddo
506 
507 !* Set scale factor, scale(i) = max( |a(i,j)| ), for each row
508  do i=1,n
509  index(i) = i ! Initialize row index list
510  scalemax = 0.0
511  do j=1,n
512  if( abs(a(i,j)) .gt. scalemax ) then
513  scalemax = abs(a(i,j))
514  endif
515  enddo
516  scale(i) = scalemax
517  enddo
518 
519 !* Loop over rows k = 1, ..., (N-1)
520  do k=1,n-1
521 !* Select pivot row from max( |a(j,k)/s(j)| )
522  ratiomax = 0.0
523  jpivot = k
524  do i=k,n
525  ratio = abs(a(index(i),k))/scale(index(i))
526  if( ratio .gt. ratiomax ) then
527  jpivot=i
528  ratiomax = ratio
529  endif
530  enddo
531 !* Perform pivoting using row index list
532  indexj = index(k)
533  if( jpivot .ne. k ) then ! Pivot
534  indexj = index(jpivot)
535  index(jpivot) = index(k) ! Swap index jPivot and k
536  index(k) = indexj
537  endif
538 !* Perform forward elimination
539  do i=k+1,n
540  coeff = a(index(i),k)/a(indexj,k)
541  do j=k+1,n
542  a(index(i),j) = a(index(i),j) - coeff*a(indexj,j)
543  enddo
544  a(index(i),k) = coeff
545  do j=1,n
546  b(index(i),j) = b(index(i),j) - a(index(i),k)*b(indexj,j)
547  enddo
548  enddo
549  enddo
550 
551 
552 !* Perform backsubstitution
553  do k=1,n
554  ainv(n,k) = b(index(n),k)/a(index(n),n)
555  do i=n-1,1,-1
556  sum = b(index(i),k)
557  do j=i+1,n
558  sum = sum - a(index(i),j)*ainv(j,k)
559  enddo
560  ainv(i,k) = sum/a(index(i),i)
561  enddo
562  enddo
563 
564  return
565 
566 END SUBROUTINE invmatrix
567 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
568 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
569 
570 
571 
572 
573 
574 
575 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
582 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
583 SUBROUTINE derivn3(N,NFUN,X,FUN,DFUN)
584 
585  use itm_types
586 
587  IMPLICIT NONE
588 
589  INTEGER :: n ! number of radial points (input)
590  INTEGER :: nfun ! number of functions
591  INTEGER :: i, ifun
592 
593  REAL (R8) :: x(n), & ! argument array (input)
594  fun(nfun,n), & ! function array (input)
595  dfun(nfun,n), & ! function derivative array (output)
596  y(n), & ! function array 1-D (local)
597  dy1(n) ! function derivative array 1-D (local)
598  REAL (R8) :: h(n),dy2(n)
599 
600 
601 
602  DO i=1,n-1
603  h(i) = x(i+1)-x(i)
604  END DO
605 
606  DO ifun = 1,nfun
607 
608  DO i = 1,n
609  y(i) = fun(ifun,i)
610  END DO
611 
612  DO i=2,n-1
613  dy1(i) = ((y(i+1)-y(i))*h(i-1)/h(i)+(y(i)-y(i-1))*h(i)/h(i-1)) &
614  /(h(i)+h(i-1))
615  dy2(i) = 2.e0_r8*((y(i-1)-y(i))/h(i-1)+(y(i+1)-y(i))/h(i)) &
616  /(h(i)+h(i-1))
617  END DO
618 
619  dy1(1) = dy1(2)-dy2(2)*h(1)
620  dy1(n) = dy1(n-1)+dy2(n-1)*h(n-1)
621 
622 
623  DO i = 1,n
624  dfun(ifun,i) = dy1(i)
625  END DO
626 
627 
628  END DO
629 
630  RETURN
631 
632 
633 END SUBROUTINE derivn3
634 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
635 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
solver
Definition: solver.f90:1
type_solver
Definition: type_solver.f90:1
fun
subroutine fun(X, F)
Definition: Ev2.f:10
solution3
subroutine solution3(SOLVER, ifail)
This subroutine solves transport equations using matrix PROGONKA method.
Definition: solution3.f90:8
mprogonka
subroutine mprogonka(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: solution3.f90:287
derivn3
subroutine derivn3(N, NFUN, X, FUN, DFUN)
These subroutines calculate first and second derivatives, DY1 and DY2, of function Y respect to argum...
Definition: solution3.f90:583
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
invmatrix
subroutine invmatrix(n, AA, Ainv)
Computes inverse of matrix Input AA - Matrix A (n by n) N - Dimension of matrix AA.
Definition: solution3.f90:456