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