1!-------------------------------------------------------------------------------
2SUBROUTINE ZirkaTest( Model,Solver,dt,TransientSimulation ) ! {{{
3!-------------------------------------------------------------------------------
4  USE DefUtils
5  USE ZirkaUtils
6
7  IMPLICIT NONE
8!------------------------------------------------------------------------------
9  TYPE(Solver_t) :: Solver
10  TYPE(Model_t) :: Model
11
12  REAL(KIND=dp) :: dt
13  LOGICAL :: TransientSimulation
14  integer :: unit, testnum
15  TYPE(ValueList_t), POINTER :: params
16  type(Variable_t), POINTER :: variable
17  logical :: found
18!------------------------------------------------------------------------------
19! Local variables
20!------------------------------------------------------------------------------
21
22  testnum = ListGetInteger(GetSolverparams(Solver), 'test number', found)
23  if (.not. found) testnum = 6
24
25  call info('ZirkaTest','Running zirkatest')
26  select case(testnum)
27  case (6); call test6()
28  case (5); call test5()
29  case (4); call test4()
30  case (3); call test3()
31  end select
32
33  OPEN (newunit=unit, FILE='TEST.PASSED')
34  write(unit, '(I0)') 1
35
36  close(unit)
37
38contains
39
40subroutine test6() ! {{{
41  use zirka
42  use defutils
43  type(MasterCurve_t) :: mc, mc2
44  real(kind=dp), allocatable :: BHasc(:,:), BHsingle(:,:), Bsingle(:), Hsingle(:)
45  real(kind=dp), allocatable :: rps(:)
46  real(kind=dp), allocatable :: H(:), dH(:), Hnobuf(:)
47  real(kind=dp), pointer :: rps_array(:,:)
48  integer  :: n, m
49  integer, parameter :: nr = 20
50  type(ZirkaABC_t), POINTER :: ABCParams
51  allocate(ABCParams )
52
53  BHasc = readBH('../HB/BH_asc_real')
54  BHsingle = readBH('../HB/BH_single_real')
55
56  allocate(H(size(BHasc,1)), dH(size(BHasc,1)), Hnobuf(size(BHasc,1)))
57
58  allocate(rps(nr))
59  rps = [(2.0*(-0.5_dp)**n, n=0,nr-1)]
60
61  mc = MasterCurve_t(InitSplineLoop(BHasc, BHsingle), &
62      ABCparams, n_cachesubsample = 2)
63
64  call mc % drive(1.8_dp)
65
66  DO n = 1,nr-1
67    CALL mc % drive(rps(n), cache=.false.)
68  END DO
69  CALL mc % drive(rps(nr), cache=.true.)
70
71  rps_array => null()
72  m = 250
73  BLOCK
74    integer :: i_r
75    rps_array => ListGetConstRealArray(GetSolverParams(), 'rps', found)
76    if (.not. found) then
77      rps = [0.0_dp, 1.8_dp, -1.0_dp, 0.5_dp, 0.25_dp, 0.4_dp, -0.4_dp, -0.8_dp, 0.9_dp]
78    else
79      rps = rps_array(:,1)
80    end if
81
82    ! rps = -rps
83    ! rps = [0.0_dp, -1.84102067_dp, 1.3_dp, -0.5_dp]
84    do i_r = 2, ubound(rps,1)
85      call mc % drive(rps(i_r-1), cache = .true.)
86        associate (B_low => rps(i_r-1),  B_up => rps(i_r))
87          call print_interval(mc, B_low, B_up, m)
88        end associate
89    end do
90  END BLOCK
91
92end subroutine ! }}}
93
94subroutine print_interval(mc, B_low, B_up, m, Hexact) ! {{{
95  type(MasterCurve_t) :: mc
96  real(kind=dp) :: B_low, B_up, H, Hnobuf, dH
97  real(kind=dp), intent(in), optional :: Hexact
98  integer :: n, m
99  !-------------------------------------------------------------------------------
100  DO n = 0, m-1
101    associate( Bhere => (B_low + n*(B_up-B_low)/m) )
102      H = mc % eval(Bhere, cached = .true., dhdb=dH)
103      Hnobuf = mc % eval(Bhere, cached = .false.)
104
105      if (present(Hexact)) then
106        print *, Bhere, H, dH, Hnobuf, Hexact
107      else
108        print *, Bhere, H, dH, Hnobuf
109      end if
110
111    end associate
112  end do
113end subroutine ! }}}
114
115subroutine test5() ! {{{
116  use zirka
117  use defutils
118  type(MasterCurve_t) :: mc, mc2
119  real(kind=dp), allocatable :: BHasc(:,:), BHsingle(:,:), Bsingle(:), Hsingle(:)
120  real(kind=dp), allocatable :: rps(:)
121  real(kind=dp), allocatable :: H(:), dH(:), Hnobuf(:)
122  integer  :: n, m
123  integer, parameter :: nr = 20
124  type(ZirkaABC_t), POINTER :: ABCParams
125  allocate(ABCParams )
126
127  BHasc = readBH('../HB/BH_asc_real')
128  BHsingle = readBH('../HB/BH_single_real')
129
130  allocate(H(size(BHasc,1)), dH(size(BHasc,1)), Hnobuf(size(BHasc,1)))
131
132  allocate(rps(nr))
133  rps = [(2.0*(-0.5_dp)**n, n=0,nr-1)]
134
135  mc = MasterCurve_t(InitSplineLoop(BHasc, BHsingle), &
136      ABCparams)
137
138  call mc % drive(1.8_dp)
139
140  DO n = 1,nr-1
141    CALL mc % drive(rps(n), cache=.false.)
142  END DO
143  CALL mc % drive(rps(nr), cache=.true.)
144
145
146  m = 250
147  associate (B_low => -0.0_dp, B_up => 1.9_dp )
148    DO n = 1,m
149      associate( Bhere => (B_low + n*(B_up-B_low)/m) )
150        H(n) = mc % eval(Bhere, cached = .true., dhdb=dH(n))
151        Hnobuf(n) = mc % eval(Bhere)
152        print *, Bhere, H(n), dH(n), Hnobuf(n)
153      end associate
154    END DO
155  end associate
156
157end subroutine ! }}}
158!-------------------------------------------------------------------------------
159subroutine test4() ! {{{
160  use zirka
161  use defutils
162  type(MasterCurve_t) :: mc, mc2
163  real(kind=dp), allocatable :: BHasc(:,:), BHsingle(:,:), Bsingle(:), Hsingle(:)
164  real(kind=dp), allocatable :: rps(:)
165  real(kind=dp) :: H
166  integer  :: n, m
167  integer, parameter :: nr = 20
168  type(ZirkaABC_t), POINTER :: ABCParams
169  real(kind=dp), pointer :: rps_array(:,:)
170  allocate(ABCParams )
171
172  BHasc = readBH('../HB/BH_asc_real')
173  BHsingle = readBH('../HB/BH_single_real')
174
175
176  allocate(rps(nr))
177  rps = [(2.0*(-0.5_dp)**n, n=0,nr-1)]
178
179  mc = MasterCurve_t(InitSplineLoop(BHasc, BHsingle), &
180      ABCparams)
181
182  rps_array => null()
183  rps_array => ListGetConstRealArray(GetSolverParams(), 'rps_4', found)
184
185  do n=1,size(rps_array,1)
186    call mc % drive(rps_array(n,1), cache = .false.)
187    call mc % printme()
188    print *, 'B = ', rps_array(n,1)
189    H = mc % eval (rps_array(n,1), cached = .false.)
190    print *, 'H = ',H
191    ! call printeval(mc, rps_array(n,1))
192  end do
193  ! call mc % drive(0.1_dp, cache = .false.)
194  ! call mc % printme()
195  ! call mc % drive(-1.8_dp, cache = .false.)
196  ! call mc % printme()
197  ! call mc % drive(0.8_dp, cache = .false.)
198  ! call mc % printme()
199  ! call mc % drive(-1.2_dp, cache = .false.)
200  ! call mc % printme()
201
202  ! deallocate(rps_array)
203
204end subroutine ! }}}
205subroutine test3() ! {{{
206  use zirka
207  use defutils
208  type(MasterCurve_t) :: mc, mc2
209  real(kind=dp), pointer :: BHasc(:,:), BHsingle(:,:)
210  real(kind=dp), pointer :: rps(:), H(:)
211  real(kind=dp), pointer :: rps_array(:,:)
212  real(kind=dp), allocatable, target :: scalar_array(:,:)
213
214  integer  :: n, m
215  integer :: nr
216  type(ZirkaABC_t), POINTER :: ABCParams
217  allocate(ABCParams )
218
219  call GetConstRealArray(model % materials(2) % values, BHasc, 'Ascending BH curve', found)
220  if(.not. Found) call fatal('ZirkaTest', 'Ascending BH curve not found')
221  call GetConstRealArray(model % materials(2) % values, BHSingle, 'Single valued BH curve', found)
222  if(.not. Found) call fatal('ZirkaTest', 'Single valued BH curve not found')
223
224  call info('ZirkaTest','Running zirkatest')
225  rps_array => ListGetConstRealArray(GetSolverParams(), 'zirka init sequence', found)
226  rps(1:size(rps_array,1)) => rps_array(:,1)
227
228  mc = MasterCurve_t(InitSplineLoop(BHasc, BHsingle), &
229      ABCparams, n_cachesubsample = 2)
230
231  nr = size(rps,1)
232
233  DO n = 1,nr-1
234    CALL mc % drive(rps(n), cache=.false.)
235  END DO
236  CALL mc % drive(rps(nr), cache=.true.)
237
238  rps => null()
239
240  m = 1
241  BLOCK
242    integer :: i_r
243    ! rps_array => ListGetConstRealArray(GetSolverParams(), 'rps_3', found)
244    scalar_array = readdat("scalars_a.dat", 7, 9)
245  call info('ZirkaTest','Running zirkatest')
246    rps(1:size(scalar_array,1)) => scalar_array(:,3)
247    H(1:size(scalar_array,1)) => scalar_array(:,4)
248
249    do i_r = 2, ubound(rps,1)
250      call mc % drive(rps(i_r-1), cache = .true.)
251        associate (B_low => rps(i_r-1),  B_up => rps(i_r))
252          call print_interval(mc, B_low, B_up, m, H(i_r-1))
253        end associate
254    end do
255    i_r = ubound(rps,1)
256    call mc % drive(rps(i_r), cache = .true.)
257    associate (B_low => rps(i_r),  B_up => rps(i_r))
258      call print_interval(mc, B_low, B_up, m, H(i_r))
259    end associate
260    variable => VariableGet(solver % variable, "testvar")
261    variable % values = abs(H(i_r)-mc % eval(rps(i_r), cached = .true.))/abs(H(i_r))
262    variable % norm = abs(H(i_r)-mc % eval(rps(i_r), cached = .true.))/abs(H(i_r))
263  END BLOCK
264
265end subroutine ! }}}
266function readBH(f) result(BH) ! {{{
267  real(kind=dp), allocatable :: BH(:,:)
268  character(len=*) :: f
269  real*8 :: lenf
270  integer :: iu, len, n
271
272  open(newunit=iu, file=f, action='read', status='old')
273
274  read(iu, *) len
275
276  allocate(BH(1:len,2))
277
278  do n = 1,len
279    read (iu,*) BH(n,1:2)
280  end do
281  close(iu)
282
283end function ! }}}
284
285function readdat(f, m, n) result(table) ! {{{
286  real(kind=dp), allocatable :: table(:,:)
287  character(len=*) :: f
288  real*8 :: lenf
289  integer :: iu, m, n, k
290
291  open(newunit=iu, file=f, action='read', status='old')
292
293  allocate(table(m,n))
294
295  do k = 1,m
296    read (iu,*) table(k,1:n)
297  end do
298  close(iu)
299
300end function ! }}}
301END SUBROUTINE ZirkaTest ! }}}
302!-------------------------------------------------------------------------------
303