ETS-Core  version:0.0.4-46-ge2d8
Core actors for the ETS-6
 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 ids_schemas
27  USE type_transport_solver_numerics
28 
29 
30  IMPLICIT NONE
31 
32  INTEGER, INTENT (INOUT) :: ifail
33 
34 ! +++ Input/Output with ETS:
35  TYPE (solver_io), 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 (IDS_REAL) :: rho(solver%nrho) !radii
47  REAL (IDS_REAL) :: rhomax, rhomax2 !AF 9.Jul.2010
48 
49  REAL (IDS_REAL) :: amix !fraction of new sollution mixed
50 
51  REAL (IDS_REAL) :: y(solver%nrho), ym(solver%nrho) !function at the current amd previous time steps
52  REAL (IDS_REAL) :: dy(solver%nrho) !derivative of function
53  REAL (IDS_REAL) :: yy(solver%nrho) !new function !AF 12.Jul.2010
54  REAL (IDS_REAL) :: dyy(solver%nrho) !derivative of new function !AF 12.Jul.2010
55 
56  REAL (IDS_REAL) :: a(solver%nrho), b(solver%nrho) !coefficients
57  REAL (IDS_REAL) :: c(solver%nrho), d(solver%nrho) !coefficients
58  REAL (IDS_REAL) :: e(solver%nrho), f(solver%nrho) !coefficients
59  REAL (IDS_REAL) :: g(solver%nrho), h !coefficients
60  REAL (IDS_REAL) :: dd(solver%nrho), de(solver%nrho) !derivatives of coefficients
61 
62  REAL (IDS_REAL) :: 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 (IDS_REAL) :: x(solver%nrho) !radii
70 
71  REAL (IDS_REAL) :: sol(solver%nrho) !solution
72  REAL (IDS_REAL) :: dsol(solver%nrho) !derivative of solution
73 
74  REAL (IDS_REAL) :: as(solver%nrho), bs(solver%nrho) !coefficients
75  REAL (IDS_REAL) :: cs(solver%nrho) !coefficients
76 
77  REAL (IDS_REAL) :: acron(solver%nrho), bcron(solver%nrho) !coefficients
78  REAL (IDS_REAL) :: ccron(solver%nrho), dcron(solver%nrho) !coefficients
79 
80 ! REAL (IDS_REAL) :: mata_in(SOLVER%NRHO,5,5), matb_in(SOLVER%NRHO,5,5) !coefficients ! AF, 10.Jul.2010
81 ! REAL (IDS_REAL) :: matc_in(SOLVER%NRHO,5,5), matd_in(SOLVER%NRHO,5) !coefficients ! AF, 10.Jul.2010
82  REAL (IDS_REAL) :: mata_in(solver%nrho,1,1), matb_in(solver%nrho,1,1) !coefficients ! AF, 10.Jul.2010
83  REAL (IDS_REAL) :: matc_in(solver%nrho,1,1), matd_in(solver%nrho,1) !coefficients ! AF, 10.Jul.2010
84 
85  ! REAL (IDS_REAL) :: VS(2), US(2), WS(2) !boundary conditions ! AF, 10.Jul.2010
86 
87  integer, parameter :: nbeq = 1 !AF 9.Jul.2010
88  REAL(IDS_REAL) :: dt_in !AF 9.Jul.2010
89  REAL(IDS_REAL) :: dx_in !AF 9.Jul.2010
90  REAL(IDS_REAL) :: sca_f_in !AF 9.Jul.2010
91  REAL(IDS_REAL), dimension(2,nbeq,3) :: t0_in,t1_in !AF 9.Jul.2010
92  REAL(IDS_REAL), dimension(2,nbeq) :: v0_in,v1_in !AF 9.Jul.2010
93 
94  REAL(IDS_REAL), 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_ids_real
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_ids_real
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_ids_real/(nrho-1) !0.01 !AF - 10.Jul.2010
244 ! 0 implicite, 1 explicite 0.5 Krank Nicholson
245  sca_f_in = 0.e0_ids_real !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_ids_real-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_IDS_REAL-AMIX) + SOL(IRHO)*AMIX !AF - 9.Jul.2010
266  ! DY(IRHO) = DY(IRHO)*(1.e0_IDS_REAL-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 ids_schemas
309 
310 
311  IMPLICIT NONE
312 
313  INTEGER :: n ! number of radial points (input)
314  INTEGER :: i
315 
316  REAL (IDS_REAL) :: x(n), & ! argument array (input)
317  y(n), & ! function array (input)
318  dy1(n) ! function derivative array (output)
319  REAL (IDS_REAL) :: h(n),dy2(n)
320 
321  DO i=1,n-1
322  h(i)=x(i+1)-x(i)
323  END DO
324 
325  DO i=2,n-1
326  dy1(i)=((y(i+1)-y(i))*h(i-1)/h(i)+(y(i)-y(i-1))*h(i)/h(i-1)) &
327  /(h(i)+h(i-1))
328  dy2(i)=2.e0_ids_real*((y(i-1)-y(i))/h(i-1)+(y(i+1)-y(i))/h(i)) &
329  /(h(i)+h(i-1))
330  END DO
331 
332  dy1(1)=dy1(2)-dy2(2)*h(1)
333  dy1(n)=dy1(n-1)+dy2(n-1)*h(n-1)
334 
335  RETURN
336 END SUBROUTINE derivn10
337 
338 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
340 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
341 
342 subroutine cos_zconversion(e,nbrho)
343  use cos_precision
344  implicit none
345 
346  integer,intent(in) :: nbrho
347  real(DP),dimension(nbrho) :: e
348 
349  if (nbrho < 4) then
350  print*, "error"
351  return
352  endif
353 
354  ! continuite au centre pour les NaN
355  if ( (abs(e(1)).gt.huge(e(1))) .or. (abs(e(1)).eq.0) ) then
356  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)
357  endif
358 
359 end subroutine cos_zconversion
360 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
361 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
362 
363 
364 !********************************************************************************
365 module cos_transportdata
366  implicit none
367  integer :: dimk,dimm,dimmat,dimbb,dimn
368  real(kind=8),dimension(:,:), pointer :: f,fp,dfpdx,dfpdx2
369  real(kind=8),dimension(:,:),allocatable :: matd,matdp
370  real(kind=8),dimension(:,:,:),pointer :: T0,T1
371  real(kind=8),dimension(:,:,:),allocatable :: mata,matb,matc, &
372  matap,matbp,matcp
373  real(kind=8),dimension(:,:), pointer :: V0,V1
374  integer,dimension(:),allocatable :: mode
375 
376  real(kind=8),dimension(:,:),allocatable :: mats
377  real(kind=8),dimension(:,:),allocatable :: matu,matv,matw,matx,maty,matz
378  real(kind=8),dimension(:),allocatable,target:: mat_coo,mat_csr
379  real(kind=8),dimension(:),allocatable,target:: vectbb
380  integer,dimension(:),allocatable,target :: cooi,cooj,csri,csrj
381  real(kind=8), pointer:: dx,dt,sca_f
382  integer :: verbose = 0
383 end module cos_transportdata
384 !module cos_pde1dsolver_interface ! AF - 9.Jul.2010, this interface is not in use
385 ! interface
386 ! subroutine cos_pde1solver(f_in,fp_in, &
387 ! mata_in,matb_in,matc_in,matd_in, &
388 ! matap_in,matbp_in,matcp_in,matdp_in, &
389 ! dimk_in,dimn_in,dimm_in,dx_in,dt_in,sca_f_in, &
390 ! T0_in,T1_in,V0_in,V1_in,mode_in,dfpdx_in,dfpdx2_in)
391 ! use cos_transportdata
392 ! implicit none
393 ! integer, intent(in) :: dimk_in,dimn_in,dimm_in
394 ! real(kind=8),intent(in),target :: dt_in,dx_in,sca_f_in
395 ! real(kind=8),dimension(dimk_in,dimm_in),target :: f_in,fp_in
396 ! real(kind=8),dimension(dimk_in,dimm_in),target,optional:: dfpdx_in,dfpdx2_in
397 ! real(kind=8),dimension(dimk_in,dimm_in),target :: matd_in,matdp_in
398 ! real(kind=8),dimension(dimk_in,dimm_in,dimm_in),target :: mata_in,matb_in,matc_in, &
399 ! matap_in,matbp_in,matcp_in
400 ! real(kind=8),dimension(2,dimm_in,3),target :: T0_in,T1_in
401 ! real(kind=8),dimension(2,dimm_in),target :: V0_in,V1_in
402 ! real(kind=8),dimension(dimm_in) :: mode_in
403 ! end subroutine cos_pde1solver
404 ! end interface
405 !end module cos_pde1dsolver_interface
406 !*********************************************************************************
407 subroutine cos_pde1solver(f_in,fp_in, &
408  mata_in,matb_in,matc_in,matd_in, &
409  matap_in,matbp_in,matcp_in,matdp_in, &
410  dimk_in,dimn_in,dimm_in,dx_in,dt_in,sca_f_in, &
411  t0_in,t1_in,v0_in,v1_in,mode_in,dfpdx_in,dfpdx2_in)
412  !*********************************************************************
413  !solve the equations :
414  ! Diff(f,t) = mata * diff(f,x,x) + matb * diff(f,x) + matc * f + matd
415  !
416  !dt_in : time step
417  !dx_in : space step
418  !
419  !f : solution at the time t
420  !fp : solution at the time t+dt
421  !dfpdx : diff(fp,x)
422  !dfpdx2 : diff(fp,x,x)
423  !
424  !dimk_in : dimension of space index
425  ! x=x0 + (k-1)*dx, k->[1..dimk]
426  !dimm_in : dimension of equation number
427  !dimn_in : dimension of the solution at a one grid point
428  !
429  !sca_f_in : scalar
430  ! 0 -> full implicit
431  ! 0.5 -> Crank Nicholson
432  ! 1 -> explicit (don't used it)
433  !
434  !
435  !The big matrix
436  !--------------
437  !
438  ! k=1 k=2 ----------------------------------------- k=K
439  ! (1,1)--------(1,m)(1,1)--------(1,m)
440  ! | b(1,m,n) | | c(1,m,n) |
441  ! | | | |
442  ! (n,1)--------(n,m)(n,1)--------(n,m)
443  ! (1,1)--------(1,m)(1,1)--------(1,m)(1,1)--------(1,m)
444  ! | a(2,m,n) | | b(2,m,n) | | c(2,m,n) |
445  ! | | | | | |
446  ! (n,m)--------(n,m)(n,m)--------(n,m)(n,m)--------(n,m)
447  ! .
448  ! .
449  ! .
450  ! .
451  ! .
452  ! .
453  ! .
454  ! .
455  ! .
456  !
457  ! (1,1)--------(1,m)(1,1)--------(1,m)
458  ! | a(K,m,n) | | b(K,m,n) |
459  ! | | | |
460  ! (n,m)--------(n,m)(n,m)--------(n,m)
461  !
462  ! where n=m if we are in predictive mode
463  ! n<m if we are in interpretative mode
464  ! the first index is the space index
465  ! the second index is the coupling coefficient between variables
466  ! the third index is the rank of the equation
467  !
468  ! T0(l,m,j) : coefficients of the boundaries condition at x0
469  ! V0(l,m) : values of the boundaries condition at x0
470  ! T0(:,:,3) * diff(f,x,x) + T0(:,:,2) * diff(f,x) + T0(:,:,1) * f = V0
471  ! where l=1 for t
472  ! l=2 for t+dt
473  ! T1(l,m,j) , V1(l,m) : idem as above but at x1
474  !*************************************************************************************
475  use cos_transportdata
476  implicit none
477  integer, intent(in) :: dimk_in,dimn_in,dimm_in
478  real(kind=8),intent(in),target :: dt_in,dx_in,sca_f_in
479  real(kind=8),dimension(dimk_in,dimm_in),target :: f_in,fp_in
480  real(kind=8),dimension(dimk_in,dimm_in),target,optional:: dfpdx_in,dfpdx2_in
481  real(kind=8),dimension(dimk_in,dimm_in),target :: matd_in,matdp_in
482  real(kind=8),dimension(dimk_in,dimm_in,dimm_in),target :: mata_in,matb_in,matc_in, &
483  matap_in,matbp_in,matcp_in
484  real(kind=8),dimension(2,dimm_in,3),target :: t0_in,t1_in
485  real(kind=8),dimension(2,dimm_in),target :: v0_in,v1_in
486 
487  real(kind=8),dimension(:,:,:), allocatable :: mat_aux
488 
489  real(kind=8),dimension(dimm_in) :: mode_in
490  integer,dimension(:),allocatable :: iwk,cooi_aux
491 
492  real(kind=8) :: dtx2,dt2x,aux
493  integer :: i,j,k,m,n
494  integer :: counteri,counterj,counter
495  integer :: interpretatif=0,condlim=0
496 
497  interface
498  integer function condlimord2_k1(interpretatif)
499  integer,intent(in) :: interpretatif
500  end function condlimord2_k1
501 
502  integer function condlimord2_kk(interpretatif)
503  integer,intent(in) :: interpretatif
504  end function condlimord2_kk
505  end interface
506 
507 
508  write(*,*) 'dans solver',sca_f_in
509  if (sca_f_in==1.0) then
510  print*,"Error, the explicit Euler algorithm must not be used"
511  stop
512  endif
513 
514  !**********
515  !allocation
516  !**********
517  dimk=dimk_in
518  dimm=dimm_in
519  dimn=dimn_in
520  sca_f=>sca_f_in
521  dx=>dx_in
522  dt=>dt_in
523  f=>f_in
524  fp=>fp_in
525  if (present(dfpdx_in)) dfpdx=>dfpdx_in
526  if (present(dfpdx2_in)) dfpdx2=>dfpdx2_in
527  allocate(mata(dimk_in,dimm_in,dimm_in),matb(dimk_in,dimm_in,dimm_in),&
528  matc(dimk_in,dimm_in,dimm_in),matd(dimk_in,dimm_in) )
529  allocate(matap(dimk_in,dimm_in,dimm_in),matbp(dimk_in,dimm_in,dimm_in),&
530  matcp(dimk_in,dimm_in,dimm_in),matdp(dimk_in,dimm_in) )
531  allocate(mode(dimm_in))
532 
533  mata=mata_in
534  matb=matb_in
535  matc=matc_in
536  matd=matd_in
537 
538  matap=matap_in
539  matbp=matbp_in
540  matcp=matcp_in
541  matdp=matdp_in
542  !print*,"matap",matap(1:10,1,1)
543 
544  mode=int(mode_in)
545  t0=>t0_in
546  t1=>t1_in
547  v0=>v0_in
548  v1=>v1_in
549  allocate(mats(dimk,dimm))
550  allocate(matu(2,dimm),matv(2,dimm),matw(2,dimm))
551  allocate(matx(2,dimm),maty(2,dimm),matz(2,dimm))
552  allocate(mat_aux(dimk,dimm,dimm))
553 
554  dimmat=(dimk-2)*dimn*dimn*3 +4*dimn*dimn
555  dimbb=dimn*dimk
556  allocate(vectbb(dimbb),mat_coo(dimmat),cooi(dimmat),cooj(dimmat))
557 
558  !**************
559  !initialisation
560  !**************
561  interpretatif=sum(mode)
562  if (interpretatif.ne.0) interpretatif=1
563  dtx2 = dt/(dx*dx)
564  dt2x = dt/(2.0_8*dx)
565  write(*,*) 'dtx2=',dtx2
566  mat_aux = matb
567  matb = -2.0_8*dtx2*mata + dt*matc
568  matc = dtx2*mata + dt2x*mat_aux
569  mata = dtx2*mata - dt2x*mat_aux
570  matd = dt*matd
571 
572  mat_aux = matbp
573  matbp = -2.0_8*dtx2*matap + dt*matcp
574  matcp = dtx2*matap + dt2x*mat_aux
575  matap = dtx2*matap - dt2x*mat_aux
576  matdp = dt*matdp
577 
578  deallocate(mat_aux)
579 
580  mats = sca_f*matd + (1.0_8-sca_f)*matdp + f;
581 
582  !*********************
583  !boundaries conditions
584  !*********************
585  condlim=sum(t0(:,:,3))
586  if ( condlimord2_k1(interpretatif) > 0) then
587  if (condlim>0) then
588  print*," boundaries condition of order 2 is not implemented"
589  stop
590  else
591  !print*,"boundaries condition of order 1 k=1"
592  call condlimord1_k1
593  endif
594  endif
595 
596  condlim=sum(t1(:,:,3))
597  if ( condlimord2_kk(interpretatif) > 0) then
598  if (condlim>0) then
599  print*," boundaries condition of order 2 is not implemented"
600  stop
601  else
602  !print*,"boundaries condition of order 1 k=dimk"
603  call condlimord1_kk
604  endif
605  endif
606 
607 
608  !*******************
609  !interpretative mode
610  !*******************
611  write(*,*) 'dimm=',dimm
612  do m = 1,dimm
613  if (mode(m)==1) cycle
614  do k=2,dimk-1
615  do n=1,dimm
616  if (mode(n)==1) then
617  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)) &
618  + (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))
619  endif
620  enddo
621  enddo
622  enddo
623 
624  !*************************
625  !filling the second member
626  !*************************
627  !for k=1
628  k=1
629  counter=0
630  do m=1,dimm
631  if (mode(m)==0) then
632  counter=counter+1
633  vectbb(counter) = mats(k,m)
634  do n=1,dimm
635  if (mode(n)==0) then
636  vectbb(counter) = vectbb(counter) + sca_f*(matb(k,n,m)*f(k,n) + matc(k,n,m)*f(k+1,n))
637  endif
638  enddo
639  endif
640  enddo
641 
642  !for k =[2,dimk-1]
643  do k=2,dimk-1
644  do m=1,dimm
645  if (mode(m)==0) then
646  counter=counter+1
647  vectbb(counter) = mats(k,m)
648  do n=1,dimm
649  if (mode(n)==0) then
650  vectbb(counter) = vectbb(counter) + sca_f*(mata(k,n,m)*f(k-1,n) + matb(k,n,m)*f(k,n) &
651  + matc(k,n,m)*f(k+1,n))
652  endif
653  enddo
654  endif
655  enddo
656  enddo
657 
658  !for k=dimk
659  k=dimk
660  do m=1,dimm
661  if (mode(m)==0) then
662  counter=counter+1
663  vectbb(counter) = mats(k,m)
664  do n=1,dimm
665  if (mode(n)==0) then
666  vectbb(counter) = vectbb(counter) + sca_f*(mata(k,n,m)*f(k-1,n) + matb(k,n,m)*f(k,n))
667  endif
668  enddo
669  endif
670  enddo
671 
672  !***********************
673  !filling the big matrix
674  !***********************
675  k=1
676  counter =0
677  counteri=0
678  counterj=0
679  write(*,*) 'counter=',counter
680  call build_matcoo(matbp,1,k,counter,counteri,counterj)
681  counteri=0
682  counterj=dimn
683  call build_matcoo(matcp,0,k,counter,counteri,counterj)
684  do k=2,dimk-1
685  counteri=(k-1)*dimn
686  counterj=(k-2)*dimn
687  call build_matcoo(matap,0,k,counter,counteri,counterj)
688  counteri=(k-1)*dimn
689  counterj=(k-1)*dimn
690  call build_matcoo(matbp,1,k,counter,counteri,counterj)
691  counteri=(k-1)*dimn
692  counterj=k*dimn
693  call build_matcoo(matcp,0,k,counter,counteri,counterj)
694  enddo
695  k=dimk
696  counteri=(k-1)*dimn
697  counterj=(k-2)*dimn
698  call build_matcoo(matap,0,k,counter,counteri,counterj)
699  counteri=(k-1)*dimn
700  counterj=(k-1)*dimn
701  call build_matcoo(matbp,1,k,counter,counteri,counterj)
702  !*****************
703  ! print the matrix
704  !*****************
705  if (verbose>0) then
706  print*,"dimk",dimk,dimn,dimm,counter
707  allocate(iwk(dimbb+1),cooi_aux(dimmat))
708  allocate(mat_csr(dimmat),csri(dimmat),csrj(dimmat))
709  endif
710  dimmat=counter
711  if (verbose>0) then
712  print*,"SPARSE BIG MATRIX:"
713  cooi_aux=cooi
714  ! call coicsr (dimbb,dimmat,1,mat_coo,cooj,cooi,iwk)
715  !call coocsr(dimbb,dimmat,mat_coo,cooi_aux,cooj,mat_csr,csrj,csri)
716  ! call dump (1,dimbb,1,mat_coo,cooj,cooi,6)
717  !call dump (1,dimbb,1,mat_csr,csrj,csri,6)
718  print*,"RHS:"
719  print*,vectbb
720  endif
721  !****************
722  ! solve de matrix
723  !****************
724  call solvepde1d
725 
726  !*****************
727  ! save de solution
728  !*****************
729  counter=0
730  do k=1,dimk
731  do m=1,dimm
732  if (mode(m)==0) then
733  counter=counter+1
734  fp(k,m) = vectbb(counter)
735  endif
736  enddo
737  enddo
738 
739  !fp(1,1)=mat_coo(1)
740  !fp(2,1)=mat_coo(4)
741  !fp(3,1)=aux
742  !fp(4,1)=dt2x
743 
744  !do m=1,dimm
745  ! if (mode(m) == 1) then
746  ! do k=1,dimk
747  ! fp(k,m) = fp_in(k,m)
748  ! enddo
749  ! endif
750  !enddo
751  if (verbose > 0) then
752  print*,"SOLUTION:"
753  print*,vectbb
754  endif
755 
756  !************
757  !finalisation
758  !************
759  deallocate(matu,matv,matw,matx,maty,matz,mats)
760  deallocate(mata,matb,matc,matd,matap,matbp,matcp,matdp)
761  deallocate(mode)
762  if (verbose==0) then
763  deallocate(mat_coo,cooi,cooj,vectbb)
764  else
765  deallocate(mat_coo,cooi,cooj,vectbb,iwk,cooi_aux,mat_csr,csri,csrj)
766  endif
767 
768  !if (dfpdx(1,1)==0.and.dfpdx2(1,1)==0) then
769  !print*,"dfpdx and dfpdx2 not available"
770  !endif
771 
772 end subroutine cos_pde1solver
773 !**********************************************************************************
774 subroutine build_matcoo(matp,option,k,counter,counteri,counterj)
775  use cos_transportdata
776  implicit none
777  integer,intent(in) :: option,k
778  integer,intent(inout) :: counter,counteri,counterj
779  integer :: counteri_aux,counterj_aux
780  real(kind=8),dimension(dimk,dimm,dimm) :: matp
781 
782  integer :: m,n
783 
784  counteri_aux=counteri
785  counterj_aux=counterj
786  do m=1,dimm
787  !print*,"mod",mode(1)
788  if (mode(m)==0) then
789  counteri=counteri+1
790  counterj=counterj_aux
791  do n=1,dimm
792  if (mode(n)==0) then
793  counterj=counterj+1
794  counter=counter+1
795  if (counteri==counterj) then
796  mat_coo(counter)=1.0_8-(1.0_8-sca_f)*matp(k,n,m)
797  else
798  mat_coo(counter)=-(1.0_8-sca_f)*matp(k,n,m)
799  endif
800  if (mat_coo(counter)==0.0) then
801  counter=counter-1
802  !print*,"vide",counter,sca_f,matp(k,n,m)
803  else
804  cooi(counter)=counteri
805  cooj(counter)=counterj
806  endif
807  endif
808  enddo
809  endif
810  enddo
811  !print*,"counter",counter
812 end subroutine build_matcoo
813 !************************************************************************************
814 subroutine solvepde1d
815  use cos_transportdata
816  implicit none
817  include 'mpif.h'
818  include 'dmumps_struc.h'
819  type(dmumps_struc) :: mumps_par
820  integer:: ierr, i
821 
822  call mpi_init(ierr)
823  ! Define a communicator for the package.
824  mumps_par%COMM = mpi_comm_world
825  ! Initialize an instance of the package
826  ! for L U factorization (sym = 0, with working host)
827  mumps_par%JOB = -1
828  mumps_par%SYM = 0
829  mumps_par%PAR = 1
830  call dmumps(mumps_par)
831 
832  ! Define problem on the host (processor 0)
833  IF ( mumps_par%MYID .eq. 0 ) THEN
834  !order of the matrix
835  mumps_par%N = dimk*dimn
836  !number of entries
837  mumps_par%NZ = dimmat
838  !allocate( mumps_par%IRN ( mumps_par%NZ ) )
839  !allocate( mumps_par%JCN ( mumps_par%NZ ) )
840  !allocate( mumps_par%A( mumps_par%NZ ) )
841  !allocate( mumps_par%RHS ( mumps_par%N ) )
842  !DO I = 1, mumps_par%NZ
843  ! mumps_par%IRN(I) = cooi(i)
844  ! mumps_par%JCN(I) = cooj(i)
845  ! mumps_par%A(I) = mat_coo(i)
846  !END DO
847  !DO I = 1, mumps_par%N
848  ! mumps_par%RHS(I) = vectbb(i)
849  !END DO
850  mumps_par%IRN=>cooi
851  mumps_par%JCN=>cooj
852  mumps_par%A=>mat_coo
853  mumps_par%RHS=>vectbb
854  END IF
855  if (verbose <= 1) then
856  mumps_par%ICNTL(4)=1
857  mumps_par%ICNTL(3)=0
858  !only because we have a problem with mexfile
859  !mumps_par%ICNTL(1)=0
860  endif
861 
862  ! do i = 1,40
863  ! print*,"icntl",i,mumps_par%ICNTL(i)
864  ! enddo
865  ! Call package for solution
866  mumps_par%JOB = 6
867  CALL dmumps(mumps_par)
868  ! Solution has been assembled on the host
869  IF ( mumps_par%MYID .eq. 0 .and. verbose > 1) THEN
870  WRITE( 6, * ) ' MUMPS Solution is ',(mumps_par%RHS(i),i=1,mumps_par%N)
871  END IF
872  !vectbb=mumps_par%RHS
873  ! Deallocate user data
874  IF ( mumps_par%MYID .eq. 0 )THEN
875  !DEALLOCATE( mumps_par%IRN )
876  !DEALLOCATE( mumps_par%JCN )
877  !DEALLOCATE( mumps_par%A )
878  !DEALLOCATE( mumps_par%RHS )
879  END IF
880  ! Destroy the instance (deallocate internal data structures)
881  mumps_par%JOB = -2
882  CALL dmumps(mumps_par)
883  CALL mpi_finalize(ierr)
884 end subroutine solvepde1d
885 !**********************************************************************************
886 integer function condlimord2_k1(interpretatif)
887 
888  !***********************
889  !boundary conditions k=1
890  !***********************
891 
892  use cos_transportdata
893  implicit none
894 
895  integer,intent(in) :: interpretatif
896  integer :: dirichlet=0,neumann=0,m,i
897  real(kind=8) :: aux
898 
899  condlimord2_k1=0
900 
901  matu = t0(:,:,3)/(dx*dx) + t0(:,:,2) / (2.0_8*dx)
902  matv = t0(:,:,1) - 2.0_8*t0(:,:,3) / (dx*dx)
903  matw = t0(:,:,3)/(dx*dx) - t0(:,:,2) / (2.0_8*dx)
904 
905  do m=1,dimm
906  if (mode(m)==1) cycle
907  if (abs(matw(1,m))>1e-16) then
908  if (verbose>1) print*,"cond1 k=1 m=",m
909  neumann=1
910  if (interpretatif==1.or.dirichlet==1) then
911  if (verbose>1) print*,"boundary conditions not compatible"
912  condlimord2_k1 = 1
913  endif
914  else
915  if (abs(matv(1,m))>1e-16) then
916  if (verbose>1) print*,"cond2 k=1 m=",m
917  dirichlet=1
918  if (neumann==1) then
919  if (verbose>1) print*,"boundary conditions not compatible"
920  condlimord2_k1 = 1
921  endif
922  else
923  if (verbose>1) print*,"boundary conditions not supported"
924  condlimord2_k1 = 1
925  endif
926  endif
927 
928  if (abs(matw(2,m))>1e-16) then
929  if (verbose>1) print*,"cond1p k=1 m=",m
930  neumann=1
931  if (interpretatif==1.or.dirichlet==1) then
932  if (verbose>1) print*,"boundary conditions not compatible"
933  condlimord2_k1 = 1
934  endif
935  else
936  if (abs(matv(2,m))>1e-16) then
937  if (verbose>1) print*,"cond2p k=1 m=",m
938  dirichlet=1
939  if (neumann==1) then
940  if (verbose>1) print*,"boundary conditions not compatible"
941  condlimord2_k1 = 1
942  endif
943  else
944  if (verbose>1) print*,"boundary conditions not supported"
945  condlimord2_k1 =1
946  endif
947  endif
948  enddo
949 
950  if (interpretatif==1) condlimord2_k1=1
951 
952  if (condlimord2_k1 .ne. 0) return
953 
954  do m=1,dimm
955  if (mode(m)==1) cycle
956  if (abs(matw(1,m))>1e-16) then
957  matu(1,m) = - matu(1,m)/matw(1,m)
958  matv(1,m) = - matv(1,m)/matw(1,m)
959  matw(1,m) = v0(1,m) /matw(1,m)
960  do i=1,dimm
961  matb(1,m,i) = matb(1,m,i) + matv(1,m)*mata(1,m,i)
962  matc(1,m,i) = matc(1,m,i) + matu(1,m)*mata(1,m,i)
963  mats(1,m) = mats(1,m) + sca_f*matw(1,m)*mata(1,m,i)
964  enddo
965  !mata(1,m,:) = 0.0
966  else
967  matx(1,m) = -matu(1,m)/matv(1,m)
968  maty(1,m) = v0(1,m) /matv(1,m)
969  matb(1,:,m) = 0.0_8
970  matc(1,:,m) = 0.0_8
971  mats(1,m) = 0.0_8
972  do i=1,dimm
973  matb(2,m,i) = matb(2,m,i) + matx(1,m)*mata(2,m,i)
974  mats(2,i) = mats(2,i) + sca_f*maty(1,m)*mata(2,m,i)
975  enddo
976  mata(2,m,:) = 0.0_8
977  endif
978 
979  if (abs(matw(2,m))>1e-16) then
980  matu(2,m) = - matu(2,m)/matw(2,m)
981  matv(2,m) = - matv(2,m)/matw(2,m)
982  matw(2,m) = v0(2,m) /matw(2,m)
983  do i=1,dimm
984  matbp(1,m,i) = matbp(1,m,i) + matv(2,m)*matap(1,m,i)
985  matcp(1,m,i) = matcp(1,m,i) + matu(2,m)*matap(1,m,i)
986  mats(1,i) = mats(1,i) + (1.0_8-sca_f)*matw(2,m)*matap(1,m,i)
987  enddo
988  !matap(1,m,:) = 0.0
989  else
990  matx(2,m) = -matu(2,m)/matv(2,m)
991  maty(2,m) = v0(2,m) /matv(2,m)
992  matbp(1,:,m) = 0.0_8
993  matcp(1,:,m) = 0.0_8
994  matcp(1,m,m) = matx(2,m)/(1.0_8 - sca_f)
995  mats(1,m) = maty(2,m)
996  do i=1,dimm
997  matbp(2,m,i) = matbp(2,m,i) + matx(2,m)*matap(2,m,i)
998  mats(2,i) = mats(2,i) + (1.0_8-sca_f)*maty(2,m)*matap(2,m,i)
999  enddo
1000  matap(2,m,:) = 0.0_8
1001  endif
1002  enddo
1003 
1004 end function condlimord2_k1
1005 !***********************************************************************************
1006 integer function condlimord2_kk(interpretatif)
1007 
1008  !***********************
1009  !boundary conditions k=dimk
1010  !***********************
1011 
1012  use cos_transportdata
1013  implicit none
1014 
1015  integer,intent(in) :: interpretatif
1016  integer :: dirichlet=0,neumann=0,m,i
1017  real(kind=8) :: aux
1018 
1019  condlimord2_kk=0
1020 
1021  matu = t1(:,:,3)/(dx*dx) + t1(:,:,2) / (2.0_8*dx)
1022  matv = t1(:,:,1) - 2.0_8*t1(:,:,3) / (dx*dx)
1023  matw = t1(:,:,3)/(dx*dx) - t1(:,:,2) / (2.0_8*dx)
1024 
1025  do m=1,dimm
1026  if (mode(m)==1) cycle
1027  if (abs(matu(1,m))>1e-16) then
1028  if (verbose>1) print*,"cond1 k=dimk m=",m
1029  neumann=1
1030  if (interpretatif==1.or.dirichlet==1) then
1031  if (verbose>1) print*,"boundary conditions not compatible"
1032  condlimord2_kk=1
1033  endif
1034  else
1035  if (abs(matv(1,m))>1e-16) then
1036  if (verbose>1) print*,"cond2 k=dimk m=",m
1037  dirichlet=1
1038  if (neumann==1) then
1039  if (verbose>1) print*,"boundary conditions not compatible"
1040  condlimord2_kk=1
1041  endif
1042  else
1043  if (verbose>1) print*,"bondary conditions not supported"
1044  condlimord2_kk=1
1045  endif
1046  endif
1047  if (abs(matu(2,m))>1e-16) then
1048  if (verbose>1) print*,"cond1p k=dimk m=",m
1049  neumann=1
1050  if (interpretatif==1.or.dirichlet==1) then
1051  if (verbose>1) print*,"boundary conditions not compatible"
1052  condlimord2_kk=1
1053  endif
1054  else
1055  if (abs(matv(2,m))>1e-16) then
1056  if (verbose>1) print*,"cond2p k=dimk m=",m
1057  if (neumann==1) then
1058  if (verbose>1) print*,"boundary conditions not compatible"
1059  condlimord2_kk=1
1060  endif
1061  else
1062  if (verbose>1) print*,"boundary conditions not supported"
1063  condlimord2_kk=1
1064  endif
1065  endif
1066  enddo
1067  if (interpretatif==1) condlimord2_kk=1
1068 
1069  if (condlimord2_kk.ne.0) return
1070 
1071  do m=1,dimm
1072  if (mode(m)==1) cycle
1073  if (abs(matu(1,m))>1e-16) then
1074  aux=matu(1,m)
1075  matu(1,m) = - matw(1,m)/aux
1076  matv(1,m) = - matv(1,m)/aux
1077  matw(1,m) = v1(1,m) /aux
1078  do i=1,dimm
1079  matb(dimk,m,i) = matb(dimk,m,i) + matv(1,m)*matc(dimk,m,i)
1080  mata(dimk,m,i) = mata(dimk,m,i) + matu(1,m)*matc(dimk,m,i)
1081  mats(dimk,i) = mats(dimk,i) + sca_f*matw(1,m)*matc(dimk,m,i)
1082  enddo
1083  !matc(dimk,m,:) = 0.0
1084  else
1085  matx(1,m) = -matw(1,m)/matv(1,m)
1086  maty(1,m) = v1(1,m) /matv(1,m)
1087  mata(dimk,:,m) = 0.0_8
1088  matb(dimk,:,m) = 0.0_8
1089  mats(dimk,m) = 0.0_8
1090  do i=1,dimm
1091  matb(dimk-1,m,i) = matb(dimk-1,m,i) + matx(1,m)*matc(dimk-1,m,i)
1092  mats(dimk-1,i) = mats(dimk-1,i) + sca_f*maty(1,m)*matc(dimk-1,m,i)
1093  enddo
1094  matc(dimk-1,m,:) = 0.0_8
1095  endif
1096  if (abs(matu(2,m))>1e-16) then
1097  aux=matu(2,m)
1098  matu(2,m) = - matw(2,m)/aux
1099  matv(2,m) = - matv(2,m)/aux
1100  matw(2,m) = v1(2,m) /aux
1101  do i=1,dimm
1102  matbp(dimk,m,i) = matbp(dimk,m,i) + matv(2,m)*matcp(dimk,m,i)
1103  matap(dimk,m,i) = matap(dimk,m,i) + matu(2,m)*matcp(dimk,m,i)
1104  mats(dimk,i) = mats(dimk,i) + (1.0_8-sca_f)*matw(2,m)*matcp(dimk,m,i)
1105  enddo
1106  !matcp(dimk,m,:) = 0.0
1107  else
1108  matx(2,m) = -matw(2,m)/matv(2,m)
1109  maty(2,m) = v1(2,m) /matv(2,m)
1110  matap(dimk,:,m) = 0.0_8
1111  matap(dimk,m,m) = matx(2,m)/(1.0_8 - sca_f)
1112  matbp(dimk,:,m) = 0.0_8
1113  mats(dimk,m) = maty(2,m)
1114  do i=1,dimm
1115  matbp(dimk-1,m,i) = matbp(dimk-1,m,i) + matx(2,m)*matcp(dimk-1,m,i)
1116  mats(dimk-1,i) = mats(dimk-1,i) + (1.0_8-sca_f)*maty(2,m)*matcp(dimk-1,m,i)
1117  enddo
1118  matcp(dimk-1,m,:) = 0.0_8
1119  endif
1120  enddo
1121 
1122 end function condlimord2_kk
1123 !**************************************************************************************
1124 subroutine condlimord1_k1
1125  !***********************
1126  !boundary conditions k=1
1127  !***********************
1128 
1129  use cos_transportdata
1130  implicit none
1131 
1132  integer :: m,i
1133  real(kind=8) :: aux
1134 
1135  matu = t0(:,:,2)
1136  matv = t0(:,:,1)*dx - t0(:,:,2)
1137  maty = v0(:,:)*dx
1138  print*,"dx",dx,t0(2,1,1), t0(2,1,2)
1139  print*,"matv",matv(2,1)
1140 
1141  do m=1,dimm
1142  if (mode(m)==1) cycle
1143  matx(1,m) = -matu(1,m)/matv(1,m)
1144  maty(1,m) = maty(1,m)/matv(1,m)
1145  matb(1,:,m) = 0.0_8
1146  matc(1,:,m) = 0.0_8
1147  mats(1,m) = 0.0_8
1148  do i=1,dimm
1149  matb(2,m,i) = matb(2,m,i) + matx(1,m)*mata(2,m,i)
1150  mats(2,i) = mats(2,i) + sca_f*maty(1,m)*mata(2,m,i)
1151  enddo
1152  mata(2,m,:) = 0.0_8
1153  if (abs(matv(2,m))>1e-16) then
1154  matx(2,m) = -matu(2,m)/matv(2,m)
1155  maty(2,m) = maty(2,m)/matv(2,m)
1156  matbp(1,:,m) = 0.0_8
1157  matcp(1,:,m) = 0.0_8
1158  matcp(1,m,m) = matx(2,m)/(1.0_8 - sca_f)
1159  mats(1,m) = maty(2,m)
1160  do i=1,dimm
1161  matbp(2,m,i) = matbp(2,m,i) + matx(2,m)*matap(2,m,i)
1162  mats(2,i) = mats(2,i) + (1.0_8-sca_f)*maty(2,m)*matap(2,m,i)
1163  enddo
1164  matap(2,m,:) = 0.0_8
1165  else
1166  print*,"Boundaries conditions error",matv(2,m)
1167  endif
1168  enddo
1169 end subroutine condlimord1_k1
1170 !**********************************************************************************************
1171 subroutine condlimord1_kk
1172  !**************************
1173  !boundary conditions k=dimk
1174  !**************************
1175 
1176  use cos_transportdata
1177  implicit none
1178 
1179  integer :: m,i
1180  real(kind=8) :: aux
1181 
1182  matv = t1(:,:,2) + dx*t1(:,:,1)
1183  matw = - t1(:,:,2)
1184  maty = v1(:,:)*dx
1185 
1186  do m=1,dimm
1187  if (mode(m)==1) cycle
1188  matx(1,m) = -matw(1,m)/matv(1,m)
1189  maty(1,m) = maty(1,m)/matv(1,m)
1190  mata(dimk,:,m) = 0.0_8
1191  matb(dimk,:,m) = 0.0_8
1192  mats(dimk,m) = 0.0_8
1193  do i=1,dimm
1194  matb(dimk-1,m,i) = matb(dimk-1,m,i) + matx(1,m)*matc(dimk-1,m,i)
1195  mats(dimk-1,i) = mats(dimk-1,i) + sca_f*maty(1,m)*matc(dimk-1,m,i)
1196  enddo
1197  matc(dimk-1,m,:) = 0.0_8
1198  if (abs(matv(2,m))>1e-16) then
1199  matx(2,m) = -matw(2,m)/matv(2,m)
1200  maty(2,m) = maty(2,m)/matv(2,m)
1201  matap(dimk,:,m) = 0.0_8
1202  matap(dimk,m,m) = matx(2,m)/(1.0_8 - sca_f)
1203  matbp(dimk,:,m) = 0.0_8
1204  mats(dimk,m) = maty(2,m)
1205  do i=1,dimm
1206  matbp(dimk-1,m,i) = matbp(dimk-1,m,i) + matx(2,m)*matcp(dimk-1,m,i)
1207  mats(dimk-1,i) = mats(dimk-1,i) + (1.0_8-sca_f)*maty(2,m)*matcp(dimk-1,m,i)
1208  enddo
1209  matcp(dimk-1,m,:) = 0.0_8
1210  else
1211  print*,"Boundaries conditions error",matv(2,m)
1212  endif
1213  enddo
1214 
1215 end subroutine condlimord1_kk
1216 
1217 #endif
build_matcoo
subroutine build_matcoo(matp, option, k, counter, counteri, counterj)
Definition: solution10.f90:774
cos_zconversion
subroutine cos_zconversion(e, nbrho)
These subroutines correct the central value (if not finite value)
Definition: solution10.f90:342
derivn10
subroutine derivn10(N, X, Y, DY1)
These subroutines calculate first and second derivatives, DY1 and DY2, of function Y respect to argum...
Definition: solution10.f90:306
cos_precision
types
Definition: cos_precision.f90:9
type_transport_solver_numerics::solver_io
The SOLVER_IO type includes variables that are used input and output to the nuermical solver of the t...
Definition: type_transport_solver_numerics.f90:63
cos_pde1solver
subroutine cos_pde1solver(f_in, fp_in, mata_in, matb_in, matc_in, matd_in, matap_in, matbp_in, matcp_in, matdp_in, dimk_in, dimn_in, dimm_in, dx_in, dt_in, sca_f_in, T0_in, T1_in, V0_in, V1_in, mode_in, dfpdx_in, dfpdx2_in)
Definition: solution10.f90:407
condlimord2_k1
integer function condlimord2_k1(interpretatif)
Definition: solution10.f90:886
cos_transportdata
Definition: solution10.f90:365
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
condlimord1_k1
subroutine condlimord1_k1
Definition: solution10.f90:1124
condlimord2_kk
integer function condlimord2_kk(interpretatif)
Definition: solution10.f90:1006
solvepde1d
subroutine solvepde1d
Definition: solution10.f90:814
type_transport_solver_numerics
The module declares types of variables used by numerical solvers.
Definition: type_transport_solver_numerics.f90:12
condlimord1_kk
subroutine condlimord1_kk
Definition: solution10.f90:1171
Generated by   doxygen 1.8.5