ETS  \$Id: Doxyfile 2162 2020-02-26 14:16:09Z g2dpc $
 All Classes Files Functions Variables Pages
solution10.F90
Go to the documentation of this file.
1 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
7 ! + + + + + + + + + NUMERICAL SOLUTION + + + + + + + + +
8 
9 SUBROUTINE solution10 (SOLVER, ifail)
10 
11 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
12 ! +
13 ! This subroutine is prepared to solve single +
14 ! transport equation in standardised form +
15 ! adopted by the ETS. +
16 ! +
17 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
18 ! Source: 1-D transport code RITM +
19 ! Developers: M.Basiuk, M Huynh +
20 ! +
21 ! Comments: based on CRONOS +
22 ! adapted for ETS +
23 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
24 
25 
26  USE type_solver
27 
28  use itm_types
29 
30  IMPLICIT NONE
31 
32  INTEGER, INTENT (INOUT) :: ifail
33 
34 ! +++ Input/Output with ETS:
35  TYPE (numerics), INTENT (INOUT) :: solver !contains all I/O quantities to numerics part
36 
37 #ifdef WANTCOS
38 
39 ! +++ Internal input/output parameters:
40  INTEGER :: idim, ndim !equation index and total number of equations to be solved
41 
42  INTEGER :: irho, nrho !radius index, number of radial points
43 
44  INTEGER :: flag !flag for equation: 0 - interpretative (not solved), 1 - predictive (solved)
45 
46  REAL (R8) :: rho(solver%nrho) !radii
47  REAL (R8) :: rhomax, rhomax2 !AF 9.Jul.2010
48 
49  REAL (R8) :: amix !fraction of new sollution mixed
50 
51  REAL (R8) :: y(solver%nrho), ym(solver%nrho) !function at the current amd previous time steps
52  REAL (R8) :: dy(solver%nrho) !derivative of function
53  REAL (R8) :: yy(solver%nrho) !new function !AF 12.Jul.2010
54  REAL (R8) :: dyy(solver%nrho) !derivative of new function !AF 12.Jul.2010
55 
56  REAL (R8) :: a(solver%nrho), b(solver%nrho) !coefficients
57  REAL (R8) :: c(solver%nrho), d(solver%nrho) !coefficients
58  REAL (R8) :: e(solver%nrho), f(solver%nrho) !coefficients
59  REAL (R8) :: g(solver%nrho), h !coefficients
60  REAL (R8) :: dd(solver%nrho), de(solver%nrho) !derivatives of coefficients
61 
62  REAL (R8) :: v(2), u(2), w(2) !boundary conditions
63 
64 
65 ! +++ Coefficients used by internal numerical solver:
66  INTEGER :: i, n ,ipsi !radius index, number of radial points
67  INTEGER :: isp !spline flag
68 
69  REAL (R8) :: x(solver%nrho) !radii
70 
71  REAL (R8) :: sol(solver%nrho) !solution
72  REAL (R8) :: dsol(solver%nrho) !derivative of solution
73 
74  REAL (R8) :: as(solver%nrho), bs(solver%nrho) !coefficients
75  REAL (R8) :: cs(solver%nrho) !coefficients
76 
77  REAL (R8) :: acron(solver%nrho), bcron(solver%nrho) !coefficients
78  REAL (R8) :: ccron(solver%nrho), dcron(solver%nrho) !coefficients
79 
80 ! REAL (R8) :: mata_in(SOLVER%NRHO,5,5), matb_in(SOLVER%NRHO,5,5) !coefficients ! AF, 10.Jul.2010
81 ! REAL (R8) :: matc_in(SOLVER%NRHO,5,5), matd_in(SOLVER%NRHO,5) !coefficients ! AF, 10.Jul.2010
82  REAL (R8) :: mata_in(solver%nrho,1,1), matb_in(solver%nrho,1,1) !coefficients ! AF, 10.Jul.2010
83  REAL (R8) :: matc_in(solver%nrho,1,1), matd_in(solver%nrho,1) !coefficients ! AF, 10.Jul.2010
84 
85  ! REAL (R8) :: VS(2), US(2), WS(2) !boundary conditions ! AF, 10.Jul.2010
86 
87  integer, parameter :: nbeq = 1 !AF 9.Jul.2010
88  REAL(R8) :: dt_in !AF 9.Jul.2010
89  REAL(R8) :: dx_in !AF 9.Jul.2010
90  REAL(R8) :: sca_f_in !AF 9.Jul.2010
91  REAL(R8), dimension(2,nbeq,3) :: t0_in,t1_in !AF 9.Jul.2010
92  REAL(R8), dimension(2,nbeq) :: v0_in,v1_in !AF 9.Jul.2010
93 
94  REAL(R8), dimension(nbeq) :: mode_in
95 
96 ! + + + + + + + + + INTERFACE PART + + + + + + + + + + +
97 
98 
99 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
100 ! +++ Set up local variables (input, obtained from ETS):
101 ! Control parameters:
102  ndim = solver%NDIM
103  nrho = solver%NRHO
104  amix = solver%AMIX
105  h = solver%H !AF, 10.Jul.2010
106 
107 
108 
109 ! +++ Solution of equations starting from 1 to NDIM:
110  equation_loop1: DO idim = 1, ndim
111 
112  flag = solver%EQ_FLAG(idim)
113 
114 
115  IF (flag.EQ.0) goto 20 !equation is not solved
116 
117 
118 
119 ! +++ Numerical coefficients obtained from physics part in form:
120 !
121 ! (A*Y-B*Y(t-1))/H + 1/C * (-D*Y' + E*Y) = F - G*Y ! AF 10.Jul.2010 - prime missing
122 ! (A*Y-B*Y(t-1))/H + 1/C * (-D*Y' + E*Y)' = F - G*Y ! AF 10.Jul.2010
123 ! conversion to solver COS
124 ! Ajj, Bjj, Cjj, Djj
125 !
126 ! dY/dt = Ajj Y'' + Bjj Y' + Cjj Y + Djj
127 !
128 
129  rho_loop1: DO irho=1,nrho
130  rho(irho) = solver%RHO(irho)
131 
132  y(irho) = solver%Y(idim,irho)
133  dy(irho) = solver%DY(idim,irho)
134  ym(irho) = solver%YM(idim,irho)
135 
136  a(irho) = solver%A(idim,irho)
137  b(irho) = solver%B(idim,irho)
138  c(irho) = solver%C(idim,irho)
139  d(irho) = solver%D(idim,irho)
140  e(irho) = solver%E(idim,irho)
141  f(irho) = solver%F(idim,irho)
142  g(irho) = solver%G(idim,irho)
143 
144  END DO rho_loop1
145 
146 
147  rhomax = rho(nrho)
148  rhomax2 = rhomax*rhomax
149  !pi2 = pi*pi ! AF 9.Jul.2010 - these are not being used
150  CALL derivn10(nrho,rho,d,dd)
151  CALL derivn10(nrho,rho,e,de)
152 
153  rho_loop2: DO irho=1,nrho
154 
155  acron(irho) = d(irho) / c(irho) / rhomax2
156  bcron(irho) = (dd(irho) / c(irho) - e(irho) / c(irho)) / rhomax
157 
158  END DO rho_loop2
159 
160 
161  rho_loop3: DO irho=1,nrho
162 
163 ! Ccron(IRHO) = -DE(IRHO) / C(IRHO) + G(IRHO) !AF, 10.Jul.2010
164  ccron(irho) = -de(irho) / c(irho) - g(irho) !AF, 10.Jul.2010
165  dcron(irho) = f(irho)
166 
167  END DO rho_loop3
168 
169  ccron = ccron + (b-a)/h !AF, 10.Jul.2010
170 
171  acron = acron/a !AF, 10.Jul.2010
172  bcron = bcron/a !AF, 10.Jul.2010
173  ccron = ccron/a !AF, 10.Jul.2010
174  dcron = dcron/a !AF, 10.Jul.2010
175 
176 ! case one equation
177  ipsi = 1
178 ! case multi equation (not avialable for the ETS
179 ! ipsi = 1 : current diffusion equation
180 ! ipsi = 2 : Heat Te
181 ! ipsi = 3 : Heat Ti
182 ! ipsi = 4 : ne density
183 ! ipsi = 5 : Vtor
184 !
185  call cos_zconversion(acron,nrho)
186 ! mata_in(1:nrho,ipsi,ipsi)=Acron(1:nbrho) !AF 9.Jul.2010
187  mata_in(1:nrho,ipsi,ipsi)=acron(1:nrho) !AF 9.Jul.2010
188 
189  call cos_zconversion(bcron,nrho)
190 ! matb_in(1:nbrho,ipsi,ipsi)=Bcron(1:nrho) !AF 9.Jul.2010
191  matb_in(1:nrho,ipsi,ipsi)=bcron(1:nrho) !AF 9.Jul.2010
192 
193  call cos_zconversion(ccron,nrho)
194 ! matc_in(1:nbrho,ipsi,ipsi) = Ccron(1:nrho) !AF 9.Jul.2010
195  matc_in(1:nrho,ipsi,ipsi) = ccron(1:nrho) !AF 9.Jul.2010
196 
197  call cos_zconversion(dcron,nrho)
198 ! matd_in(:,ipsi) = Dcron(1:nrho) !AF 9.Jul.2010
199  matd_in(1:nrho,ipsi) = dcron(1:nrho) !AF 9.Jul.2010
200 
201 
202 
203 ! boundary conditions - AF 10.Jul.2010
204 
205 u = solver%U(idim,:)
206 v = solver%V(idim,:)
207 w = solver%W(idim,:)
208 
209 ! rho = 0
210 t0_in(:,1,3) = 0.e0_r8
211 t0_in(:,1,2) = v(1)
212 t0_in(:,1,1) = u(1)
213 v0_in(:,1) = w(1)
214 ! rho = 1
215 t1_in(:,1,3) = 0.e0_r8
216 t1_in(:,1,2) = v(2)
217 t1_in(:,1,1) = u(2)
218 v1_in(:,1) = w(2)
219 
220 ! AF 10.Jul.2010: 1 for t, 2 for t+dt
221 ! do i=1,2
222 
223  ! poloidal magnetic flux
224 ! T0_in(i,1,3) = 0
225 ! T0_in(i,1,2) = 1
226 ! T0_in(i,1,1) = 0
227 ! V0_in(i,1) = W(0)
228 
229 ! T1_in(i,1,3) = 0
230 ! T1_in(i,1,2) = 1
231 ! T1_in(i,1,1) = 0
232 ! V1_in(i,1) = W(2)
233 
234 ! end do
235 
236 ! AF - end
237 
238 
239 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
240 ! +++ Solution of transport equation:
241  !nbeq = 1 !AF - 9.Jul.2010 - moved above to be used in the dimensions of the boundary parameters T0_in, etc.
242  dt_in = h !0.01 !AF - 10.Jul.2010
243  dx_in = 1.e0_r8/(nrho-1) !0.01 !AF - 10.Jul.2010
244 ! 0 implicite, 1 explicite 0.5 Krank Nicholson
245  sca_f_in = 0.e0_r8 !AF - 9.Jul.2010
246 
247  mode_in = 1-flag !AF - 12.Jul.2010
248 
249  call cos_pde1solver(ym,yy, & !AF - 12.Jul.2010 - to have the mixing option
250  !call cos_pde1solver(YM,Y, & !AF - 10.Jul.2010, 12.Jul.2010
251  !call cos_pde1solver(Y,DY, & !AF - 10.Jul.2010
252  mata_in,matb_in,matc_in,matd_in, &
253  mata_in,matb_in,matc_in,matd_in, &
254  nrho,nbeq,nbeq,dx_in,dt_in,sca_f_in, &
255  t0_in,t1_in,v0_in,v1_in,mode_in) ! AF - 12.Jul.2010
256 
257  y = y*(1.e0_r8-amix) + yy*amix ! AF - 12.Jul.2010
258  CALL derivn10(nrho,rho,y,dy) ! AF - 12.Jul.2010
259  solver%Y(idim,:) = y
260  solver%DY(idim,:) = dy
261 
262  ! RHO_LOOP4: DO IRHO = 1,NRHO ! AF - 12.Jul.2010
263 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
264 ! +++ New function and derivative:
265  ! Y(IRHO) = Y(IRHO)*(1.e0_R8-AMIX) + SOL(IRHO)*AMIX !AF - 9.Jul.2010
266  ! DY(IRHO) = DY(IRHO)*(1.e0_R8-AMIX) + DSOL(IRHO)*AMIX !AF - 9.Jul.2010
267 
268 
269 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
270 ! +++ Return solution to ETS:
271  !SOLVER%Y(IDIM,IRHO) = Y(IRHO) ! AF - 12.Jul.2010
272  !SOLVER%DY(IDIM,IRHO) = DY(IRHO)! AF - 12.Jul.2010
273 
274  !END DO RHO_LOOP4 ! AF - 12.Jul.2010
275 
276 
277 20 CONTINUE
278 
279 
280  END DO equation_loop1
281 
282 #endif
283 
284  RETURN
285 
286 
287 
288 END SUBROUTINE solution10
289 
290 #ifdef WANTCOS
291 
292 
293 
294 ! + + + + + + + + + NUMERICAL PART + + + + + + + + + + +
295 
296 
297 
298 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
305 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
306 SUBROUTINE derivn10(N,X,Y,DY1) !AF 10.Jul.2010 - this is actually the same as DERIVN1 in solver1 - should this clone remain here for independence from solver 1?
307 
308  use itm_types
309 
310  IMPLICIT NONE
311 
312  INTEGER :: n ! number of radial points (input)
313  INTEGER :: i
314 
315  REAL (R8) :: x(n), & ! argument array (input)
316  y(n), & ! function array (input)
317  dy1(n) ! function derivative array (output)
318  REAL (R8) :: h(n),dy2(n)
319 
320  DO i=1,n-1
321  h(i)=x(i+1)-x(i)
322  END DO
323 
324  DO i=2,n-1
325  dy1(i)=((y(i+1)-y(i))*h(i-1)/h(i)+(y(i)-y(i-1))*h(i)/h(i-1)) &
326  /(h(i)+h(i-1))
327  dy2(i)=2.e0_r8*((y(i-1)-y(i))/h(i-1)+(y(i+1)-y(i))/h(i)) &
328  /(h(i)+h(i-1))
329  END DO
330 
331  dy1(1)=dy1(2)-dy2(2)*h(1)
332  dy1(n)=dy1(n-1)+dy2(n-1)*h(n-1)
333 
334  RETURN
335 END SUBROUTINE derivn10
336 
337 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
339 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
340 
341 subroutine cos_zconversion(e,nbrho)
342  use cos_precision
343  implicit none
344 
345  integer,intent(in) :: nbrho
346  real(DP),dimension(nbrho) :: e
347 
348  if (nbrho < 4) then
349  print*, "error"
350  return
351  endif
352 
353  ! continuite au centre pour les NaN
354  if ( (abs(e(1)).gt.huge(e(1))) .or. (abs(e(1)).eq.0) ) then
355  e(1) = (61.0_8/46.0_8) * e(2) - (9.0_8/23.0_8) * e(3) + (3.0_8/46.0_8) * e(4)
356  endif
357 
358 end subroutine cos_zconversion
359 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
360 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
361 
362 
363 !********************************************************************************
364 module cos_transportdata
365  implicit none
366  integer :: dimk,dimm,dimmat,dimbb,dimn
367  real(kind=8),dimension(:,:), pointer :: f,fp,dfpdx,dfpdx2
368  real(kind=8),dimension(:,:),allocatable :: matd,matdp
369  real(kind=8),dimension(:,:,:),pointer :: t0,t1
370  real(kind=8),dimension(:,:,:),allocatable :: mata,matb,matc, &
371  matap,matbp,matcp
372  real(kind=8),dimension(:,:), pointer :: v0,v1
373  integer,dimension(:),allocatable :: mode
374 
375  real(kind=8),dimension(:,:),allocatable :: mats
376  real(kind=8),dimension(:,:),allocatable :: matu,matv,matw,matx,maty,matz
377  real(kind=8),dimension(:),allocatable,target:: mat_coo,mat_csr
378  real(kind=8),dimension(:),allocatable,target:: vectbb
379  integer,dimension(:),allocatable,target :: cooi,cooj,csri,csrj
380  real(kind=8), pointer:: dx,dt,sca_f
381  integer :: verbose = 0
382 end module cos_transportdata
383 !module cos_pde1dsolver_interface ! AF - 9.Jul.2010, this interface is not in use
384 ! interface
385 ! subroutine cos_pde1solver(f_in,fp_in, &
386 ! mata_in,matb_in,matc_in,matd_in, &
387 ! matap_in,matbp_in,matcp_in,matdp_in, &
388 ! dimk_in,dimn_in,dimm_in,dx_in,dt_in,sca_f_in, &
389 ! T0_in,T1_in,V0_in,V1_in,mode_in,dfpdx_in,dfpdx2_in)
390 ! use cos_transportdata
391 ! implicit none
392 ! integer, intent(in) :: dimk_in,dimn_in,dimm_in
393 ! real(kind=8),intent(in),target :: dt_in,dx_in,sca_f_in
394 ! real(kind=8),dimension(dimk_in,dimm_in),target :: f_in,fp_in
395 ! real(kind=8),dimension(dimk_in,dimm_in),target,optional:: dfpdx_in,dfpdx2_in
396 ! real(kind=8),dimension(dimk_in,dimm_in),target :: matd_in,matdp_in
397 ! real(kind=8),dimension(dimk_in,dimm_in,dimm_in),target :: mata_in,matb_in,matc_in, &
398 ! matap_in,matbp_in,matcp_in
399 ! real(kind=8),dimension(2,dimm_in,3),target :: T0_in,T1_in
400 ! real(kind=8),dimension(2,dimm_in),target :: V0_in,V1_in
401 ! real(kind=8),dimension(dimm_in) :: mode_in
402 ! end subroutine cos_pde1solver
403 ! end interface
404 !end module cos_pde1dsolver_interface
405 !*********************************************************************************
406 subroutine cos_pde1solver(f_in,fp_in, &
407  mata_in,matb_in,matc_in,matd_in, &
408  matap_in,matbp_in,matcp_in,matdp_in, &
409  dimk_in,dimn_in,dimm_in,dx_in,dt_in,sca_f_in, &
410  t0_in,t1_in,v0_in,v1_in,mode_in,dfpdx_in,dfpdx2_in)
411  !*********************************************************************
412  !solve the equations :
413  ! Diff(f,t) = mata * diff(f,x,x) + matb * diff(f,x) + matc * f + matd
414  !
415  !dt_in : time step
416  !dx_in : space step
417  !
418  !f : solution at the time t
419  !fp : solution at the time t+dt
420  !dfpdx : diff(fp,x)
421  !dfpdx2 : diff(fp,x,x)
422  !
423  !dimk_in : dimension of space index
424  ! x=x0 + (k-1)*dx, k->[1..dimk]
425  !dimm_in : dimension of equation number
426  !dimn_in : dimension of the solution at a one grid point
427  !
428  !sca_f_in : scalar
429  ! 0 -> full implicit
430  ! 0.5 -> Crank Nicholson
431  ! 1 -> explicit (don't used it)
432  !
433  !
434  !The big matrix
435  !--------------
436  !
437  ! k=1 k=2 ----------------------------------------- k=K
438  ! (1,1)--------(1,m)(1,1)--------(1,m)
439  ! | b(1,m,n) | | c(1,m,n) |
440  ! | | | |
441  ! (n,1)--------(n,m)(n,1)--------(n,m)
442  ! (1,1)--------(1,m)(1,1)--------(1,m)(1,1)--------(1,m)
443  ! | a(2,m,n) | | b(2,m,n) | | c(2,m,n) |
444  ! | | | | | |
445  ! (n,m)--------(n,m)(n,m)--------(n,m)(n,m)--------(n,m)
446  ! .
447  ! .
448  ! .
449  ! .
450  ! .
451  ! .
452  ! .
453  ! .
454  ! .
455  !
456  ! (1,1)--------(1,m)(1,1)--------(1,m)
457  ! | a(K,m,n) | | b(K,m,n) |
458  ! | | | |
459  ! (n,m)--------(n,m)(n,m)--------(n,m)
460  !
461  ! where n=m if we are in predictive mode
462  ! n<m if we are in interpretative mode
463  ! the first index is the space index
464  ! the second index is the coupling coefficient between variables
465  ! the third index is the rank of the equation
466  !
467  ! T0(l,m,j) : coefficients of the boundaries condition at x0
468  ! V0(l,m) : values of the boundaries condition at x0
469  ! T0(:,:,3) * diff(f,x,x) + T0(:,:,2) * diff(f,x) + T0(:,:,1) * f = V0
470  ! where l=1 for t
471  ! l=2 for t+dt
472  ! T1(l,m,j) , V1(l,m) : idem as above but at x1
473  !*************************************************************************************
474  use cos_transportdata
475  implicit none
476  integer, intent(in) :: dimk_in,dimn_in,dimm_in
477  real(kind=8),intent(in),target :: dt_in,dx_in,sca_f_in
478  real(kind=8),dimension(dimk_in,dimm_in),target :: f_in,fp_in
479  real(kind=8),dimension(dimk_in,dimm_in),target,optional:: dfpdx_in,dfpdx2_in
480  real(kind=8),dimension(dimk_in,dimm_in),target :: matd_in,matdp_in
481  real(kind=8),dimension(dimk_in,dimm_in,dimm_in),target :: mata_in,matb_in,matc_in, &
482  matap_in,matbp_in,matcp_in
483  real(kind=8),dimension(2,dimm_in,3),target :: t0_in,t1_in
484  real(kind=8),dimension(2,dimm_in),target :: v0_in,v1_in
485 
486  real(kind=8),dimension(:,:,:), allocatable :: mat_aux
487 
488  real(kind=8),dimension(dimm_in) :: mode_in
489  integer,dimension(:),allocatable :: iwk,cooi_aux
490 
491  real(kind=8) :: dtx2,dt2x,aux
492  integer :: i,j,k,m,n
493  integer :: counteri,counterj,counter
494  integer :: interpretatif=0,condlim=0
495 
496  interface
497  integer function condlimord2_k1(interpretatif)
498  integer,intent(in) :: interpretatif
499  end function condlimord2_k1
500 
501  integer function condlimord2_kk(interpretatif)
502  integer,intent(in) :: interpretatif
503  end function condlimord2_kk
504  end interface
505 
506 
507  write(*,*) 'dans solver',sca_f_in
508  if (sca_f_in==1.0) then
509  print*,"Error, the explicit Euler algorithm must not be used"
510  stop
511  endif
512 
513  !**********
514  !allocation
515  !**********
516  dimk=dimk_in
517  dimm=dimm_in
518  dimn=dimn_in
519  sca_f=>sca_f_in
520  dx=>dx_in
521  dt=>dt_in
522  f=>f_in
523  fp=>fp_in
524  if (present(dfpdx_in)) dfpdx=>dfpdx_in
525  if (present(dfpdx2_in)) dfpdx2=>dfpdx2_in
526  allocate(mata(dimk_in,dimm_in,dimm_in),matb(dimk_in,dimm_in,dimm_in),&
527  matc(dimk_in,dimm_in,dimm_in),matd(dimk_in,dimm_in) )
528  allocate(matap(dimk_in,dimm_in,dimm_in),matbp(dimk_in,dimm_in,dimm_in),&
529  matcp(dimk_in,dimm_in,dimm_in),matdp(dimk_in,dimm_in) )
530  allocate(mode(dimm_in))
531 
532  mata=mata_in
533  matb=matb_in
534  matc=matc_in
535  matd=matd_in
536 
537  matap=matap_in
538  matbp=matbp_in
539  matcp=matcp_in
540  matdp=matdp_in
541  !print*,"matap",matap(1:10,1,1)
542 
543  mode=int(mode_in)
544  t0=>t0_in
545  t1=>t1_in
546  v0=>v0_in
547  v1=>v1_in
548  allocate(mats(dimk,dimm))
549  allocate(matu(2,dimm),matv(2,dimm),matw(2,dimm))
550  allocate(matx(2,dimm),maty(2,dimm),matz(2,dimm))
551  allocate(mat_aux(dimk,dimm,dimm))
552 
553  dimmat=(dimk-2)*dimn*dimn*3 +4*dimn*dimn
554  dimbb=dimn*dimk
555  allocate(vectbb(dimbb),mat_coo(dimmat),cooi(dimmat),cooj(dimmat))
556 
557  !**************
558  !initialisation
559  !**************
560  interpretatif=sum(mode)
561  if (interpretatif.ne.0) interpretatif=1
562  dtx2 = dt/(dx*dx)
563  dt2x = dt/(2.0_8*dx)
564  write(*,*) 'dtx2=',dtx2
565  mat_aux = matb
566  matb = -2.0_8*dtx2*mata + dt*matc
567  matc = dtx2*mata + dt2x*mat_aux
568  mata = dtx2*mata - dt2x*mat_aux
569  matd = dt*matd
570 
571  mat_aux = matbp
572  matbp = -2.0_8*dtx2*matap + dt*matcp
573  matcp = dtx2*matap + dt2x*mat_aux
574  matap = dtx2*matap - dt2x*mat_aux
575  matdp = dt*matdp
576 
577  deallocate(mat_aux)
578 
579  mats = sca_f*matd + (1.0_8-sca_f)*matdp + f;
580 
581  !*********************
582  !boundaries conditions
583  !*********************
584  condlim=sum(t0(:,:,3))
585  if ( condlimord2_k1(interpretatif) > 0) then
586  if (condlim>0) then
587  print*," boundaries condition of order 2 is not implemented"
588  stop
589  else
590  !print*,"boundaries condition of order 1 k=1"
591  call condlimord1_k1
592  endif
593  endif
594 
595  condlim=sum(t1(:,:,3))
596  if ( condlimord2_kk(interpretatif) > 0) then
597  if (condlim>0) then
598  print*," boundaries condition of order 2 is not implemented"
599  stop
600  else
601  !print*,"boundaries condition of order 1 k=dimk"
602  call condlimord1_kk
603  endif
604  endif
605 
606 
607  !*******************
608  !interpretative mode
609  !*******************
610  write(*,*) 'dimm=',dimm
611  do m = 1,dimm
612  if (mode(m)==1) cycle
613  do k=2,dimk-1
614  do n=1,dimm
615  if (mode(n)==1) then
616  mats(k,m)=mats(k,m) + sca_f*(mata(k,n,m)*f(k-1,n)+matb(k,n,m)*f(k,n)+matc(k,n,m)*f(k+1,n)) &
617  + (1.0_8-sca_f)*(matap(k,n,m)*fp(k-1,n)+matbp(k,n,m)*fp(k,n)+matcp(k,n,m)*fp(k+1,n))
618  endif
619  enddo
620  enddo
621  enddo
622 
623  !*************************
624  !filling the second member
625  !*************************
626  !for k=1
627  k=1
628  counter=0
629  do m=1,dimm
630  if (mode(m)==0) then
631  counter=counter+1
632  vectbb(counter) = mats(k,m)
633  do n=1,dimm
634  if (mode(n)==0) then
635  vectbb(counter) = vectbb(counter) + sca_f*(matb(k,n,m)*f(k,n) + matc(k,n,m)*f(k+1,n))
636  endif
637  enddo
638  endif
639  enddo
640 
641  !for k =[2,dimk-1]
642  do k=2,dimk-1
643  do m=1,dimm
644  if (mode(m)==0) then
645  counter=counter+1
646  vectbb(counter) = mats(k,m)
647  do n=1,dimm
648  if (mode(n)==0) then
649  vectbb(counter) = vectbb(counter) + sca_f*(mata(k,n,m)*f(k-1,n) + matb(k,n,m)*f(k,n) &
650  + matc(k,n,m)*f(k+1,n))
651  endif
652  enddo
653  endif
654  enddo
655  enddo
656 
657  !for k=dimk
658  k=dimk
659  do m=1,dimm
660  if (mode(m)==0) then
661  counter=counter+1
662  vectbb(counter) = mats(k,m)
663  do n=1,dimm
664  if (mode(n)==0) then
665  vectbb(counter) = vectbb(counter) + sca_f*(mata(k,n,m)*f(k-1,n) + matb(k,n,m)*f(k,n))
666  endif
667  enddo
668  endif
669  enddo
670 
671  !***********************
672  !filling the big matrix
673  !***********************
674  k=1
675  counter =0
676  counteri=0
677  counterj=0
678  write(*,*) 'counter=',counter
679  call build_matcoo(matbp,1,k,counter,counteri,counterj)
680  counteri=0
681  counterj=dimn
682  call build_matcoo(matcp,0,k,counter,counteri,counterj)
683  do k=2,dimk-1
684  counteri=(k-1)*dimn
685  counterj=(k-2)*dimn
686  call build_matcoo(matap,0,k,counter,counteri,counterj)
687  counteri=(k-1)*dimn
688  counterj=(k-1)*dimn
689  call build_matcoo(matbp,1,k,counter,counteri,counterj)
690  counteri=(k-1)*dimn
691  counterj=k*dimn
692  call build_matcoo(matcp,0,k,counter,counteri,counterj)
693  enddo
694  k=dimk
695  counteri=(k-1)*dimn
696  counterj=(k-2)*dimn
697  call build_matcoo(matap,0,k,counter,counteri,counterj)
698  counteri=(k-1)*dimn
699  counterj=(k-1)*dimn
700  call build_matcoo(matbp,1,k,counter,counteri,counterj)
701  !*****************
702  ! print the matrix
703  !*****************
704  if (verbose>0) then
705  print*,"dimk",dimk,dimn,dimm,counter
706  allocate(iwk(dimbb+1),cooi_aux(dimmat))
707  allocate(mat_csr(dimmat),csri(dimmat),csrj(dimmat))
708  endif
709  dimmat=counter
710  if (verbose>0) then
711  print*,"SPARSE BIG MATRIX:"
712  cooi_aux=cooi
713  ! call coicsr (dimbb,dimmat,1,mat_coo,cooj,cooi,iwk)
714  !call coocsr(dimbb,dimmat,mat_coo,cooi_aux,cooj,mat_csr,csrj,csri)
715  ! call dump (1,dimbb,1,mat_coo,cooj,cooi,6)
716  !call dump (1,dimbb,1,mat_csr,csrj,csri,6)
717  print*,"RHS:"
718  print*,vectbb
719  endif
720  !****************
721  ! solve de matrix
722  !****************
723  call solvepde1d
724 
725  !*****************
726  ! save de solution
727  !*****************
728  counter=0
729  do k=1,dimk
730  do m=1,dimm
731  if (mode(m)==0) then
732  counter=counter+1
733  fp(k,m) = vectbb(counter)
734  endif
735  enddo
736  enddo
737 
738  !fp(1,1)=mat_coo(1)
739  !fp(2,1)=mat_coo(4)
740  !fp(3,1)=aux
741  !fp(4,1)=dt2x
742 
743  !do m=1,dimm
744  ! if (mode(m) == 1) then
745  ! do k=1,dimk
746  ! fp(k,m) = fp_in(k,m)
747  ! enddo
748  ! endif
749  !enddo
750  if (verbose > 0) then
751  print*,"SOLUTION:"
752  print*,vectbb
753  endif
754 
755  !************
756  !finalisation
757  !************
758  deallocate(matu,matv,matw,matx,maty,matz,mats)
759  deallocate(mata,matb,matc,matd,matap,matbp,matcp,matdp)
760  deallocate(mode)
761  if (verbose==0) then
762  deallocate(mat_coo,cooi,cooj,vectbb)
763  else
764  deallocate(mat_coo,cooi,cooj,vectbb,iwk,cooi_aux,mat_csr,csri,csrj)
765  endif
766 
767  !if (dfpdx(1,1)==0.and.dfpdx2(1,1)==0) then
768  !print*,"dfpdx and dfpdx2 not available"
769  !endif
770 
771 end subroutine cos_pde1solver
772 !**********************************************************************************
773 subroutine build_matcoo(matp,option,k,counter,counteri,counterj)
774  use cos_transportdata
775  implicit none
776  integer,intent(in) :: option,k
777  integer,intent(inout) :: counter,counteri,counterj
778  integer :: counteri_aux,counterj_aux
779  real(kind=8),dimension(dimk,dimm,dimm) :: matp
780 
781  integer :: m,n
782 
783  counteri_aux=counteri
784  counterj_aux=counterj
785  do m=1,dimm
786  !print*,"mod",mode(1)
787  if (mode(m)==0) then
788  counteri=counteri+1
789  counterj=counterj_aux
790  do n=1,dimm
791  if (mode(n)==0) then
792  counterj=counterj+1
793  counter=counter+1
794  if (counteri==counterj) then
795  mat_coo(counter)=1.0_8-(1.0_8-sca_f)*matp(k,n,m)
796  else
797  mat_coo(counter)=-(1.0_8-sca_f)*matp(k,n,m)
798  endif
799  if (mat_coo(counter)==0.0) then
800  counter=counter-1
801  !print*,"vide",counter,sca_f,matp(k,n,m)
802  else
803  cooi(counter)=counteri
804  cooj(counter)=counterj
805  endif
806  endif
807  enddo
808  endif
809  enddo
810  !print*,"counter",counter
811 end subroutine build_matcoo
812 !************************************************************************************
813 subroutine solvepde1d
814  use cos_transportdata
815  implicit none
816  include 'mpif.h'
817  include 'dmumps_struc.h'
818  type(dmumps_struc) :: mumps_par
819  integer:: ierr, i
820 
821  call mpi_init(ierr)
822  ! Define a communicator for the package.
823  mumps_par%COMM = mpi_comm_world
824  ! Initialize an instance of the package
825  ! for L U factorization (sym = 0, with working host)
826  mumps_par%JOB = -1
827  mumps_par%SYM = 0
828  mumps_par%PAR = 1
829  call dmumps(mumps_par)
830 
831  ! Define problem on the host (processor 0)
832  IF ( mumps_par%MYID .eq. 0 ) THEN
833  !order of the matrix
834  mumps_par%N = dimk*dimn
835  !number of entries
836  mumps_par%NZ = dimmat
837  !allocate( mumps_par%IRN ( mumps_par%NZ ) )
838  !allocate( mumps_par%JCN ( mumps_par%NZ ) )
839  !allocate( mumps_par%A( mumps_par%NZ ) )
840  !allocate( mumps_par%RHS ( mumps_par%N ) )
841  !DO I = 1, mumps_par%NZ
842  ! mumps_par%IRN(I) = cooi(i)
843  ! mumps_par%JCN(I) = cooj(i)
844  ! mumps_par%A(I) = mat_coo(i)
845  !END DO
846  !DO I = 1, mumps_par%N
847  ! mumps_par%RHS(I) = vectbb(i)
848  !END DO
849  mumps_par%IRN=>cooi
850  mumps_par%JCN=>cooj
851  mumps_par%A=>mat_coo
852  mumps_par%RHS=>vectbb
853  END IF
854  if (verbose <= 1) then
855  mumps_par%ICNTL(4)=1
856  mumps_par%ICNTL(3)=0
857  !only because we have a problem with mexfile
858  !mumps_par%ICNTL(1)=0
859  endif
860 
861  ! do i = 1,40
862  ! print*,"icntl",i,mumps_par%ICNTL(i)
863  ! enddo
864  ! Call package for solution
865  mumps_par%JOB = 6
866  CALL dmumps(mumps_par)
867  ! Solution has been assembled on the host
868  IF ( mumps_par%MYID .eq. 0 .and. verbose > 1) THEN
869  WRITE( 6, * ) ' MUMPS Solution is ',(mumps_par%RHS(i),i=1,mumps_par%N)
870  END IF
871  !vectbb=mumps_par%RHS
872  ! Deallocate user data
873  IF ( mumps_par%MYID .eq. 0 )THEN
874  !DEALLOCATE( mumps_par%IRN )
875  !DEALLOCATE( mumps_par%JCN )
876  !DEALLOCATE( mumps_par%A )
877  !DEALLOCATE( mumps_par%RHS )
878  END IF
879  ! Destroy the instance (deallocate internal data structures)
880  mumps_par%JOB = -2
881  CALL dmumps(mumps_par)
882  CALL mpi_finalize(ierr)
883 end subroutine solvepde1d
884 !**********************************************************************************
885 integer function condlimord2_k1(interpretatif)
886 
887  !***********************
888  !boundary conditions k=1
889  !***********************
890 
891  use cos_transportdata
892  implicit none
893 
894  integer,intent(in) :: interpretatif
895  integer :: dirichlet=0,neumann=0,m,i
896  real(kind=8) :: aux
897 
898  condlimord2_k1=0
899 
900  matu = t0(:,:,3)/(dx*dx) + t0(:,:,2) / (2.0_8*dx)
901  matv = t0(:,:,1) - 2.0_8*t0(:,:,3) / (dx*dx)
902  matw = t0(:,:,3)/(dx*dx) - t0(:,:,2) / (2.0_8*dx)
903 
904  do m=1,dimm
905  if (mode(m)==1) cycle
906  if (abs(matw(1,m))>1e-16) then
907  if (verbose>1) print*,"cond1 k=1 m=",m
908  neumann=1
909  if (interpretatif==1.or.dirichlet==1) then
910  if (verbose>1) print*,"boundary conditions not compatible"
911  condlimord2_k1 = 1
912  endif
913  else
914  if (abs(matv(1,m))>1e-16) then
915  if (verbose>1) print*,"cond2 k=1 m=",m
916  dirichlet=1
917  if (neumann==1) then
918  if (verbose>1) print*,"boundary conditions not compatible"
919  condlimord2_k1 = 1
920  endif
921  else
922  if (verbose>1) print*,"boundary conditions not supported"
923  condlimord2_k1 = 1
924  endif
925  endif
926 
927  if (abs(matw(2,m))>1e-16) then
928  if (verbose>1) print*,"cond1p k=1 m=",m
929  neumann=1
930  if (interpretatif==1.or.dirichlet==1) then
931  if (verbose>1) print*,"boundary conditions not compatible"
932  condlimord2_k1 = 1
933  endif
934  else
935  if (abs(matv(2,m))>1e-16) then
936  if (verbose>1) print*,"cond2p k=1 m=",m
937  dirichlet=1
938  if (neumann==1) then
939  if (verbose>1) print*,"boundary conditions not compatible"
940  condlimord2_k1 = 1
941  endif
942  else
943  if (verbose>1) print*,"boundary conditions not supported"
944  condlimord2_k1 =1
945  endif
946  endif
947  enddo
948 
949  if (interpretatif==1) condlimord2_k1=1
950 
951  if (condlimord2_k1 .ne. 0) return
952 
953  do m=1,dimm
954  if (mode(m)==1) cycle
955  if (abs(matw(1,m))>1e-16) then
956  matu(1,m) = - matu(1,m)/matw(1,m)
957  matv(1,m) = - matv(1,m)/matw(1,m)
958  matw(1,m) = v0(1,m) /matw(1,m)
959  do i=1,dimm
960  matb(1,m,i) = matb(1,m,i) + matv(1,m)*mata(1,m,i)
961  matc(1,m,i) = matc(1,m,i) + matu(1,m)*mata(1,m,i)
962  mats(1,m) = mats(1,m) + sca_f*matw(1,m)*mata(1,m,i)
963  enddo
964  !mata(1,m,:) = 0.0
965  else
966  matx(1,m) = -matu(1,m)/matv(1,m)
967  maty(1,m) = v0(1,m) /matv(1,m)
968  matb(1,:,m) = 0.0_8
969  matc(1,:,m) = 0.0_8
970  mats(1,m) = 0.0_8
971  do i=1,dimm
972  matb(2,m,i) = matb(2,m,i) + matx(1,m)*mata(2,m,i)
973  mats(2,i) = mats(2,i) + sca_f*maty(1,m)*mata(2,m,i)
974  enddo
975  mata(2,m,:) = 0.0_8
976  endif
977 
978  if (abs(matw(2,m))>1e-16) then
979  matu(2,m) = - matu(2,m)/matw(2,m)
980  matv(2,m) = - matv(2,m)/matw(2,m)
981  matw(2,m) = v0(2,m) /matw(2,m)
982  do i=1,dimm
983  matbp(1,m,i) = matbp(1,m,i) + matv(2,m)*matap(1,m,i)
984  matcp(1,m,i) = matcp(1,m,i) + matu(2,m)*matap(1,m,i)
985  mats(1,i) = mats(1,i) + (1.0_8-sca_f)*matw(2,m)*matap(1,m,i)
986  enddo
987  !matap(1,m,:) = 0.0
988  else
989  matx(2,m) = -matu(2,m)/matv(2,m)
990  maty(2,m) = v0(2,m) /matv(2,m)
991  matbp(1,:,m) = 0.0_8
992  matcp(1,:,m) = 0.0_8
993  matcp(1,m,m) = matx(2,m)/(1.0_8 - sca_f)
994  mats(1,m) = maty(2,m)
995  do i=1,dimm
996  matbp(2,m,i) = matbp(2,m,i) + matx(2,m)*matap(2,m,i)
997  mats(2,i) = mats(2,i) + (1.0_8-sca_f)*maty(2,m)*matap(2,m,i)
998  enddo
999  matap(2,m,:) = 0.0_8
1000  endif
1001  enddo
1002 
1003 end function condlimord2_k1
1004 !***********************************************************************************
1005 integer function condlimord2_kk(interpretatif)
1006 
1007  !***********************
1008  !boundary conditions k=dimk
1009  !***********************
1010 
1011  use cos_transportdata
1012  implicit none
1013 
1014  integer,intent(in) :: interpretatif
1015  integer :: dirichlet=0,neumann=0,m,i
1016  real(kind=8) :: aux
1017 
1018  condlimord2_kk=0
1019 
1020  matu = t1(:,:,3)/(dx*dx) + t1(:,:,2) / (2.0_8*dx)
1021  matv = t1(:,:,1) - 2.0_8*t1(:,:,3) / (dx*dx)
1022  matw = t1(:,:,3)/(dx*dx) - t1(:,:,2) / (2.0_8*dx)
1023 
1024  do m=1,dimm
1025  if (mode(m)==1) cycle
1026  if (abs(matu(1,m))>1e-16) then
1027  if (verbose>1) print*,"cond1 k=dimk m=",m
1028  neumann=1
1029  if (interpretatif==1.or.dirichlet==1) then
1030  if (verbose>1) print*,"boundary conditions not compatible"
1031  condlimord2_kk=1
1032  endif
1033  else
1034  if (abs(matv(1,m))>1e-16) then
1035  if (verbose>1) print*,"cond2 k=dimk m=",m
1036  dirichlet=1
1037  if (neumann==1) then
1038  if (verbose>1) print*,"boundary conditions not compatible"
1039  condlimord2_kk=1
1040  endif
1041  else
1042  if (verbose>1) print*,"bondary conditions not supported"
1043  condlimord2_kk=1
1044  endif
1045  endif
1046  if (abs(matu(2,m))>1e-16) then
1047  if (verbose>1) print*,"cond1p k=dimk m=",m
1048  neumann=1
1049  if (interpretatif==1.or.dirichlet==1) then
1050  if (verbose>1) print*,"boundary conditions not compatible"
1051  condlimord2_kk=1
1052  endif
1053  else
1054  if (abs(matv(2,m))>1e-16) then
1055  if (verbose>1) print*,"cond2p k=dimk m=",m
1056  if (neumann==1) then
1057  if (verbose>1) print*,"boundary conditions not compatible"
1058  condlimord2_kk=1
1059  endif
1060  else
1061  if (verbose>1) print*,"boundary conditions not supported"
1062  condlimord2_kk=1
1063  endif
1064  endif
1065  enddo
1066  if (interpretatif==1) condlimord2_kk=1
1067 
1068  if (condlimord2_kk.ne.0) return
1069 
1070  do m=1,dimm
1071  if (mode(m)==1) cycle
1072  if (abs(matu(1,m))>1e-16) then
1073  aux=matu(1,m)
1074  matu(1,m) = - matw(1,m)/aux
1075  matv(1,m) = - matv(1,m)/aux
1076  matw(1,m) = v1(1,m) /aux
1077  do i=1,dimm
1078  matb(dimk,m,i) = matb(dimk,m,i) + matv(1,m)*matc(dimk,m,i)
1079  mata(dimk,m,i) = mata(dimk,m,i) + matu(1,m)*matc(dimk,m,i)
1080  mats(dimk,i) = mats(dimk,i) + sca_f*matw(1,m)*matc(dimk,m,i)
1081  enddo
1082  !matc(dimk,m,:) = 0.0
1083  else
1084  matx(1,m) = -matw(1,m)/matv(1,m)
1085  maty(1,m) = v1(1,m) /matv(1,m)
1086  mata(dimk,:,m) = 0.0_8
1087  matb(dimk,:,m) = 0.0_8
1088  mats(dimk,m) = 0.0_8
1089  do i=1,dimm
1090  matb(dimk-1,m,i) = matb(dimk-1,m,i) + matx(1,m)*matc(dimk-1,m,i)
1091  mats(dimk-1,i) = mats(dimk-1,i) + sca_f*maty(1,m)*matc(dimk-1,m,i)
1092  enddo
1093  matc(dimk-1,m,:) = 0.0_8
1094  endif
1095  if (abs(matu(2,m))>1e-16) then
1096  aux=matu(2,m)
1097  matu(2,m) = - matw(2,m)/aux
1098  matv(2,m) = - matv(2,m)/aux
1099  matw(2,m) = v1(2,m) /aux
1100  do i=1,dimm
1101  matbp(dimk,m,i) = matbp(dimk,m,i) + matv(2,m)*matcp(dimk,m,i)
1102  matap(dimk,m,i) = matap(dimk,m,i) + matu(2,m)*matcp(dimk,m,i)
1103  mats(dimk,i) = mats(dimk,i) + (1.0_8-sca_f)*matw(2,m)*matcp(dimk,m,i)
1104  enddo
1105  !matcp(dimk,m,:) = 0.0
1106  else
1107  matx(2,m) = -matw(2,m)/matv(2,m)
1108  maty(2,m) = v1(2,m) /matv(2,m)
1109  matap(dimk,:,m) = 0.0_8
1110  matap(dimk,m,m) = matx(2,m)/(1.0_8 - sca_f)
1111  matbp(dimk,:,m) = 0.0_8
1112  mats(dimk,m) = maty(2,m)
1113  do i=1,dimm
1114  matbp(dimk-1,m,i) = matbp(dimk-1,m,i) + matx(2,m)*matcp(dimk-1,m,i)
1115  mats(dimk-1,i) = mats(dimk-1,i) + (1.0_8-sca_f)*maty(2,m)*matcp(dimk-1,m,i)
1116  enddo
1117  matcp(dimk-1,m,:) = 0.0_8
1118  endif
1119  enddo
1120 
1121 end function condlimord2_kk
1122 !**************************************************************************************
1123 subroutine condlimord1_k1
1124  !***********************
1125  !boundary conditions k=1
1126  !***********************
1127 
1128  use cos_transportdata
1129  implicit none
1130 
1131  integer :: m,i
1132  real(kind=8) :: aux
1133 
1134  matu = t0(:,:,2)
1135  matv = t0(:,:,1)*dx - t0(:,:,2)
1136  maty = v0(:,:)*dx
1137  print*,"dx",dx,t0(2,1,1), t0(2,1,2)
1138  print*,"matv",matv(2,1)
1139 
1140  do m=1,dimm
1141  if (mode(m)==1) cycle
1142  matx(1,m) = -matu(1,m)/matv(1,m)
1143  maty(1,m) = maty(1,m)/matv(1,m)
1144  matb(1,:,m) = 0.0_8
1145  matc(1,:,m) = 0.0_8
1146  mats(1,m) = 0.0_8
1147  do i=1,dimm
1148  matb(2,m,i) = matb(2,m,i) + matx(1,m)*mata(2,m,i)
1149  mats(2,i) = mats(2,i) + sca_f*maty(1,m)*mata(2,m,i)
1150  enddo
1151  mata(2,m,:) = 0.0_8
1152  if (abs(matv(2,m))>1e-16) then
1153  matx(2,m) = -matu(2,m)/matv(2,m)
1154  maty(2,m) = maty(2,m)/matv(2,m)
1155  matbp(1,:,m) = 0.0_8
1156  matcp(1,:,m) = 0.0_8
1157  matcp(1,m,m) = matx(2,m)/(1.0_8 - sca_f)
1158  mats(1,m) = maty(2,m)
1159  do i=1,dimm
1160  matbp(2,m,i) = matbp(2,m,i) + matx(2,m)*matap(2,m,i)
1161  mats(2,i) = mats(2,i) + (1.0_8-sca_f)*maty(2,m)*matap(2,m,i)
1162  enddo
1163  matap(2,m,:) = 0.0_8
1164  else
1165  print*,"Boundaries conditions error",matv(2,m)
1166  endif
1167  enddo
1168 end subroutine condlimord1_k1
1169 !**********************************************************************************************
1170 subroutine condlimord1_kk
1171  !**************************
1172  !boundary conditions k=dimk
1173  !**************************
1174 
1175  use cos_transportdata
1176  implicit none
1177 
1178  integer :: m,i
1179  real(kind=8) :: aux
1180 
1181  matv = t1(:,:,2) + dx*t1(:,:,1)
1182  matw = - t1(:,:,2)
1183  maty = v1(:,:)*dx
1184 
1185  do m=1,dimm
1186  if (mode(m)==1) cycle
1187  matx(1,m) = -matw(1,m)/matv(1,m)
1188  maty(1,m) = maty(1,m)/matv(1,m)
1189  mata(dimk,:,m) = 0.0_8
1190  matb(dimk,:,m) = 0.0_8
1191  mats(dimk,m) = 0.0_8
1192  do i=1,dimm
1193  matb(dimk-1,m,i) = matb(dimk-1,m,i) + matx(1,m)*matc(dimk-1,m,i)
1194  mats(dimk-1,i) = mats(dimk-1,i) + sca_f*maty(1,m)*matc(dimk-1,m,i)
1195  enddo
1196  matc(dimk-1,m,:) = 0.0_8
1197  if (abs(matv(2,m))>1e-16) then
1198  matx(2,m) = -matw(2,m)/matv(2,m)
1199  maty(2,m) = maty(2,m)/matv(2,m)
1200  matap(dimk,:,m) = 0.0_8
1201  matap(dimk,m,m) = matx(2,m)/(1.0_8 - sca_f)
1202  matbp(dimk,:,m) = 0.0_8
1203  mats(dimk,m) = maty(2,m)
1204  do i=1,dimm
1205  matbp(dimk-1,m,i) = matbp(dimk-1,m,i) + matx(2,m)*matcp(dimk-1,m,i)
1206  mats(dimk-1,i) = mats(dimk-1,i) + (1.0_8-sca_f)*maty(2,m)*matcp(dimk-1,m,i)
1207  enddo
1208  matcp(dimk-1,m,:) = 0.0_8
1209  else
1210  print*,"Boundaries conditions error",matv(2,m)
1211  endif
1212  enddo
1213 
1214 end subroutine condlimord1_kk
1215 
1216 #endif
cos_precision
types
Definition: cos_precision.f90:9
solver
Definition: solver.f90:1
type_solver
Definition: type_solver.f90:1
solution10
subroutine solution10(SOLVER, ifail)
This subroutine is prepared to solve single transport equation in standardised form adopted by the ET...
Definition: solution10.F90:9
cos_zconversion
subroutine cos_zconversion(e, nbrho)
Definition: cos_zconversion.f90:1
type_solver::numerics
Definition: type_solver.f90:22