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