1! (This is a demo program for the alternative workflow PSD, or 2! Preconditioner and Solver Decoupled, contributed by Neil Hodge. 3! To execute this program, it must be built with an MPI library, 4! and the Fortran 90 interface and the SA-AMG preconditioner must 5! be enabled.) 6 7PROGRAM lis_driver 8 9 IMPLICIT NONE 10 11! per the Intel fortran user forums, including a *.h file 12! via the preprocessor mechanism works properly, versus using 13! the fortran include mechanism 14 15#include "lisf.h" 16 17 INTERFACE 18 SUBROUTINE plotxy(plotdata,filename,datalength, & 19 grid,limits,title1,xlabel,ylabel,series_labels, & 20 legend,lines,points,colors,widths,output_animate) 21 REAL(kind(0.d0)),INTENT(IN) :: plotdata(:,:) 22 CHARACTER(LEN=*),INTENT(IN) :: filename 23 INTEGER,INTENT(IN),OPTIONAL :: datalength(:) 24 INTEGER,INTENT(IN),OPTIONAL :: grid 25 REAL(kind(0.d0)),INTENT(IN),OPTIONAL :: limits(:,:) 26 CHARACTER(LEN=*),INTENT(IN),OPTIONAL :: title1 27 CHARACTER(LEN=*),INTENT(IN),OPTIONAL :: xlabel,ylabel 28 CHARACTER(LEN=*),INTENT(IN),OPTIONAL :: series_labels(:) 29 INTEGER,INTENT(IN),OPTIONAL :: legend 30 INTEGER,INTENT(IN),OPTIONAL :: lines(:) 31 INTEGER,INTENT(IN),OPTIONAL :: points(:) 32 INTEGER,INTENT(IN),OPTIONAL :: colors(:) 33 INTEGER,INTENT(IN),OPTIONAL :: widths(:) 34 INTEGER,INTENT(IN),OPTIONAL :: output_animate 35 END SUBROUTINE 36 END INTERFACE 37 38 !====================================================================== 39 ! parameters 40 !====================================================================== 41 INTEGER,PARAMETER :: singI=kind(0) 42 INTEGER,PARAMETER :: fullR=kind(0.d0) 43 44 REAL(fullR),PARAMETER :: ZERO = 0.0_fullR 45 REAL(fullR),PARAMETER :: ROOT_EPS = 0.0000000596046448_fullR ! 2^(-24) = 5.96046448E-8 46 REAL(fullR),PARAMETER :: ONE = 1.0_fullR 47 REAL(fullR),PARAMETER :: TWO = 2.0_fullR 48 REAL(fullR),PARAMETER :: THREE = 3.0_fullR 49 REAL(fullR),PARAMETER :: FOUR = 4.0_fullR 50 REAL(fullR),PARAMETER :: FIVE = 5.0_fullR 51 REAL(fullR),PARAMETER :: SIX = 6.0_fullR 52 REAL(fullR),PARAMETER :: SEVEN = 7.0_fullR 53 REAL(fullR),PARAMETER :: EIGHT = 8.0_fullR 54 REAL(fullR),PARAMETER :: NINE = 9.0_fullR 55 REAL(fullR),PARAMETER :: TEN = 10.0_fullR 56 REAL(fullR),PARAMETER :: HALF = 0.500000000000000_fullR 57 REAL(fullR),PARAMETER :: THIRD = 0.333333333333333_fullR 58 REAL(fullR),PARAMETER :: FOURTH = 0.250000000000000_fullR 59 REAL(fullR),PARAMETER :: TENTH = 0.100000000000000_fullR 60 REAL(fullR),PARAMETER :: PI = 3.141592653589793_fullR 61 62 !====================================================================== 63 ! regular old variables 64 !====================================================================== 65 INTEGER(singI) :: ierr,MPIsize,rank 66 LIS_INTEGER :: i,j,gn,n,is,ie,iter 67 LIS_INTEGER,POINTER :: proc_n(:) 68 LIS_VECTOR :: b,x 69 LIS_MATRIX :: A 70 LIS_PRECON :: precon 71 LIS_SOLVER :: solver 72 CHARACTER(LEN=1024) :: options 73 REAL(fullR),POINTER :: RHS(:),LHS(:,:) 74 LOGICAL,POINTER :: LHSassem(:,:) 75 INTEGER(singI),POINTER :: HomeProc(:) 76 LOGICAL :: LinSysDefFlag 77 REAL(fullR) :: ResNorm0,ResNorm,AbsResTol,RelResTol 78 REAL(fullR) :: IncNorm0,IncNorm,AbsIncTol,RelIncTol 79 INTEGER(singI) :: IterCount,MaxIterCount 80 INTEGER(singI) :: nel,nnp,inp,DOF1,DOF2,iProc 81 REAL(fullR) :: h 82 REAL(fullR),POINTER :: xhat(:),That(:),dT(:),plotdata(:,:) 83 INTEGER(singI) :: PlotCount 84 INTEGER(singI),POINTER :: lps(:),widths(:) 85 CHARACTER(LEN=1024) :: title,PlotFilename 86 INTEGER(singI) :: NonlinLHSUpdateFreq,NumberLHSUpdates 87 88 proc_n=>NULL() 89 LHS=>NULL() 90 RHS=>NULL() 91 LHSassem=>NULL() 92 HomeProc=>NULL() 93 xhat=>NULL() 94 That=>NULL() 95 dT=>NULL() 96 plotdata=>NULL() 97 lps=>NULL() 98 widths=>NULL() 99 100 101 !====================================================================== 102 ! proceed 103 !====================================================================== 104 CALL MPI_Init(ierr) 105! CALL ErrorCheck("MPI_Init",ierr) 106 107 CALL MPI_Comm_size(MPI_COMM_WORLD,MPIsize,ierr) 108! CALL ErrorCheck("MPI_Comm_rank",ierr) 109 CALL MPI_Comm_rank(MPI_COMM_WORLD,rank,ierr) 110! CALL ErrorCheck("MPI_Comm_rank",ierr) 111 112 IF ( MPIsize.LT.4 ) THEN 113 IF ( rank.EQ.0 ) THEN 114 WRITE(UNIT=*,FMT='(a)') "number of processes must be 4 or more" 115 CALL lis_finalize(ierr) 116 ENDIF 117 STOP 118 ENDIF 119 120 ! do it! 121 CALL lis_initialize(ierr) 122! CALL ErrorCheck("lis_initialize",ierr) 123 124 ! mesh definition 125 nel=100 126 nnp=nel+1 127 128 ! allocate data 129 ALLOCATE(xhat(nnp)) 130 ALLOCATE(That(nnp)) 131 ALLOCATE(RHS(nnp)) 132 ALLOCATE(LHS(nnp,nnp)) 133 ALLOCATE(LHSassem(nnp,nnp)) 134 IF (rank==0) THEN 135 ALLOCATE(plotdata(nnp,2)) 136 ALLOCATE(lps(1)) 137 lps=1 138 ALLOCATE(widths(1)) 139 widths=3 140 END IF 141 142 ! parallelism data 143 ALLOCATE(proc_n(0:MPIsize-1)) 144 proc_n=0 145 !====================================== 146 ! DOFs on all processes 147 !====================================== 148! proc_n(0)=INT(REAL(nnp,fullR)/TEN,singI) 149! proc_n(1)=2*proc_n(0) 150! proc_n(2)=3*proc_n(0) 151! proc_n(3)=nnp-SUM(proc_n(0:2)) 152 !====================================== 153 ! alternate cases, testing out various versions of n_local=0 154 !====================================== 155! proc_n(0)=0 156! proc_n(1)=INT(nnp/4,singI) 157! proc_n(2)=proc_n(1) 158! proc_n(3)=nnp-SUM(proc_n(0:2)) 159 !====================================== 160 proc_n(0)=0 161 proc_n(1)=INT(REAL(nnp,fullR)/THREE,singI) 162 proc_n(2)=0 163 proc_n(3)=nnp-SUM(proc_n(0:2)) 164 ! output DOFs per process 165 IF (rank==0) THEN 166 DO i=0,MPIsize-1 167 WRITE(UNIT=*,FMT='(a,i0,a,i0)') "proc_n(",i,")=",proc_n(i) 168 END DO 169 WRITE(UNIT=*,FMT='(a)') "" 170 END IF 171 172 ! initial condition 173 h=ONE/REAL(nel,fullR) 174 DO inp=1,nnp 175 xhat(inp)=REAL((inp-1),fullR)*h 176 That(inp)=202.0+202.0*COS(FOUR*THIRD*SEVEN*xhat(inp)) 177 IF (rank==0) plotdata(inp,1)=xhat(inp) 178 END DO 179 ! Dirichlet BC 180 That(1)=101.0 181 That(nnp)=808.0 182 ! initial state plot 183 IF (rank==0) THEN 184 DO inp=1,nnp 185 plotdata(inp,2)=That(inp) 186 END DO 187 WRITE(UNIT=title,FMT='(a,i0)') "nonlinear iteration ",0 188 PlotCount=1 189 WRITE(UNIT=PlotFilename,FMT='(a,i0)') "blah",PlotCount 190 CALL plotxy(plotdata=plotdata,filename=TRIM(PlotFilename), & 191 lines=lps,points=lps,widths=widths, & 192 title1=TRIM(title),output_animate=1) 193 ! write initial state twice to deal with ffmpeg bug . . . 194 PlotCount=PlotCount+1 195 WRITE(UNIT=PlotFilename,FMT='(a,i0)') "blah",PlotCount 196 CALL plotxy(plotdata=plotdata,filename=TRIM(PlotFilename), & 197 lines=lps,points=lps,widths=widths, & 198 title1=TRIM(title),output_animate=1) 199 END IF 200 201 ! convergence and iteration-related criteria 202 AbsResTol=1.0d-06 203 RelResTol=1.0d-06 204 AbsIncTol=1.0d-05 205 RelIncTol=1.0d-05 206 MaxIterCount=50 207 NonlinLHSUpdateFreq=3 208 NumberLHSUpdates=0 209 210 ! update linear system 211 ! do this outside the nonlinear loop, just for the first iteration 212 CALL UpdateLinearSystemThermal(nel,xhat,That,RHS,LHS,LHSassem) 213 IncNorm0=ZERO 214 ResNorm0=ZERO 215 gn=SUM(proc_n) 216 DO i=1,gn 217 ResNorm0=ResNorm0+RHS(i)**2 218 END DO 219 ResNorm0=SQRT(ResNorm0) 220 IF (rank==0) THEN 221 WRITE(UNIT=*,FMT='(a,es22.15)') "initial residual norm:",ResNorm0 222 WRITE(UNIT=*,FMT='(a,es22.15)') "absolute residual tolerance:",AbsResTol 223 WRITE(UNIT=*,FMT='(a,es22.15)') "relative residual tolerance:",RelResTol 224 WRITE(UNIT=*,FMT='(a,es22.15)') "absolute increment tolerance:",AbsIncTol 225 WRITE(UNIT=*,FMT='(a,es22.15)') "relative increment tolerance:",RelIncTol 226 WRITE(UNIT=*,FMT='(a,i0)') "MaxIterCount: ",MaxIterCount 227 WRITE(UNIT=*,FMT='(a,i0)') "NonlinLHSUpdateFreq: ",NonlinLHSUpdateFreq 228 WRITE(UNIT=*,FMT='(a)') "" 229 END IF 230 ResNorm=ResNorm0 231 232 ! HomeProc indicates the sender in Bcast below 233 ALLOCATE(HomeProc(gn)) 234 DOF2=0 235 DO iProc=0,(MPIsize-1) 236 DOF1=DOF2+1 237 DOF2=SUM(proc_n(0:iProc)) 238 HomeProc(DOF1:DOF2)=iProc 239 END DO 240 ALLOCATE(dT(gn)) 241 242 ! calculate solution 243 LinSysDefFlag=.FALSE. 244 IterCount=0 245 DO 246 IterCount=IterCount+1 247 248 CALL MPI_Barrier(MPI_COMM_WORLD,ierr) 249 CALL pauser(HALF) 250 IF (rank==0) THEN 251 WRITE(UNIT=*,FMT='(a)') "" 252 WRITE(UNIT=*,FMT='(a)') "======================================================================" 253 WRITE(UNIT=*,FMT='(a,i0)') "Start of nonlinear iteration ",IterCount 254 WRITE(UNIT=*,FMT='(a)') "======================================================================" 255 END IF 256 257 ! define once, and assemble 258 IF (.NOT.LinSysDefFlag) THEN 259 CALL lis_matrix_create(MPI_COMM_WORLD,A,ierr) 260! CALL ErrorCheck("lis_matrix_set_size",ierr) 261 CALL lis_matrix_set_size(A,proc_n(rank),0,ierr) 262! CALL ErrorCheck("lis_matrix_set_size",ierr) 263 CALL lis_matrix_get_size(A,n,gn,ierr) 264! CALL ErrorCheck("lis_matrix_get_size",ierr) 265 CALL lis_matrix_get_range(A,is,ie,ierr) 266! CALL ErrorCheck("lis_matrix_get_range",ierr) 267 END IF 268 269 IF (IterCount==1) THEN 270 DO i=is,ie-1 271 DO j=1,gn 272 IF (.NOT.LHSassem(i,j)) CYCLE 273 CALL lis_matrix_set_value(LIS_INS_VALUE,i,j,LHS(i,j),A,ierr) 274! CALL ErrorCheck("lis_matrix_set_value",ierr) 275 END DO 276 END DO 277 ELSE 278 IF (NonlinLHSUpdateFreq>0 .AND. MOD(IterCount,NonlinLHSUpdateFreq)==0) THEN 279 NumberLHSUpdates=NumberLHSUpdates+1 280 DO i=is,ie-1 281 DO j=1,gn 282 IF (.NOT.LHSassem(i,j)) CYCLE 283 CALL lis_matrix_psd_set_value(LIS_INS_VALUE,i,j,LHS(i,j),A,ierr) 284! CALL ErrorCheck("lis_matrix_psd_set_value",ierr) 285 END DO 286 END DO 287 CALL lis_matrix_psd_reset_scale(A,ierr) 288! CALL ErrorCheck("lis_matrix_psd_reset_scale",ierr) 289 END IF 290 END IF 291 292 IF (.NOT.LinSysDefFlag) THEN 293 CALL lis_matrix_set_type(A,LIS_MATRIX_CSR,ierr) 294! CALL ErrorCheck("lis_matrix_set_type",ierr) 295 CALL lis_matrix_assemble(A,ierr) 296! CALL ErrorCheck("lis_matrix_assemble",ierr) 297 CALL lis_vector_duplicate(A,b,ierr) 298! CALL ErrorCheck("lis_vector_duplicate",ierr) 299 CALL lis_vector_duplicate(A,x,ierr) 300! CALL ErrorCheck("lis_vector_duplicate",ierr) 301 END IF 302 303 ! update the RHS _every_ nonlinear iteration 304 DO i=is,ie-1 305 CALL lis_vector_set_value(LIS_INS_VALUE,i,RHS(i),b,ierr) 306! CALL ErrorCheck("lis_vector_set_value",ierr) 307 END DO 308 CALL lis_vector_psd_reset_scale(b,ierr) 309! CALL ErrorCheck("lis_vector_psd_reset_scale",ierr) 310 311 IF (.NOT.LinSysDefFlag) THEN 312 CALL lis_solver_create(solver,ierr) 313! CALL ErrorCheck("lis_solver_create",ierr) 314! WRITE(UNIT=options,FMT='(a)') "-p none -i gmres -print out -scale none" 315! WRITE(UNIT=options,FMT='(a)') "-p none -i gmres -print out -scale jacobi" 316! WRITE(UNIT=options,FMT='(a)') "-p ilu -i gmres -print out -scale none" 317 WRITE(UNIT=options,FMT='(a)') "-p saamg -saamg_unsym true -i gmres -print out -scale none" 318 WRITE(UNIT=options,FMT='(a,a,es22.15)') TRIM(options)," -tol ", & 319 & MAX(MIN(AbsResTol/ResNorm,ONE),MIN(RelResTol*ResNorm0/ResNorm,ONE))/TEN 320 IF (rank==0) WRITE(UNIT=*,FMT='(a,a,a)') "solver options: """,TRIM(options),"""" 321 CALL lis_solver_set_option(TRIM(options),solver,ierr) 322! CALL ErrorCheck("lis_solver_set_option",ierr) 323 ! be sure to create and set options for solver before 324 ! calling precon_create . . . 325 ! A has to be set in "solver" _before_ precon_create is called . . . 326 CALL lis_solver_set_matrix(A,solver,ierr) 327! CALL ErrorCheck("lis_solver_set_matrix",ierr) 328 CALL lis_precon_psd_create(solver,precon,ierr) 329! CALL ErrorCheck("lis_precon_psd_create",ierr) 330 CALL lis_precon_psd_update(solver,precon,ierr) 331! CALL ErrorCheck("lis_precon_psd_update",ierr) 332 LinSysDefFlag=.TRUE. 333 ELSE 334 WRITE(UNIT=options,FMT='(a,es22.15)') "-tol ",MAX(MIN(AbsResTol/ResNorm,ONE),MIN(RelResTol*ResNorm0/ResNorm,ONE))/TEN 335 IF (rank==0) WRITE(UNIT=*,FMT='(a,a,a)') "solver options: """,TRIM(options),"""" 336 CALL lis_solver_set_option(TRIM(options),solver,ierr) 337! CALL ErrorCheck("lis_solver_set_option",ierr) 338 END IF 339 340 IF (NonlinLHSUpdateFreq>0 .AND. MOD(IterCount,NonlinLHSUpdateFreq)==0) THEN 341 IF (rank==0) WRITE(UNIT=*,FMT='(a)') "updating preconditioner now . . ." 342 CALL lis_precon_psd_update(solver,precon,ierr) 343! CALL ErrorCheck("lis_precon_psd_update",ierr) 344 END IF 345 CALL lis_solve_kernel(A,b,x,solver,precon,ierr) 346! CALL ErrorCheck("lis_solve_kernel",ierr) 347 CALL lis_solver_get_iter(solver,iter,ierr) 348! CALL ErrorCheck("lis_solver_get_iter",ierr) 349 IF (rank==0) WRITE(UNIT=*,FMT='(a,i0)') "number of iterations = ",iter 350 351 ResNorm=ZERO 352 DO i=is,ie-1 353 CALL lis_vector_get_value(x,i,dT(i),ierr) 354! CALL ErrorCheck("lis_vector_get_value",ierr) 355 dT(i)=-dT(i) 356 END DO 357 CALL MPI_Barrier(MPI_COMM_WORLD,ierr) 358 CALL pauser(HALF) 359 IF (rank==0) WRITE(UNIT=*,FMT='(a)') "" 360 361 ! make sure increment is consistent everywhere, since 362 ! (1) system is nonlinear in the independent variables, and 363 ! (2) to ensure calculation of convergence criterion is correct 364 ! use 0-based MPI rank values 365 IF (rank==0) WRITE(UNIT=*,FMT='(a)') "calling MPI_Bcast now . . ." 366 DO i=1,gn 367 CALL MPI_Bcast(dT(i),1,MPI_DOUBLE_PRECISION,HomeProc(i),MPI_COMM_WORLD,ierr) 368 END DO 369 370 ! update solution vector 371 That=That+dT 372 373 ! plot solution 374 IF (rank==0) THEN 375 DO inp=1,nnp 376 plotdata(inp,2)=That(inp) 377 END DO 378 WRITE(UNIT=title,FMT='(a,i0)') "nonlinear iteration ",IterCount 379 PlotCount=PlotCount+1 380 WRITE(UNIT=PlotFilename,FMT='(a,i0)') "blah",PlotCount 381 CALL plotxy(plotdata=plotdata,filename=TRIM(PlotFilename), & 382 lines=lps,points=lps,widths=widths, & 383 title1=TRIM(title),output_animate=1) 384 END IF 385 386 ! update linear system 387 CALL UpdateLinearSystemThermal(nel,xhat,That,RHS,LHS,LHSassem) 388 389 ! check convergence 390 ResNorm=ZERO 391 IncNorm=ZERO 392 DO i=1,gn 393 ResNorm=ResNorm+RHS(i)**2 394 IncNorm=IncNorm+dT(i)**2 395 END DO 396 ResNorm=SQRT(ResNorm) 397 IncNorm=SQRT(IncNorm) 398 IF (IncNorm0==ZERO) IncNorm0=IncNorm 399 IF (rank==0) WRITE(UNIT=*,FMT='(a)') "" 400 IF (rank==0) WRITE(UNIT=*,FMT='(a,es22.15)') "absolute residual 2-norm=",ResNorm 401 IF (rank==0) WRITE(UNIT=*,FMT='(a,es22.15)') "relative residual 2-norm=",ResNorm/ResNorm0 402 IF (rank==0) WRITE(UNIT=*,FMT='(a,es22.15)') "absolute increment 2-norm=",IncNorm 403 IF (rank==0) WRITE(UNIT=*,FMT='(a,es22.15)') "relative increment 2-norm=",IncNorm/IncNorm0 404 IF (rank==0) WRITE(UNIT=*,FMT='(a,i0)') "number of LHS updates: ",NumberLHSUpdates 405 IF (ResNorm<MAX(AbsResTol,RelResTol*ResNorm0) .AND. IncNorm<MAX(AbsIncTol,RelIncTol*IncNorm0)) THEN 406 IF (rank==0) WRITE(UNIT=*,FMT='(a)') "" 407 IF (rank==0) WRITE(UNIT=*,FMT='(a)') "nonlinear termination criterion ""convergence"" satisfied, stopping" 408 EXIT 409 END IF 410 411 IF (IterCount==MaxIterCount) THEN 412 IF (rank==0) WRITE(UNIT=*,FMT='(a)') "" 413 IF (rank==0) WRITE(UNIT=*,FMT='(a,i0,a)') & 414 & "nonlinear termination criterion num_iter=", & 415 & MaxIterCount," satisfied, stopping" 416 EXIT 417 END IF 418 END DO 419 420 IF (rank==0) WRITE(UNIT=*,FMT='(a)') "" 421 422 CALL lis_precon_destroy(precon,ierr) 423 CALL lis_solver_destroy(solver,ierr) 424 CALL lis_vector_destroy(x,ierr) 425 CALL lis_vector_destroy(b,ierr) 426 CALL lis_matrix_destroy(A,ierr) 427 428 DEALLOCATE(dT) 429 DEALLOCATE(HomeProc) 430 431 IF (rank==0) THEN 432 DEALLOCATE(widths) 433 DEALLOCATE(lps) 434 DEALLOCATE(plotdata) 435 END IF 436 DEALLOCATE(LHSassem) 437 DEALLOCATE(LHS) 438 DEALLOCATE(RHS) 439 DEALLOCATE(That) 440 DEALLOCATE(xhat) 441 442 DEALLOCATE(proc_n) 443 444 CALL lis_finalize(ierr) 445! CALL ErrorCheck("lis_finalize",ierr) 446 447 CALL MPI_Finalize(ierr) 448! CALL ErrorCheck("MPI_Finalize",ierr) 449 450END PROGRAM lis_driver 451 452 453!========================================================================== 454SUBROUTINE pauser(pause_time) 455 456 REAL(kind(0.d0)),INTENT(IN) :: pause_time 457 458 INTEGER,PARAMETER :: singI=kind(0) 459 INTEGER,PARAMETER :: fullR=kind(0.d0) 460 INTEGER(singI) :: time_data(8) 461 REAL(fullR) :: start_time,end_time,dt 462 463 CALL DATE_AND_TIME(values=time_data) 464 start_time=time_data(5)*3600 + time_data(6)*60 + time_data(7) 465 DO 466 CALL DATE_AND_TIME(values=time_data) 467 end_time=time_data(5)*3600 + time_data(6)*60 + time_data(7) 468 dt=end_time-start_time 469 IF (dt>pause_time) EXIT 470 END DO 471 472END SUBROUTINE pauser 473 474 475!========================================================================== 476!SUBROUTINE ErrorCheck(FuncName,ierr) 477 478! CHARACTER(LEN=*) :: FuncName 479! INTEGER(kind(0)),INTENT(IN) :: ierr 480 481! IF (ierr/=0) THEN 482! WRITE(UNIT=*,FMT='(a,a,a,i0,a)') "ERROR: ",TRIM(FuncName)," returned error ",ierr,", quitting . . ." 483! STOP 484! END IF 485 486!END SUBROUTINE ErrorCheck 487 488 489!========================================================================== 490SUBROUTINE UpdateLinearSystemThermal(nel,xhat,That,RHS,LHS,LHSassem) 491 492 IMPLICIT NONE 493 494 INTEGER(kind(0)),INTENT(IN) :: nel 495 REAL(kind(0.d0)),INTENT(IN) :: xhat(nel+1),That(nel+1) 496 REAL(kind(0.d0)),INTENT(OUT) :: RHS(nel+1),LHS(nel+1,nel+1) 497 LOGICAL,INTENT(OUT) :: LHSassem(nel+1,nel+1) 498 499 INTEGER,PARAMETER :: singI=kind(0) 500 INTEGER,PARAMETER :: fullR=kind(0.d0) 501 502 REAL(fullR),PARAMETER :: ZERO = 0.0_fullR 503 REAL(fullR),PARAMETER :: ROOT_EPS = 0.0000000596046448_fullR ! 2^(-24) = 5.96046448E-8 504 REAL(fullR),PARAMETER :: ONE = 1.0_fullR 505 REAL(fullR),PARAMETER :: TWO = 2.0_fullR 506 REAL(fullR),PARAMETER :: THREE = 3.0_fullR 507 REAL(fullR),PARAMETER :: FOUR = 4.0_fullR 508 REAL(fullR),PARAMETER :: FIVE = 5.0_fullR 509 REAL(fullR),PARAMETER :: SIX = 6.0_fullR 510 REAL(fullR),PARAMETER :: SEVEN = 7.0_fullR 511 REAL(fullR),PARAMETER :: EIGHT = 8.0_fullR 512 REAL(fullR),PARAMETER :: NINE = 9.0_fullR 513 REAL(fullR),PARAMETER :: TEN = 10.0_fullR 514 REAL(fullR),PARAMETER :: HALF = 0.500000000000000_fullR 515 REAL(fullR),PARAMETER :: THIRD = 0.333333333333333_fullR 516 REAL(fullR),PARAMETER :: FOURTH = 0.250000000000000_fullR 517 REAL(fullR),PARAMETER :: TENTH = 0.100000000000000_fullR 518 REAL(fullR),PARAMETER :: PI = 3.141592653589793_fullR 519 520 LOGICAL :: BuildAssemFlag 521 REAL(fullR) :: xi(2),N(2,2) 522 INTEGER(singI) :: nnp,inp,jnp,iel,iip,DOF1,DOF2 523 REAL(fullR) :: xhat_elem(2),xbar_elem,That_elem(2) 524 REAL(fullR) :: dN_dxi(2,2),dN_dx(2) 525 REAL(fulLR) :: Tip,k,dk_dT,dT_dx 526 REAL(fullR) :: Relem(2),Kelem(2,2) 527 528 529 BuildAssemFlag=.FALSE. 530 IF (.NOT.LHSassem(1,1)) BuildAssemFlag=.TRUE. 531 532 ! initialization 533 nnp=nel+1 534 DO inp=1,nnp 535 RHS(inp)=ZERO 536 DO jnp=1,nnp 537 IF (BuildAssemFlag) LHSassem(inp,jnp)=.FALSE. 538 LHS(inp,jnp)=ZERO 539 END DO 540 END DO 541 542 xi(1)=-SQRT(THIRD) 543 xi(2)=+SQRT(THIRD) 544 545 CALL CalcN(xi,N) 546 CALL CalcdNdxi(dN_dxi) 547 548 DO iel=1,nel 549 ! local element data 550 DOF1=iel 551 DOF2=iel+1 552 553 xhat_elem(1)=xhat(DOF1) 554 xhat_elem(2)=xhat(DOF2) 555 xbar_elem=(xhat_elem(1)+xhat_elem(2))/TWO 556 That_elem(1)=That(DOF1) 557 That_elem(2)=That(DOF2) 558 Relem=ZERO 559 Kelem=ZERO 560 DO iip=1,2 561 dN_dx(1)=ONE/(dN_dxi(1,iip)*xhat_elem(1)+dN_dxi(2,iip)*xhat_elem(2))*dN_dxi(1,iip) 562 dN_dx(2)=ONE/(dN_dxi(1,iip)*xhat_elem(1)+dN_dxi(2,iip)*xhat_elem(2))*dN_dxi(2,iip) 563 ! Temperature at the current integration point 564 Tip=N(1,iip)*That_elem(1)+N(2,iip)*That_elem(2) 565 ! Calculate temperature-dependent material properties 566 CALL CalcCond(Tip,xbar_elem,k,dk_dT) 567 dT_dx=dN_dx(1)*That_elem(1)+dN_dx(2)*That_elem(2) 568 569 Relem=Relem+dN_dx*k*dT_dx 570 571 Kelem(1,1)=Kelem(1,1)+dN_dx(1)*k*dN_dx(1)+dN_dx(1)*dT_dx*dk_dT*N(1,iip) 572 Kelem(1,2)=Kelem(1,2)+dN_dx(1)*k*dN_dx(2)+dN_dx(1)*dT_dx*dk_dT*N(2,iip) 573 Kelem(2,1)=Kelem(2,1)+dN_dx(2)*k*dN_dx(1)+dN_dx(2)*dT_dx*dk_dT*N(1,iip) 574 Kelem(2,2)=Kelem(2,2)+dN_dx(2)*k*dN_dx(2)+dN_dx(2)*dT_dx*dk_dT*N(2,iip) 575 END DO 576 577 ! assemble 578 RHS(DOF1)=RHS(DOF1)+Relem(1) 579 RHS(DOF2)=RHS(DOF2)+Relem(2) 580 581 IF (BuildAssemFlag) THEN 582 LHSassem(DOF1,DOF1)=.TRUE. 583 LHSassem(DOF1,DOF2)=.TRUE. 584 LHSassem(DOF2,DOF1)=.TRUE. 585 LHSassem(DOF2,DOF2)=.TRUE. 586 END IF 587 LHS(DOF1,DOF1)=LHS(DOF1,DOF1)+Kelem(1,1) 588 LHS(DOF1,DOF2)=LHS(DOF1,DOF2)+Kelem(1,2) 589 LHS(DOF2,DOF1)=LHS(DOF2,DOF1)+Kelem(2,1) 590 LHS(DOF2,DOF2)=LHS(DOF2,DOF2)+Kelem(2,2) 591 END DO 592 593 ! apply dirichlet BC 594 RHS(1)=ZERO 595 RHS(nnp)=ZERO 596 DO inp=1,nnp 597 LHS(1,inp)=ZERO 598 LHS(inp,1)=ZERO 599 LHS(nnp,inp)=ZERO 600 LHS(inp,nnp)=ZERO 601 END DO 602 LHS(1,1)=ONE 603 LHS(nnp,nnp)=ONE 604 605END SUBROUTINE UpdateLinearSystemThermal 606 607 608!========================================================================== 609SUBROUTINE CalcN(xi,N) 610 611 IMPLICIT NONE 612 613 REAL(kind(0.d0)),INTENT(IN) :: xi(2) 614 REAL(kind(0.d0)),INTENT(OUT) :: N(2,2) 615 616 INTEGER,PARAMETER :: singI=kind(0) 617 INTEGER,PARAMETER :: fullR=kind(0.d0) 618 619 REAL(fullR),PARAMETER :: ZERO = 0.0_fullR 620 REAL(fullR),PARAMETER :: ROOT_EPS = 0.0000000596046448_fullR ! 2^(-24) = 5.96046448E-8 621 REAL(fullR),PARAMETER :: ONE = 1.0_fullR 622 REAL(fullR),PARAMETER :: TWO = 2.0_fullR 623 REAL(fullR),PARAMETER :: THREE = 3.0_fullR 624 REAL(fullR),PARAMETER :: FOUR = 4.0_fullR 625 REAL(fullR),PARAMETER :: FIVE = 5.0_fullR 626 REAL(fullR),PARAMETER :: SIX = 6.0_fullR 627 REAL(fullR),PARAMETER :: SEVEN = 7.0_fullR 628 REAL(fullR),PARAMETER :: EIGHT = 8.0_fullR 629 REAL(fullR),PARAMETER :: NINE = 9.0_fullR 630 REAL(fullR),PARAMETER :: TEN = 10.0_fullR 631 REAL(fullR),PARAMETER :: HALF = 0.500000000000000_fullR 632 REAL(fullR),PARAMETER :: THIRD = 0.333333333333333_fullR 633 REAL(fullR),PARAMETER :: FOURTH = 0.250000000000000_fullR 634 REAL(fullR),PARAMETER :: TENTH = 0.100000000000000_fullR 635 REAL(fullR),PARAMETER :: PI = 3.141592653589793_fullR 636 637 ! ip=1 638 N(1,1)=HALF*(ONE-xi(1)) 639 N(2,1)=HALF*(ONE+xi(1)) 640 641 ! ip=2 642 N(1,2)=HALF*(ONE-xi(2)) 643 N(2,2)=HALF*(ONE+xi(2)) 644 645END SUBROUTINE CalcN 646 647 648!========================================================================== 649! only need the integration points themselves here for higher order elements 650!========================================================================== 651SUBROUTINE CalcdNdxi(dN_dxi) 652 653 IMPLICIT NONE 654 655 REAL(kind(0.d0)),INTENT(OUT) :: dN_dxi(2,2) 656 657 INTEGER,PARAMETER :: singI=kind(0) 658 INTEGER,PARAMETER :: fullR=kind(0.d0) 659 660 REAL(fullR),PARAMETER :: ZERO = 0.0_fullR 661 REAL(fullR),PARAMETER :: ROOT_EPS = 0.0000000596046448_fullR ! 2^(-24) = 5.96046448E-8 662 REAL(fullR),PARAMETER :: ONE = 1.0_fullR 663 REAL(fullR),PARAMETER :: TWO = 2.0_fullR 664 REAL(fullR),PARAMETER :: THREE = 3.0_fullR 665 REAL(fullR),PARAMETER :: FOUR = 4.0_fullR 666 REAL(fullR),PARAMETER :: FIVE = 5.0_fullR 667 REAL(fullR),PARAMETER :: SIX = 6.0_fullR 668 REAL(fullR),PARAMETER :: SEVEN = 7.0_fullR 669 REAL(fullR),PARAMETER :: EIGHT = 8.0_fullR 670 REAL(fullR),PARAMETER :: NINE = 9.0_fullR 671 REAL(fullR),PARAMETER :: TEN = 10.0_fullR 672 REAL(fullR),PARAMETER :: HALF = 0.500000000000000_fullR 673 REAL(fullR),PARAMETER :: THIRD = 0.333333333333333_fullR 674 REAL(fullR),PARAMETER :: FOURTH = 0.250000000000000_fullR 675 REAL(fullR),PARAMETER :: TENTH = 0.100000000000000_fullR 676 REAL(fullR),PARAMETER :: PI = 3.141592653589793_fullR 677 678 ! ip=1 679 dN_dxi(1,1)=-HALF 680 dN_dxi(2,1)=+HALF 681 682 ! ip=2 683 dN_dxi(1,2)=-HALF 684 dN_dxi(2,2)=+HALF 685 686END SUBROUTINE CalcdNdxi 687 688 689!========================================================================== 690SUBROUTINE CalcCond(T,xbar,k,dk_dT) 691 692 IMPLICIT NONE 693 694 REAL(kind(0.d0)),INTENT(IN) :: T,xbar 695 REAL(kind(0.d0)),INTENT(OUT) :: k,dk_dT 696 697 INTEGER,PARAMETER :: singI=kind(0) 698 INTEGER,PARAMETER :: fullR=kind(0.d0) 699 700 REAL(fullR),PARAMETER :: ZERO = 0.0_fullR 701 REAL(fullR),PARAMETER :: ROOT_EPS = 0.0000000596046448_fullR ! 2^(-24) = 5.96046448E-8 702 REAL(fullR),PARAMETER :: ONE = 1.0_fullR 703 REAL(fullR),PARAMETER :: TWO = 2.0_fullR 704 REAL(fullR),PARAMETER :: THREE = 3.0_fullR 705 REAL(fullR),PARAMETER :: FOUR = 4.0_fullR 706 REAL(fullR),PARAMETER :: FIVE = 5.0_fullR 707 REAL(fullR),PARAMETER :: SIX = 6.0_fullR 708 REAL(fullR),PARAMETER :: SEVEN = 7.0_fullR 709 REAL(fullR),PARAMETER :: EIGHT = 8.0_fullR 710 REAL(fullR),PARAMETER :: NINE = 9.0_fullR 711 REAL(fullR),PARAMETER :: TEN = 10.0_fullR 712 REAL(fullR),PARAMETER :: HALF = 0.500000000000000_fullR 713 REAL(fullR),PARAMETER :: THIRD = 0.333333333333333_fullR 714 REAL(fullR),PARAMETER :: FOURTH = 0.250000000000000_fullR 715 REAL(fullR),PARAMETER :: TENTH = 0.100000000000000_fullR 716 REAL(fullR),PARAMETER :: PI = 3.141592653589793_fullR 717 718 REAL(fullR) :: TabT(16),T1,T2 719 REAL(fullR) :: Tabk(16),k1,k2 720 INTEGER(singI) :: NumInts,iInt 721 LOGICAL :: FoundFlag 722 723 TabT(1)=-1.0e+10 724 TabT(2)=273.0 725 TabT(3)=431.6 726 TabT(4)=590.0 727 TabT(5)=748.7 728 TabT(6)=907.2 729 TabT(7)=1065.8 730 TabT(8)=1224.3 731 TabT(9)=1382.9 732 TabT(10)=1541.4 733 TabT(11)=1650.0 734 TabT(12)=1660.0 735 TabT(13)=1690.0 736 TabT(14)=1700.0 737 TabT(15)=3200.0 738 TabT(16)=1.0e+10 739 740 Tabk(1)=3.0e-02 741 Tabk(2)=3.0e-02 742 Tabk(3)=3.0e-02 743 Tabk(4)=3.0e-02 744 Tabk(5)=3.0e-02 745 Tabk(6)=3.0e-02 746 Tabk(7)=3.0e-02 747 Tabk(8)=3.0e-02 748 Tabk(9)=3.0e-02 749 Tabk(10)=3.0e-02 750 Tabk(11)=3.0e-02 751 Tabk(12)=3.0e-02 752 Tabk(13)=3.0e-02 753 Tabk(14)=3.0e-02 754 Tabk(15)=3.0e-02 755 Tabk(16)=3.0e-02 756 !============================================= 757 ! scaling to induce ill-conditioning . . . 758 !============================================= 759 IF (xbar<=HALF) THEN 760! Tabk=Tabk*(TENTH**2) 761 ELSE IF (HALF<xbar) THEN 762! Tabk=Tabk*(TEN**2) 763 END IF 764 !============================================= 765 ! scaling to induce nonlinearity . . . 766 !============================================= 767 Tabk(1)=Tabk(1)*ONE 768 Tabk(2)=Tabk(2)*TWO 769 Tabk(3)=Tabk(3)*FIVE 770 Tabk(4)=Tabk(4)*TEN 771 Tabk(5)=Tabk(5)*TEN*TWO 772 Tabk(6)=Tabk(6)*TEN*FIVE 773 Tabk(7)=Tabk(7)*TEN*EIGHT 774 Tabk(8)=Tabk(8)*TEN*TEN 775 Tabk(9)=Tabk(9)*TEN*NINE 776 Tabk(10)=Tabk(10)*TEN*EIGHT 777 Tabk(11)=Tabk(11)*TEN*SEVEN 778 Tabk(12)=Tabk(12)*TEN*SIX 779 Tabk(13)=Tabk(13)*TEN*FIVE 780 Tabk(14)=Tabk(14)*TEN*FOUR 781 Tabk(15)=Tabk(15)*TEN*THREE 782 Tabk(16)=Tabk(16)*TEN*TWO 783 784 785 NumInts=15 786 FoundFlag=.FALSE. 787 DO iInt=1,NumInts 788 T1=TabT(iInt) 789 T2=TabT(iInt+1) 790 IF (T<T1) CYCLE 791 IF (T>T2) CYCLE 792 793 FoundFlag=.TRUE. 794 k1=Tabk(iInt) 795 k2=Tabk(iInt+1) 796 k=k1 + (T-T1)/(T2-T1)*(k2-k1) 797 dk_dT=(k2-k1)/(T2-T1) 798 EXIT 799 END DO 800 801 IF (.NOT.FoundFlag) THEN 802 WRITE(UNIT=*,FMT='(a)') "ERROR: conductivities not calculated, quitting . . ." 803 STOP 804 END IF 805 806END SUBROUTINE CalcCond 807 808 809!========================================================================== 810! Plot (x,y) data using gnuplot 811! * legend removal 812! * series names 813! * axes names 814! * title 815! 816! NOTE: There is a "bug" (probably) in gnuplot. Basically, "term wxt", which is 817! required if one wants to do fancy zooming in and out of plots, does not work 818! for the way I originally had my plot files structured. Thus, I had to re-work 819! them, so that the data is in a seperate file. Lame . . . 820! 821! Some of the behavior of this subroutine, as well as the options shown below 822! require some explanation . . . 823! 824! * the data is assumed to be in one of two formats, the first being 825! 826! x y1 y2 . . . yn 827! ======================= 828! # # # # 829! # # # # 830! # # # # 831! . . . . . . . . . . . 832! 833! For this type of data, it is assumed that each data series y1, y2, . 834! . ., yn has a value corresponding to each value of x. 835! 836! Alternately, for data series that do not share a common domain, the data 837! should be formatted like 838! 839! x1 y1 840! ======== 841! # # 842! # # 843! # # 844! x2 y2 845! ======== 846! # # 847! # # 848! # # 849! x3 y3 850! ======== 851! # # 852! # # 853! # # 854! . 855! . 856! . 857! xn yn 858! ======== 859! # # 860! # # 861! # # 862! 863! In this case, the variable "datalength" should also be passed in, 864! which defines the length of each data series. 865! 866! * color codes: 867! 1=black 868! 2=blue 869! 3=red 870! 4=green 871! 872!========================================================================== 873SUBROUTINE plotxy(plotdata,filename,datalength, & 874 grid,limits,title1,xlabel,ylabel,series_labels, & 875 legend,lines,points,colors,widths,output_animate) 876 877 REAL(kind(0.d0)),INTENT(IN) :: plotdata(:,:) 878 CHARACTER(*),INTENT(IN) :: filename 879 INTEGER,INTENT(IN),OPTIONAL :: datalength(:) 880 INTEGER,INTENT(IN),OPTIONAL :: grid 881 REAL(kind(0.d0)),INTENT(IN),OPTIONAL :: limits(:,:) 882 CHARACTER(LEN=*),INTENT(IN),OPTIONAL :: title1 883 CHARACTER(LEN=*),INTENT(IN),OPTIONAL :: xlabel,ylabel 884 CHARACTER(LEN=*),INTENT(IN),OPTIONAL :: series_labels(:) 885 INTEGER,INTENT(IN),OPTIONAL :: legend 886 INTEGER,INTENT(IN),OPTIONAL :: lines(:) 887 INTEGER,INTENT(IN),OPTIONAL :: points(:) 888 INTEGER,INTENT(IN),OPTIONAL :: colors(:) 889 INTEGER,INTENT(IN),OPTIONAL :: widths(:) 890 INTEGER,INTENT(IN),OPTIONAL :: output_animate 891 892 INTEGER,PARAMETER :: singI=kind(0) 893 INTEGER,PARAMETER :: fullR=kind(0.d0) 894 INTEGER(singI),PARAMETER :: PLOT_FID=3 895 896 REAL(fullR),PARAMETER :: ZERO = 0.0_fullR 897 REAL(fullR),PARAMETER :: ROOT_EPS = 0.0000000596046448_fullR ! 2^(-24) = 5.96046448E-8 898 REAL(fullR),PARAMETER :: ONE = 1.0_fullR 899 REAL(fullR),PARAMETER :: TWO = 2.0_fullR 900 REAL(fullR),PARAMETER :: THREE = 3.0_fullR 901 REAL(fullR),PARAMETER :: FOUR = 4.0_fullR 902 REAL(fullR),PARAMETER :: FIVE = 5.0_fullR 903 REAL(fullR),PARAMETER :: SIX = 6.0_fullR 904 REAL(fullR),PARAMETER :: SEVEN = 7.0_fullR 905 REAL(fullR),PARAMETER :: EIGHT = 8.0_fullR 906 REAL(fullR),PARAMETER :: NINE = 9.0_fullR 907 REAL(fullR),PARAMETER :: TEN = 10.0_fullR 908 REAL(fullR),PARAMETER :: HALF = 0.500000000000000_fullR 909 REAL(fullR),PARAMETER :: THIRD = 0.333333333333333_fullR 910 REAL(fullR),PARAMETER :: FOURTH = 0.250000000000000_fullR 911 REAL(fullR),PARAMETER :: TENTH = 0.100000000000000_fullR 912 REAL(fullR),PARAMETER :: PI = 3.141592653589793_fullR 913 914 INTEGER(singI) :: length,nseries,i 915 REAL(fullR) :: minx,maxx 916 CHARACTER(LEN=1024) :: textline 917 INTEGER(singI) :: llegend,loutput_animate 918 INTEGER(singI),POINTER :: llines(:),lpoints(:),lcolors(:),lwidths(:) 919 INTEGER(singI) :: linestat,pointstat,colorstat,widthstat 920 INTEGER(singI) :: iseries,ndatalength,idatalength,startrow,endrow 921 CHARACTER(LEN=1024),POINTER :: filename_complete(:) 922 923 INTEGER :: debug=0 924 925 IF (debug>=1) THEN 926 WRITE(UNIT=*,FMT='(a)') "status of optional parameters:" 927 IF (PRESENT(datalength)) THEN 928 WRITE(UNIT=*,FMT='(a)') " ""datalength"" exists" 929 ELSE 930 WRITE(UNIT=*,FMT='(a)') " ""datalength"" does not exist" 931 END IF 932 IF (PRESENT(grid)) THEN 933 WRITE(UNIT=*,FMT='(a)') " ""grid"" exists" 934 ELSE 935 WRITE(UNIT=*,FMT='(a)') " ""grid"" does not exist" 936 END IF 937 IF (PRESENT(limits)) THEN 938 WRITE(UNIT=*,FMT='(a)') " ""limits"" exists" 939 ELSE 940 WRITE(UNIT=*,FMT='(a)') " ""limits"" does not exist" 941 END IF 942 IF (PRESENT(title1)) THEN 943 WRITE(UNIT=*,FMT='(a)') " ""title"" exists" 944 ELSE 945 WRITE(UNIT=*,FMT='(a)') " ""title"" does not exist" 946 END IF 947 IF (PRESENT(xlabel)) THEN 948 WRITE(UNIT=*,FMT='(a)') " ""xlabel"" exists" 949 ELSE 950 WRITE(UNIT=*,FMT='(a)') " ""xlabel"" does not exist" 951 END IF 952 IF (PRESENT(ylabel)) THEN 953 WRITE(UNIT=*,FMT='(a)') " ""ylabel"" exists" 954 ELSE 955 WRITE(UNIT=*,FMT='(a)') " ""ylabel"" does not exist" 956 END IF 957 IF (PRESENT(series_labels)) THEN 958 WRITE(UNIT=*,FMT='(a)') " ""series_labels"" exists" 959 ELSE 960 WRITE(UNIT=*,FMT='(a)') " ""series_labels"" does not exist" 961 END IF 962 IF (PRESENT(legend)) THEN 963 WRITE(UNIT=*,FMT='(a)') " ""legend"" exists" 964 ELSE 965 WRITE(UNIT=*,FMT='(a)') " ""legend"" does not exist" 966 END IF 967 IF (PRESENT(lines)) THEN 968 WRITE(UNIT=*,FMT='(a)') " ""lines"" exists" 969 ELSE 970 WRITE(UNIT=*,FMT='(a)') " ""lines"" does not exist" 971 END IF 972 IF (PRESENT(points)) THEN 973 WRITE(UNIT=*,FMT='(a)') " ""points"" exists" 974 ELSE 975 WRITE(UNIT=*,FMT='(a)') " ""points"" does not exist" 976 END IF 977 IF (PRESENT(colors)) THEN 978 WRITE(UNIT=*,FMT='(a)') " ""colors"" exists" 979 ELSE 980 WRITE(UNIT=*,FMT='(a)') " ""colors"" does not exist" 981 END IF 982 IF (PRESENT(widths)) THEN 983 WRITE(UNIT=*,FMT='(a)') " ""widths"" exists" 984 ELSE 985 WRITE(UNIT=*,FMT='(a)') " ""widths"" does not exist" 986 END IF 987 IF (PRESENT(output_animate)) THEN 988 WRITE(UNIT=*,FMT='(a)') " ""output_animate"" exists" 989 ELSE 990 WRITE(UNIT=*,FMT='(a)') " ""output_animate"" does not exist" 991 END IF 992 END IF 993 994 length=SIZE(plotdata,1) 995 nseries=SIZE(plotdata,2)-1 996 IF (PRESENT(datalength)) THEN 997 ndatalength=SIZE(datalength) 998 ELSE 999 ndatalength=1 1000 END IF 1001 1002 OPEN(UNIT=PLOT_FID,FILE=TRIM(filename)//".plot",ACTION="WRITE") 1003 IF (debug>=1) THEN 1004 WRITE(UNIT=*,FMT='(a)') "writing plot commands" 1005 END IF 1006 1007 loutput_animate=0 1008 IF (PRESENT(output_animate)) THEN 1009 loutput_animate=output_animate 1010 END IF 1011 1012 ! 1013 ! write to a file, as opposed to the screen 1014 ! 1015 ! note that it is not guaranteed that all plot types will be available for 1016 ! all gnuplot installations 1017 ! 1018 IF (loutput_animate==1) THEN 1019 WRITE(UNIT=PLOT_FID,FMT='(a,a,a)') "set terminal png truecolor size 1600,1200; set output '",TRIM(filename),".png'" 1020 ELSE 1021 WRITE(UNIT=PLOT_FID,FMT='(a)') "set terminal wxt size 1600,1200" 1022! WRITE(UNIT=PLOT_FID,FMT='(a)') "set terminal X11" 1023 END IF 1024 1025 ! 1026 ! set min and max value for x and y axes 1027 ! 1028 IF (.not.PRESENT(limits)) THEN 1029 minx=MINVAL(plotdata(:,1)) 1030 maxx=MAXVAL(plotdata(:,1)) 1031 WRITE(UNIT=PLOT_FID,FMT='(a,es22.15,a,es22.15,a)') "set xrange [",(minx-(maxx-minx)/TEN),":",(maxx+(maxx-minx)/TEN),"]" 1032 minx = HUGE(minx) 1033 maxx = -HUGE(maxx) 1034 DO i = 1,nseries 1035 minx=MIN(minx,MINVAL(plotdata(:,i+1))) 1036 maxx=MAX(maxx,MAXVAL(plotdata(:,i+1))) 1037 END DO 1038 WRITE(UNIT=PLOT_FID,FMT='(a,es22.15,a,es22.15,a)') "set yrange [",(minx-(maxx-minx)/TEN),":",(maxx+(maxx-minx)/TEN),"]" 1039 ELSE 1040 minx=limits(1,1) 1041 maxx=limits(1,2) 1042 WRITE(UNIT=PLOT_FID,FMT='(a,es22.15,a,es22.15,a)') "set xrange [",(minx),":",(maxx),"]" 1043 minx=limits(2,1) 1044 maxx=limits(2,2) 1045 WRITE(UNIT=PLOT_FID,FMT='(a,es22.15,a,es22.15,a)') "set yrange [",(minx),":",(maxx),"]" 1046 END IF 1047 1048 ! 1049 ! set the aspect ratio 1050 ! 1051 WRITE(UNIT=PLOT_FID,FMT='(a)') "set size square" 1052 1053 ! 1054 ! turn on the grid 1055 ! 1056 IF (PRESENT(grid)) THEN 1057 IF (grid==1) WRITE(UNIT=PLOT_FID,FMT='(a)') "set grid" 1058 ELSE 1059 WRITE(UNIT=PLOT_FID,FMT='(a)') "set grid" 1060 END IF 1061 1062 ! 1063 ! add a title 1064 ! 1065 IF (PRESENT(title1)) THEN 1066 WRITE(UNIT=PLOT_FID,FMT='(a,a,a)') "set title """,TRIM(title1),"""" 1067 END IF 1068 1069 ! 1070 ! add axis labels 1071 ! 1072 IF (PRESENT(xlabel)) WRITE(UNIT=PLOT_FID,FMT='(a,a,a)') "set xlabel """,TRIM(xlabel),"""" 1073 IF (PRESENT(ylabel)) WRITE(UNIT=PLOT_FID,FMT='(a,a,a)') "set ylabel """,TRIM(ylabel),"""" 1074 1075 ! 1076 ! add various other stuff 1077 ! 1078 llegend=0 1079 IF (PRESENT(legend)) THEN 1080 llegend=legend 1081 END IF 1082 ALLOCATE(llines(1)) 1083 llines=0 1084 IF (PRESENT(lines)) THEN 1085 DEALLOCATE(llines) 1086 ALLOCATE(llines(SIZE(lines,1))) 1087 llines=lines 1088 END IF 1089 ALLOCATE(lpoints(1)) 1090 lpoints=0 1091 IF (PRESENT(points)) THEN 1092 DEALLOCATE(lpoints) 1093 ALLOCATE(lpoints(SIZE(points,1))) 1094 lpoints=points 1095 END IF 1096 ALLOCATE(lcolors(1)) 1097 lcolors=1 1098 IF (PRESENT(colors)) THEN 1099 DEALLOCATE(lcolors) 1100 ALLOCATE(lcolors(SIZE(colors,1))) 1101 lcolors=colors 1102 END IF 1103 ALLOCATE(lwidths(1)) 1104 lwidths=1 1105 IF (PRESENT(widths)) THEN 1106 DEALLOCATE(lwidths) 1107 ALLOCATE(lwidths(SIZE(widths,1))) 1108 lwidths=widths 1109 END IF 1110 1111 IF (llegend==0) THEN 1112 WRITE(UNIT=PLOT_FID,FMT='(a,a,a)') "set nokey" 1113 END IF 1114 1115 ! create the names of the data files and store them 1116 IF (nseries>1) THEN 1117 ALLOCATE(filename_complete(nseries)) 1118 ELSE IF (ndatalength>=1) THEN 1119 ALLOCATE(filename_complete(ndatalength)) 1120 END IF 1121 IF (nseries>1) THEN 1122 DO iseries=1,nseries 1123 WRITE(UNIT=filename_complete(iseries),FMT='(a,a,i0.1,a)') TRIM(filename),"_data",iseries,".plot" 1124 IF (debug>=1) WRITE(UNIT=*,FMT='(a,a,a)') "creating plot data file named """,TRIM(filename_complete(iseries)),"""" 1125 END DO 1126 ELSE IF (nseries==1) THEN 1127 DO idatalength=1,ndatalength 1128 WRITE(UNIT=filename_complete(idatalength),FMT='(a,a,i0.1,a)') TRIM(filename),"_data",idatalength,".plot" 1129 IF (debug>=1) WRITE(UNIT=*,FMT='(a,a,a)') "creating plot data file named """,TRIM(filename_complete(idatalength)),"""" 1130 END DO 1131 END IF 1132 1133 IF (nseries>1) THEN 1134 DO iseries=1,nseries 1135 IF (iseries==1) THEN 1136 WRITE(UNIT=textline,FMT='(a,a,a)') "plot """,TRIM(filename_complete(iseries)),"""" 1137 ELSE 1138 WRITE(UNIT=textline,FMT='(a,a,a)') " """,TRIM(filename_complete(iseries)),"""" 1139 END IF 1140 IF (PRESENT(series_labels)) THEN 1141 WRITE(UNIT=textline,FMT='(a,a,a,a)') TRIM(textline), " title """,TRIM(series_labels(iseries)),"""" 1142 END IF 1143 IF (SIZE(llines,1)==1) THEN 1144 linestat=llines(1) 1145 ELSE 1146 linestat=llines(iseries) 1147 END IF 1148 IF (SIZE(lpoints,1)==1) THEN 1149 pointstat=lpoints(1) 1150 ELSE 1151 pointstat=lpoints(iseries) 1152 END IF 1153 IF ((linestat==1).AND.(pointstat==0)) THEN 1154 WRITE(UNIT=textline,FMT='(a,a)') TRIM(textline), " with lines" 1155 ELSE IF ((linestat==0).AND.(pointstat==1)) THEN 1156 WRITE(UNIT=textline,FMT='(a,a)') TRIM(textline), " with points" 1157 ELSE IF ((linestat==1).AND.(pointstat==1)) THEN 1158 WRITE(UNIT=textline,FMT='(a,a)') TRIM(textline), " with linespoints" 1159 END IF 1160 colorstat=0 1161 IF (SIZE(lcolors,1)==1) THEN 1162 colorstat=lcolors(1) 1163 ELSE 1164 colorstat=lcolors(iseries) 1165 END IF 1166 IF (colorstat==1) THEN 1167 WRITE(UNIT=textline,FMT='(a,a)') TRIM(textline), " linecolor -1" 1168 ELSE IF (colorstat==2) THEN 1169 WRITE(UNIT=textline,FMT='(a,a)') TRIM(textline), " linecolor 3" 1170 ELSE IF (colorstat==3) THEN 1171 WRITE(UNIT=textline,FMT='(a,a)') TRIM(textline), " linecolor 1" 1172 ELSE IF (colorstat==4) THEN 1173 WRITE(UNIT=textline,FMT='(a,a)') TRIM(textline), " linecolor 2" 1174 END IF 1175 widthstat=0 1176 IF (SIZE(lwidths,1)==1) THEN 1177 widthstat=lwidths(1) 1178 ELSE 1179 widthstat=lwidths(iseries) 1180 END IF 1181 IF (widthstat>1) THEN 1182 WRITE(UNIT=textline,FMT='(a,a,i0.1)') TRIM(textline), " linewidth ",widthstat 1183 END IF 1184 IF ((nseries>1) .AND. (iseries<nseries)) THEN 1185 WRITE(UNIT=textline,FMT='(a,a)') TRIM(textline), ", \" 1186 END IF 1187 WRITE(UNIT=PLOT_FID,FMT='(a)') TRIM(textline) 1188 END DO 1189 ELSE IF (nseries==1) THEN 1190 DO idatalength=1,ndatalength 1191 IF (idatalength==1) THEN 1192 WRITE(UNIT=textline,FMT='(a,a,a)') "plot """,TRIM(filename_complete(idatalength)),"""" 1193 ELSE 1194 WRITE(UNIT=textline,FMT='(a,a,a)') " """,TRIM(filename_complete(idatalength)),"""" 1195 END IF 1196 IF (PRESENT(series_labels)) THEN 1197 WRITE(UNIT=textline,FMT='(a,a,a,a)') TRIM(textline), " title """,TRIM(series_labels(idatalength)),"""" 1198 END IF 1199 IF (SIZE(llines,1)==1) THEN 1200 linestat=llines(1) 1201 ELSE 1202 linestat=llines(idatalength) 1203 END IF 1204 IF (SIZE(lpoints,1)==1) THEN 1205 pointstat=lpoints(1) 1206 ELSE 1207 pointstat=lpoints(idatalength) 1208 END IF 1209 IF ((linestat==1).AND.(pointstat==0)) THEN 1210 WRITE(UNIT=textline,FMT='(a,a)') TRIM(textline), " with lines" 1211 ELSE IF ((linestat==0).AND.(pointstat==1)) THEN 1212 WRITE(UNIT=textline,FMT='(a,a)') TRIM(textline), " with points" 1213 ELSE IF ((linestat==1).AND.(pointstat==1)) THEN 1214 WRITE(UNIT=textline,FMT='(a,a)') TRIM(textline), " with linespoints" 1215 END IF 1216 colorstat=0 1217 IF (SIZE(lcolors,1)==1) THEN 1218 colorstat=lcolors(1) 1219 ELSE 1220 colorstat=lcolors(idatalength) 1221 END IF 1222 IF (colorstat==1) THEN 1223 WRITE(UNIT=textline,FMT='(a,a)') TRIM(textline), " linecolor -1" 1224 ELSE IF (colorstat==2) THEN 1225 WRITE(UNIT=textline,FMT='(a,a)') TRIM(textline), " linecolor 3" 1226 ELSE IF (colorstat==3) THEN 1227 WRITE(UNIT=textline,FMT='(a,a)') TRIM(textline), " linecolor 1" 1228 ELSE IF (colorstat==4) THEN 1229 WRITE(UNIT=textline,FMT='(a,a)') TRIM(textline), " linecolor 2" 1230 END IF 1231 widthstat=0 1232 IF (SIZE(lwidths,1)==1) THEN 1233 widthstat=lwidths(1) 1234 ELSE 1235 widthstat=lwidths(idatalength) 1236 END IF 1237 IF (widthstat>1) THEN 1238 WRITE(UNIT=textline,FMT='(a,a,i0.1)') TRIM(textline), " linewidth ",widthstat 1239 END IF 1240 IF ((ndatalength>1) .AND. (idatalength<ndatalength)) THEN 1241 WRITE(UNIT=textline,FMT='(a,a)') TRIM(textline), ", \" 1242 END IF 1243 WRITE(UNIT=PLOT_FID,FMT='(a)') TRIM(textline) 1244 END DO 1245 END IF 1246 1247 1248 IF (loutput_animate==0) THEN 1249 WRITE(UNIT=PLOT_FID,FMT='(a)') "pause -1" 1250 1251 ! do it all over again, to create figures for journal articles 1252 WRITE(UNIT=PLOT_FID,FMT='(a)') "" 1253 WRITE(UNIT=PLOT_FID,FMT='(a)') "# comment the following line to create a journal-formatted .eps file" 1254 WRITE(UNIT=PLOT_FID,FMT='(a)') "exit" 1255 WRITE(UNIT=PLOT_FID,FMT='(a)') "" 1256 1257 WRITE(UNIT=PLOT_FID,FMT='(a,a,a)') & 1258 "set terminal postscript enhanced color; set output '",TRIM(filename),".eps'" 1259 WRITE(UNIT=PLOT_FID,FMT='(a)') & 1260 "# use the following command to convert the .eps file to a .pdf file:" 1261 WRITE(UNIT=PLOT_FID,FMT='(a,a,a,a,a)') & 1262 "# ps2pdf -dEPSCrop ",TRIM(filename),".eps ",TRIM(filename),".pdf" 1263 1264 ! 1265 ! set min and max value for x and y axes 1266 ! 1267 IF (.NOT.PRESENT(limits)) THEN 1268 minx=MINVAL(plotdata(:,1)) 1269 maxx=MAXVAL(plotdata(:,1)) 1270 WRITE(UNIT=PLOT_FID,FMT='(a,es22.15,a,es22.15,a)') "set xrange [",(minx-(maxx-minx)/TEN),":",(maxx+(maxx-minx)/TEN),"]" 1271 minx = HUGE(minx) 1272 maxx = -HUGE(maxx) 1273 DO i = 1,nseries 1274 minx=MIN(minx,MINVAL(plotdata(:,i+1))) 1275 maxx=MAX(maxx,MAXVAL(plotdata(:,i+1))) 1276 END DO 1277 WRITE(UNIT=PLOT_FID,FMT='(a,es22.15,a,es22.15,a)') "set yrange [",(minx-(maxx-minx)/TEN),":",(maxx+(maxx-minx)/TEN),"]" 1278 ELSE 1279 minx=limits(1,1) 1280 maxx=limits(1,2) 1281 WRITE(UNIT=PLOT_FID,FMT='(a,es22.15,a,es22.15,a)') "set xrange [",(minx),":",(maxx),"]" 1282 minx=limits(2,1) 1283 maxx=limits(2,2) 1284 WRITE(UNIT=PLOT_FID,FMT='(a,es22.15,a,es22.15,a)') "set yrange [",(minx),":",(maxx),"]" 1285 END IF 1286 1287 ! 1288 ! set the aspect ratio 1289 ! 1290 WRITE(UNIT=PLOT_FID,FMT='(a)') "set size square" 1291 1292 ! 1293 ! turn on the grid 1294 ! 1295 IF (PRESENT(grid)) THEN 1296 IF (grid==1) THEN 1297 WRITE(UNIT=PLOT_FID,FMT='(a)') "set grid" 1298 END IF 1299 ELSE 1300 WRITE(UNIT=PLOT_FID,FMT='(a)') "set grid" 1301 END IF 1302 1303 ! 1304 ! add a title 1305 ! 1306 IF (PRESENT(title1)) THEN 1307 WRITE(UNIT=PLOT_FID,FMT='(a)') "set title" 1308 END IF 1309 1310 ! 1311 ! add axis labels 1312 ! 1313 WRITE(UNIT=PLOT_FID,FMT='(a)') "set xlabel" 1314 WRITE(UNIT=PLOT_FID,FMT='(a)') "set ylabel" 1315 1316 ! 1317 ! add various other stuff 1318 ! 1319 llegend=0 1320 IF (PRESENT(legend)) THEN 1321 llegend=legend 1322 END IF 1323 ALLOCATE(llines(1)) 1324 llines=0 1325 IF (PRESENT(lines)) THEN 1326 DEALLOCATE(llines) 1327 ALLOCATE(llines(SIZE(lines,1))) 1328 llines=lines 1329 END IF 1330 ALLOCATE(lpoints(1)) 1331 lpoints=0 1332 IF (PRESENT(points)) THEN 1333 DEALLOCATE(lpoints) 1334 ALLOCATE(lpoints(SIZE(points,1))) 1335 lpoints=points 1336 END IF 1337 ALLOCATE(lcolors(1)) 1338 lcolors=1 1339 IF (PRESENT(colors)) THEN 1340 DEALLOCATE(lcolors) 1341 ALLOCATE(lcolors(SIZE(colors,1))) 1342 lcolors=colors 1343 END IF 1344 ALLOCATE(lwidths(1)) 1345 lwidths=1 1346 IF (PRESENT(widths)) THEN 1347 DEALLOCATE(lwidths) 1348 ALLOCATE(lwidths(SIZE(widths,1))) 1349 lwidths=widths 1350 END IF 1351 1352 IF (llegend==0) THEN 1353 WRITE(UNIT=PLOT_FID,FMT='(a,a,a)') "set nokey" 1354 END IF 1355 1356 ! create the names of the data files and store them 1357 IF (nseries>1) THEN 1358 ALLOCATE(filename_complete(nseries)) 1359 ELSE IF (ndatalength>=1) THEN 1360 ALLOCATE(filename_complete(ndatalength)) 1361 END IF 1362 IF (nseries>1) THEN 1363 DO iseries=1,nseries 1364 WRITE(UNIT=filename_complete(iseries),FMT='(a,a,i0.1,a)') TRIM(filename),"_data",iseries,".plot" 1365 IF (debug>=1) THEN 1366 WRITE(UNIT=*,FMT='(a,a,a)') "creating plot data file named """,TRIM(filename_complete(iseries)),"""" 1367 END IF 1368 END DO 1369 ELSE IF (nseries==1) THEN 1370 DO idatalength=1,ndatalength 1371 WRITE(UNIT=filename_complete(idatalength),FMT='(a,a,i0.1,a)') TRIM(filename),"_data",idatalength,".plot" 1372 IF (debug>=1) THEN 1373 WRITE(UNIT=*,FMT='(a,a,a)') "creating plot data file named """,TRIM(filename_complete(idatalength)),"""" 1374 END IF 1375 END DO 1376 END IF 1377 1378 IF (nseries>1) THEN 1379 DO iseries=1,nseries 1380 IF (iseries==1) THEN 1381 WRITE(UNIT=textline,FMT='(a,a,a)') "plot """,TRIM(filename_complete(iseries)),"""" 1382 ELSE 1383 WRITE(UNIT=textline,FMT='(a,a,a)') " """,TRIM(filename_complete(iseries)),"""" 1384 END IF 1385 IF (PRESENT(series_labels)) THEN 1386 WRITE(UNIT=textline,FMT='(a,a,a,a)') TRIM(textline), " title """,TRIM(series_labels(iseries)),"""" 1387 END IF 1388 IF (SIZE(llines,1)==1) THEN 1389 linestat=llines(1) 1390 ELSE 1391 linestat=llines(iseries) 1392 END IF 1393 IF (SIZE(lpoints,1)==1) THEN 1394 pointstat=lpoints(1) 1395 ELSE 1396 pointstat=lpoints(iseries) 1397 END IF 1398 WRITE(UNIT=textline,FMT='(a,a)') TRIM(textline), " with lines" 1399 WRITE(UNIT=textline,FMT='(a,a)') TRIM(textline), " linecolor 3" 1400 WRITE(UNIT=textline,FMT='(a,a)') TRIM(textline), " linewidth 4" 1401 IF ((nseries>1) .AND. (iseries<nseries)) THEN 1402 WRITE(UNIT=textline,FMT='(a,a)') TRIM(textline), ", \" 1403 END IF 1404 WRITE(UNIT=PLOT_FID,FMT='(a)') TRIM(textline) 1405 END DO 1406 ELSE IF (nseries==1) THEN 1407 DO idatalength=1,ndatalength 1408 IF (idatalength==1) THEN 1409 WRITE(UNIT=textline,FMT='(a,a,a)') "plot """,TRIM(filename_complete(idatalength)),"""" 1410 ELSE 1411 WRITE(UNIT=textline,FMT='(a,a,a)') " """,TRIM(filename_complete(idatalength)),"""" 1412 END IF 1413 IF (PRESENT(series_labels)) THEN 1414 WRITE(UNIT=textline,FMT='(a,a,a,a)') TRIM(textline), " title """,TRIM(series_labels(idatalength)),"""" 1415 END IF 1416 IF (SIZE(llines,1)==1) THEN 1417 linestat=llines(1) 1418 ELSE 1419 linestat=llines(idatalength) 1420 END IF 1421 IF (SIZE(lpoints,1)==1) THEN 1422 pointstat=lpoints(1) 1423 ELSE 1424 pointstat=lpoints(idatalength) 1425 END IF 1426 WRITE(UNIT=textline,FMT='(a,a)') TRIM(textline), " with lines" 1427 WRITE(UNIT=textline,FMT='(a,a)') TRIM(textline), " linecolor 3" 1428 WRITE(UNIT=textline,FMT='(a,a)') TRIM(textline), " linewidth 4" 1429 WRITE(UNIT=PLOT_FID,FMT='(a)') TRIM(textline) 1430 END DO 1431 END IF 1432 1433 WRITE(UNIT=PLOT_FID,FMT='(a)') "" 1434 CLOSE(UNIT=PLOT_FID) 1435 1436 END IF ! loutput_animate 1437 1438 ! 1439 ! plot the data itself, taking into account the continuity type 1440 ! 1441 IF (nseries>1) THEN 1442 DO iseries=1,nseries 1443 IF (debug>=1) THEN 1444 WRITE(UNIT=*,FMT='(a,a,a)') "writing plot data to file """,TRIM(filename_complete(iseries)),"""" 1445 END IF 1446 OPEN(UNIT=PLOT_FID,FILE=TRIM(filename_complete(iseries))) 1447 1448 WRITE(UNIT=PLOT_FID,FMT='(a)') "# X1 X2" 1449 linestat=0 1450 DO i = 1,length 1451 WRITE(UNIT=PLOT_FID,FMT='(e22.15,a,e22.15)') plotdata(i,1)," ",plotdata(i,iseries+1) 1452 END DO 1453 1454 END DO 1455 1456 ELSE IF (nseries==1) THEN 1457 DO idatalength=1,ndatalength 1458 IF (debug>=1) THEN 1459 WRITE(UNIT=*,FMT='(a,a,a)') "writing plot data to file """,TRIM(filename_complete(idatalength)),"""" 1460 END IF 1461 OPEN(UNIT=PLOT_FID,FILE=TRIM(filename_complete(idatalength))) 1462 1463 IF (idatalength==1) THEN 1464 WRITE(UNIT=PLOT_FID,FMT='(a)') "# X1 X2" 1465 ELSE 1466 WRITE(UNIT=PLOT_FID,FMT='(a)') "" 1467 END IF 1468 IF (PRESENT(datalength)) THEN 1469 startrow=SUM(datalength(1:idatalength-1))+1 1470 endrow=startrow+datalength(idatalength)-1 1471 ELSE 1472 startrow=1 1473 endrow=length 1474 END IF 1475 linestat=0 1476 DO i = startrow,endrow 1477 WRITE(UNIT=PLOT_FID,FMT='(e22.15,a,e22.15)') plotdata(i,1)," ",plotdata(i,2) 1478 END DO 1479 1480 END DO 1481 1482 END IF 1483 1484 CLOSE(UNIT=PLOT_FID) 1485 1486 DEALLOCATE(llines) 1487 DEALLOCATE(lpoints) 1488 DEALLOCATE(filename_complete) 1489 1490END SUBROUTINE plotxy 1491 1492 1493