diff --git a/code/DAMASK_spectral.f90 b/code/DAMASK_spectral.f90 index b4c249d16..cafcde85b 100644 --- a/code/DAMASK_spectral.f90 +++ b/code/DAMASK_spectral.f90 @@ -77,10 +77,14 @@ program DAMASK_spectral use math use mesh, only : & - mesh_spectral_getResolution, & - mesh_spectral_getDimension, & - mesh_spectral_getHomogenization - + homog, & + res, & + res1_red, & + mesh_NcpElems, & + wgt, & + geomdim, & + virt_dim + use CPFEM, only: & CPFEM_general, & CPFEM_initAll @@ -96,8 +100,7 @@ program DAMASK_spectral rotation_tol, & itmax,& itmin, & - memory_efficient, & - divergence_correction, & + memory_efficient, & DAMASK_NumThreadsInt, & fftw_planner_flag, & fftw_timelimit @@ -132,10 +135,7 @@ program DAMASK_spectral N_l = 0_pInt, & N_t = 0_pInt, & N_n = 0_pInt, & - N_Fdot = 0_pInt, & - Npoints,& ! number of Fourier points - homog, & ! homogenization scheme used - res1_red ! to store res(1)/2 +1 + N_Fdot = 0_pInt character(len=1024) :: & line @@ -159,11 +159,6 @@ program DAMASK_spectral type(bc_type), allocatable, dimension(:) :: bc - - real(pReal) :: wgt - real(pReal), dimension(3) :: geomdim = 0.0_pReal, virt_dim = 0.0_pReal ! physical dimension of volume element per direction - integer(pInt), dimension(3) :: res = 1_pInt ! resolution (number of Fourier points) in each direction - !-------------------------------------------------------------------------------------------------- ! stress, stiffness and compliance average etc. real(pReal), dimension(3,3) :: & @@ -368,15 +363,6 @@ program DAMASK_spectral end select enddo; enddo 101 close(myUnit) - -!-------------------------------------------------------------------------------------------------- -! get resolution, dimension, homogenization and variables derived from resolution - res = mesh_spectral_getResolution() - geomdim = mesh_spectral_getDimension() - homog = mesh_spectral_getHomogenization() - res1_red = res(1)/2_pInt + 1_pInt ! size of complex array in first dimension (c2r, r2c) - Npoints = res(1)*res(2)*res(3) - wgt = 1.0_pReal/real(Npoints, pReal) !-------------------------------------------------------------------------------------------------- ! output of geometry @@ -541,8 +527,8 @@ program DAMASK_spectral call IO_read_jobBinaryFile(777,'F_aim_lastInc',trim(getSolverJobName()),size(F_aim_lastInc)) read (777,rec=1) F_aim_lastInc close (777) - coordinates = 0.0 ! change it later!!! + CPFEM_mode = 2_pInt endif ielem = 0_pInt do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt, res(1) @@ -561,13 +547,6 @@ program DAMASK_spectral !-------------------------------------------------------------------------------------------------- ! calculation of discrete angular frequencies, ordered as in FFTW (wrap around) - if (divergence_correction) then - do i = 1_pInt, 3_pInt - if (i /= minloc(geomdim,1) .and. i /= maxloc(geomdim,1)) virt_dim = geomdim/geomdim(i) - enddo - else - virt_dim = geomdim - endif do k = 1_pInt, res(3) k_s(3) = k - 1_pInt @@ -632,7 +611,7 @@ program DAMASK_spectral write(538) 'increments', bc(1:N_Loadcases)%incs ! one entry per loadcase write(538) 'startingIncrement', restartInc - 1_pInt ! start with writing out the previous inc write(538) 'eoh' ! end of header - write(538) materialpoint_results(1_pInt:materialpoint_sizeResults,1,1_pInt:Npoints) ! initial (non-deformed or read-in) results + write(538) materialpoint_results(1_pInt:materialpoint_sizeResults,1,1_pInt:mesh_NcpElems) ! initial (non-deformed or read-in) results if (debugGeneral) write(6,'(a)') 'Header of result file written out' endif @@ -661,7 +640,7 @@ program DAMASK_spectral !################################################################################################## do inc = 1_pInt, bc(loadcase)%incs totalIncsCounter = totalIncsCounter + 1_pInt - + !-------------------------------------------------------------------------------------------------- ! forwarding time timeinc_old = timeinc @@ -711,7 +690,7 @@ program DAMASK_spectral F_lastInc(i,j,k,1:3,1:3) = temp33_Real Favg = Favg + F(i,j,k,1:3,1:3) enddo; enddo; enddo - deltaF_aim = guessmode *(Favg*wgt -F_aim) ! average correction in case of guessing to + deltaF_aim = guessmode *(Favg*wgt -math_rotate_backward33(F_aim,bc(loadcase)%rotation)) ! average correction in case of guessing to do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt, res(1) F(i,j,k,1:3,1:3) = F(i,j,k,1:3,1:3) - deltaF_aim ! correct in case avg of F is not F_aim enddo; enddo; enddo @@ -806,9 +785,9 @@ program DAMASK_spectral restartWrite = .false. ! for test of regridding - ! if( bc(loadcase)%restartFrequency > 0_pInt .and. & - ! mod(inc-1,bc(loadcase)%restartFrequency) == 0_pInt .and. & - ! restartInc/=inc) call quit(-1*(restartInc+1)) ! trigger exit to regrid + ! if( bc(loadcase)%restartFrequency > 0_pInt .and. & + ! mod(inc-1,bc(loadcase)%restartFrequency) == 0_pInt .and. & + ! restartInc/=inc) call quit(-1*(restartInc+1)) ! trigger exit to regrid !-------------------------------------------------------------------------------------------------- ! copy one component of the stress field to to a single FT and check for mismatch @@ -981,7 +960,7 @@ program DAMASK_spectral endif deltaF_fourier(1,1,1,1:3,1:3) = cmplx((F_aim_lab_lastIter - F_aim_lab) & ! assign (negative) average deformation gradient change to zero frequency (real part) - * real(Npoints,pReal),0.0_pReal,pReal) ! singular point at xi=(0.0,0.0,0.0) i.e. i=j=k=1 + * real(mesh_NcpElems,pReal),0.0_pReal,pReal) ! singular point at xi=(0.0,0.0,0.0) i.e. i=j=k=1 !-------------------------------------------------------------------------------------------------- ! comparing 1 and 3x3 inverse FT results @@ -1077,7 +1056,7 @@ program DAMASK_spectral if (mod(inc,bc(loadcase)%outputFrequency) == 0_pInt) then ! at output frequency write(6,'(a)') '' write(6,'(a)') '... writing results to file ......................................' - write(538) materialpoint_results(1_pInt:materialpoint_sizeResults,1,1_pInt:Npoints) ! write result to file + write(538) materialpoint_results(1_pInt:materialpoint_sizeResults,1,1_pInt:mesh_NcpElems) ! write result to file flush(538) endif @@ -1092,9 +1071,6 @@ program DAMASK_spectral call IO_write_jobBinaryFile(777,'convergedSpectralDefgrad_lastInc',size(F_lastInc)) ! writing F_lastInc field to file write (777,rec=1) F_lastInc close (777) - call IO_write_jobBinaryFile(777,'C',size(C)) - write (777,rec=1) C - close(777) call IO_write_jobBinaryFile(777,'F_aim',size(F_aim)) write (777,rec=1) F_aim close(777) diff --git a/code/DAMASK_spectral_Driver.f90 b/code/DAMASK_spectral_Driver.f90 index ca078b32d..4ceadb829 100644 --- a/code/DAMASK_spectral_Driver.f90 +++ b/code/DAMASK_spectral_Driver.f90 @@ -123,16 +123,6 @@ program DAMASK_spectral_Driver write(6,'(a)') ' <<<+- DAMASK_spectral_Driver init -+>>>' write(6,'(a)') ' $Id$' #include "compilation_info.f90" - write(6,'(a,a)') ' Working Directory: ',trim(getSolverWorkingDirectoryName()) - write(6,'(a,a)') ' Solver Job Name: ',trim(getSolverJobName()) - write(6,'(a)') '' - write(6,'(a,a)') ' geometry file: ',trim(geometryFile) - write(6,'(a)') '' - write(6,'(a,3(i12 ))') ' resolution a b c:', res - write(6,'(a,3(f12.5))') ' dimension x y z:', geomdim - write(6,'(a,i5)') ' homogenization: ', homog - write(6,'(a,a)') '','' - write(6,'(a,a)') ' Loadcase file: ',trim(loadCaseFile) write(6,'(a)') '' !-------------------------------------------------------------------------------------------------- ! reading basic information from load case file and allocate data structure containing load cases diff --git a/code/DAMASK_spectral_SolverAL.f90 b/code/DAMASK_spectral_SolverAL.f90 index be9a3cf1a..3270abf33 100644 --- a/code/DAMASK_spectral_SolverAL.f90 +++ b/code/DAMASK_spectral_SolverAL.f90 @@ -21,6 +21,19 @@ module DAMASK_spectral_SolverAL solutionState implicit none +#ifdef PETSC +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#endif + character (len=*), parameter, public :: & DAMASK_spectral_SolverAL_label = 'AL' @@ -37,7 +50,7 @@ module DAMASK_spectral_SolverAL ! PETSc data SNES, private :: snes DM, private :: da - Vec, private :: x,r + Vec, private :: x,residual PetscMPIInt, private :: rank integer(pInt), private :: iter PetscInt, private :: xs,xm,gxs,gxm @@ -95,18 +108,20 @@ module DAMASK_spectral_SolverAL use mesh, only: & res, & geomdim - + use math, only: & + math_invSym3333 + implicit none integer(pInt) :: i,j,k - real(pReal), dimension(3,3) :: temp33_Real PetscErrorCode ierr_psc - + PetscObject dummy + call Utilities_init() write(6,'(a)') '' write(6,'(a)') ' <<<+- DAMASK_spectral_solverAL init -+>>>' write(6,'(a)') ' $Id: DAMASK_spectral_SolverAL.f90 1654 2012-08-03 09:25:48Z MPIE\m.diehl $' - #include "compilation_info.f90" +#include "compilation_info.f90" write(6,'(a)') '' allocate (F ( res(1), res(2),res(3),3,3), source = 0.0_pReal) @@ -157,7 +172,7 @@ module DAMASK_spectral_SolverAL coordinates = 0.0 ! change it later!!! endif - call constitutiveResponse(coordinates,F,F_lastInc,temperature,0.0_pReal,& + call Utilities_constitutiveResponse(coordinates,F,F_lastInc,temperature,0.0_pReal,& P,C,P_av,.false.,math_I3) !-------------------------------------------------------------------------------------------------- @@ -187,11 +202,11 @@ module DAMASK_spectral_SolverAL DMDA_STENCIL_BOX,res(1),res(2),res(3),PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE, & 18,1,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,da,ierr_psc) call DMCreateGlobalVector(da,x,ierr_psc) - call VecDuplicate(x,r,ierr_psc) - call DMDASetLocalFunction(da,AL_FormRHS,ierr_psc) + call VecDuplicate(x,residual,ierr_psc) + !call DMDASetLocalFunction(da,AL_FormRHS,ierr_psc) call SNESSetDM(snes,da,ierr_psc) - call SNESSetFunction(snes,r,SNESDMDAComputeFunction,da,ierr_psc) + call SNESSetFunction(snes,residual,AL_FormRHS,dummy,ierr_psc) call SNESSetConvergenceTest(snes,AL_converged,dummy,PETSC_NULL_FUNCTION,ierr_psc) call PetscOptionsInsertString(PetSc_options,ierr_psc) call SNESSetFromOptions(snes,ierr_psc) @@ -239,17 +254,11 @@ module DAMASK_spectral_SolverAL real(pReal), intent(in) :: timeinc, timeinc_old, temperature_bc, guessmode type(boundaryCondition), intent(in) :: P_BC,F_BC real(pReal), dimension(3,3), intent(in) :: rotation_BC - - - + SNESConvergedReason reason real(pReal), dimension(3,3) :: deltaF_aim, & - F_aim_lab, & - F_aim_lab_lastIter + F_aim_lab !-------------------------------------------------------------------------------------------------- ! loop variables, convergence etc. - real(pReal) :: err_div, err_stress - integer(pInt) :: iter, row, column, i, j, k - real(pReal) :: defgradDet, defgradDetMax, defgradDetMin real(pReal), dimension(3,3) :: temp33_Real !-------------------------------------------------------------------------------------------------- @@ -268,7 +277,7 @@ module DAMASK_spectral_SolverAL write (777,rec=1) C close(777) endif - + AL_solution%converged =.false. !-------------------------------------------------------------------------------------------------- ! winding forward of deformation aim in loadcase system if (F_BC%myType=='l') then ! calculate deltaF_aim from given L and current F @@ -306,7 +315,10 @@ module DAMASK_spectral_SolverAL call AL_InitialGuess(xx_psc) call VecRestoreArrayF90(x,xx_psc,ierr_psc) call SNESSolve(snes,PETSC_NULL_OBJECT,x,ierr_psc) - + call SNESGetConvergedReason(snes,reason,ierr_psc) + print*, 'reason', reason + if (reason > 0 ) AL_solution%converged = .true. + end function AL_solution ! ------------------------------------------------------------------- @@ -314,7 +326,7 @@ module DAMASK_spectral_SolverAL subroutine AL_InitialGuess(xx_psc) implicit none - #include + ! Input/output variables: @@ -347,6 +359,55 @@ module DAMASK_spectral_SolverAL return end subroutine AL_InitialGuess + + subroutine AL_FormRHS(snes_local,X_local,F_local,dummy,ierr_psc) + + ! Input/output variables: + SNES snes_local + Vec X_local,F_local + PetscErrorCode ierr_psc + PetscObject dummy + DM da_local + +! Declarations for use with local arrays: + PetscScalar,pointer :: lx_v(:),lf_v(:) + Vec localX + +! Scatter ghost points to local vector, using the 2-step process +! DMGlobalToLocalBegin(), DMGlobalToLocalEnd(). +! By placing code between these two statements, computations can +! be done while messages are in transition. + call SNESGetDM(snes_local,da_local,ierr_psc) + call DMGetLocalVector(da_local,localX,ierr_psc) + call DMGlobalToLocalBegin(da_local,X_local,INSERT_VALUES, & + & localX,ierr_psc) + call DMGlobalToLocalEnd(da_local,X_local,INSERT_VALUES,localX,ierr_psc) + +! Get a pointer to vector data. +! - For default PETSc vectors, VecGetArray90() returns a pointer to +! the data array. Otherwise, the routine is implementation dependent. +! - You MUST call VecRestoreArrayF90() when you no longer need access to +! the array. +! - Note that the interface to VecGetArrayF90() differs from VecGetArray(), +! and is useable from Fortran-90 Only. + + call VecGetArrayF90(localX,lx_v,ierr_psc) + call VecGetArrayF90(F_local,lf_v,ierr_psc) + +! Compute function over the locally owned part of the grid + call AL_FormRHS_local(lx_v,lf_v,dummy,ierr_psc) + +! Restore vectors + call VecRestoreArrayF90(localX,lx_v,ierr_psc) + call VecRestoreArrayF90(F_local,lf_v,ierr_psc) + +! Insert values into global vector + + call DMRestoreLocalVector(da_local,localX,ierr_psc) + + return + + end subroutine AL_FormRHS ! --------------------------------------------------------------------- ! @@ -360,7 +421,7 @@ module DAMASK_spectral_SolverAL ! Notes: ! This routine uses standard Fortran-style computations over a 3-dim array. ! - subroutine AL_FormRHS(in,x_scal,f_scal,dummy,ierr_psc) + subroutine AL_FormRHS_local(x_scal,f_scal,dummy,ierr_psc) use numerics, only: & itmax, & @@ -370,7 +431,8 @@ module DAMASK_spectral_SolverAL math_transpose33, & math_mul3333xx33 use mesh, only: & - res + res, & + wgt use DAMASK_spectral_Utilities, only: & field_real, & Utilities_forwardFFT, & @@ -379,13 +441,10 @@ module DAMASK_spectral_SolverAL Utilities_constitutiveResponse implicit none - #include integer(pInt) :: i,j,k - Input/output variables: - DMDALocalInfo in(DMDA_LOCAL_INFO_SIZE) - PetscScalar x_scal(0:17,XG_RANGE,YG_RANGE,ZG_RANGE) - PetscScalar f_scal(0:17,X_RANGE,Y_RANGE,Z_RANGE) + PetscScalar x_scal(0:17,gxs:gxs+gxm,gys:gys+gym,gzs:gzs+gzm) + PetscScalar f_scal(0:17,xs:xs+xm,ys:ys+ym,zs:zs+zm) real(pReal), dimension (3,3) :: temp33_real PetscObject dummy PetscErrorCode ierr_psc @@ -399,7 +458,7 @@ module DAMASK_spectral_SolverAL write(6,'(3(a,i6.6))') ' @ Iter. ',itmin,' < ',iter,' < ',itmax write(6,'(a,/,3(3(f12.7,1x)/))',advance='no') 'deformation gradient aim =',& math_transpose33(F_aim) - + do k=gzs,gzs+gzm; do j=gys,gys+gym; do i=gxs,gxs+gxm F(i,j,k,1,1) = x_scal(0,i,j,k) F(i,j,k,1,2) = x_scal(1,i,j,k) @@ -420,10 +479,10 @@ module DAMASK_spectral_SolverAL F_lambda(i,j,k,3,2) = x_scal(16,i,j,k) F_lambda(i,j,k,3,3) = x_scal(17,i,j,k) enddo; enddo; enddo - + !-------------------------------------------------------------------------------------------------- ! evaluate constitutive response - call constitutiveResponse(coordinates,F,F_lastInc,temperature,params%timeinc,& + call Utilities_constitutiveResponse(coordinates,F_lastInc,F,temperature,params%timeinc,& P,C,P_av,ForwardData,params%rotation_BC) ForwardData = .False. @@ -432,8 +491,6 @@ module DAMASK_spectral_SolverAL F_aim = F_aim - math_mul3333xx33(S, ((P_av - params%P_BC))) !S = 0.0 for no bc err_stress = maxval(mask_stress * (P_av - params%P_BC)) ! mask = 0.0 for no bc - - F_aim_lab = math_rotate_backward33(F_aim,params%rotation_BC) ! boundary conditions from load frame into lab (Fourier) frame !-------------------------------------------------------------------------------------------------- ! doing Fourier transform @@ -441,27 +498,28 @@ module DAMASK_spectral_SolverAL do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt, res(1) field_real(i,j,k,1:3,1:3) = math_mul3333xx33(C_scale,F_lambda(i,j,k,1:3,1:3)-F(i,j,k,1:3,1:3)) enddo; enddo; enddo + + call Utilities_forwardFFT() - call Utilities_fourierConvolution(F_aim_lab) + call Utilities_fourierConvolution(math_rotate_backward33(F_aim,params%rotation_BC)) call Utilities_backwardFFT() err_f = 0.0_pReal err_p = 0.0_pReal - do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt, res(1) temp33_real = field_real(i,j,k,1:3,1:3) - F(i,j,k,1:3,1:3) err_f = err_f + sum(temp33_real*temp33_real) temp33_real = F_lambda(i,j,k,1:3,1:3) - & - math_mul3333xx33(S_scale,P(i,j,k,1:3,1:3)) + math_I3 + (math_mul3333xx33(S_scale,P(i,j,k,1:3,1:3)) + math_I3) err_p = err_p + sum(temp33_real*temp33_real) - enddo; enddo; enddo - - err_f = wgt*sqrt(err_f)/sum((F_aim-math_I3)*(F_aim-math_I3))) - err_p = wgt*sqrt(err_p)/sum((F_aim-math_I3)*(F_aim-math_I3))) - - do k=zs,ze; do j=ys,ye; do i=xs,xe + enddo; enddo; enddo + err_f = wgt*sqrt(err_f) + err_p = wgt*sqrt(err_p) + + + do k=gzs,gzs+gzm; do j=gys,gys+gym; do i=gxs,gxs+gxm temp33_real = math_mul3333xx33(S_scale,P(i,j,k,1:3,1:3)) + math_I3 - F_lambda(i,j,k,1:3,1:3) & + F(i,j,k,1:3,1:3) - field_real(i,j,k,1:3,1:3) f_scal(0,i,j,k) = temp33_real(1,1) @@ -485,12 +543,12 @@ module DAMASK_spectral_SolverAL enddo; enddo; enddo return - end subroutine AL_FormRHS + end subroutine AL_FormRHS_local ! --------------------------------------------------------------------- ! User defined convergence check ! - subroutine AL_converged(snes,it,xnorm,snorm,fnorm,reason,dummy,ierr_psc) + subroutine AL_converged(snes_local,it,xnorm,snorm,fnorm,reason,dummy,ierr_psc) use numerics, only: & itmax, & @@ -501,10 +559,9 @@ module DAMASK_spectral_SolverAL err_stress_tolabs implicit none - #include ! Input/output variables: - SNES snes + SNES snes_local PetscInt it PetscReal xnorm, snorm, fnorm SNESConvergedReason reason @@ -512,13 +569,15 @@ module DAMASK_spectral_SolverAL PetscErrorCode ierr_psc logical :: Converged - Converged = (iter < itmax) .and. (iter > itmin) .and. & + Converged = (iter > itmin .and. & all([ err_f/sqrt(sum((F_aim-math_I3)*(F_aim-math_I3)))/err_f_tol, & err_p/sqrt(sum((F_aim-math_I3)*(F_aim-math_I3)))/err_p_tol, & - err_stress/min(maxval(abs(P_av))*err_stress_tolrel,err_stress_tolabs)] < 1.0_pReal) + err_stress/min(maxval(abs(P_av))*err_stress_tolrel,err_stress_tolabs)] < 1.0_pReal)) if (Converged) then reason = 1 + elseif (iter > itmax) then + reason = -1 else reason = 0 endif @@ -527,14 +586,17 @@ module DAMASK_spectral_SolverAL write(6,'(a,es14.7)') 'error F = ', err_f/sqrt(sum((F_aim-math_I3)*(F_aim-math_I3)))/err_f_tol write(6,'(a,es14.7)') 'error P = ', err_p/sqrt(sum((F_aim-math_I3)*(F_aim-math_I3)))/err_p_tol return + end subroutine AL_converged subroutine AL_destroy() - + use DAMASK_spectral_Utilities, only: & + Utilities_destroy implicit none + PetscErrorCode ierr_psc call VecDestroy(x,ierr_psc) - call VecDestroy(r,ierr_psc) + call VecDestroy(residual,ierr_psc) call SNESDestroy(snes,ierr_psc) call DMDestroy(da,ierr_psc) call PetscFinalize(ierr_psc) diff --git a/code/DAMASK_spectral_interface.f90 b/code/DAMASK_spectral_interface.f90 index 302a89db2..7e09892db 100644 --- a/code/DAMASK_spectral_interface.f90 +++ b/code/DAMASK_spectral_interface.f90 @@ -178,15 +178,19 @@ subroutine DAMASK_interface_init(loadCaseParameterIn,geometryParameterIn) write(6,'(a,2(i2.2,a),i2.2)') ' Time: ',dateAndTime(5),':',& dateAndTime(6),':',& dateAndTime(7) - write(6,'(a,a)') ' Host name: ', trim(hostName) - write(6,'(a,a)') ' User name: ', trim(userName) - write(6,'(a,a)') ' Path separator: ', getPathSep() - write(6,'(a,a)') ' Command line call: ', trim(commandLine) - write(6,'(a,a)') ' Geometry parameter: ', trim(geometryParameter) - write(6,'(a,a)') ' Loadcase parameter: ', trim(loadcaseParameter) + write(6,'(a,a)') ' Host name: ', trim(hostName) + write(6,'(a,a)') ' User name: ', trim(userName) + write(6,'(a,a)') ' Path separator: ', getPathSep() + write(6,'(a,a)') ' Command line call: ', trim(commandLine) + write(6,'(a,a)') ' Geometry parameter: ', trim(geometryParameter) + write(6,'(a,a)') ' Loadcase parameter: ', trim(loadcaseParameter) + write(6,'(a,a)') ' Geometry file: ', trim(geometryFile) + write(6,'(a,a)') ' Loadcase file: ', trim(loadCaseFile) + write(6,'(a,a)') ' Working Directory: ', trim(getSolverWorkingDirectoryName()) + write(6,'(a,a)') ' Solver Job Name: ', trim(getSolverJobName()) if (SpectralRestartInc > 1_pInt) write(6,'(a,i6.6)') & - ' Restart at increment: ', spectralRestartInc - write(6,'(a,l1)') ' Append to result file: ', appendToOutFile + ' Restart at increment: ', spectralRestartInc + write(6,'(a,l1,/)') ' Append to result file: ', appendToOutFile end subroutine DAMASK_interface_init diff --git a/code/FEsolving.f90 b/code/FEsolving.f90 index aa7d9028e..e81da9537 100644 --- a/code/FEsolving.f90 +++ b/code/FEsolving.f90 @@ -100,9 +100,14 @@ subroutine FE_init character(len=1024) :: line integer(pInt), dimension(1_pInt+2_pInt*maxNchunks) :: positions #endif +!$OMP CRITICAL (write2out) + write(6,*) + write(6,*) '<<<+- FEsolving init -+>>>' + write(6,*) '$Id$' +#include "compilation_info.f90" +!$OMP END CRITICAL (write2out) modelName = getSolverJobName() - #ifdef Spectral restartInc = spectralRestartInc if(restartInc <= 0_pInt) then @@ -169,13 +174,11 @@ subroutine FE_init #endif 200 close(fileunit) endif - + ! the following array are allocated by mesh.f90 and need to be deallocated in case of regridding + if (allocated(calcMode)) deallocate(calcMode) + if (allocated(FEsolving_execIP)) deallocate(FEsolving_execIP) #endif !$OMP CRITICAL (write2out) - write(6,*) - write(6,*) '<<<+- FEsolving init -+>>>' - write(6,*) '$Id$' -#include "compilation_info.f90" if (iand(debug_level(debug_FEsolving),debug_levelBasic) /= 0_pInt) then write(6,*) 'restart writing: ', restartWrite write(6,*) 'restart reading: ', restartRead diff --git a/code/Makefile b/code/Makefile index c79daa086..49a64cc65 100644 --- a/code/Makefile +++ b/code/Makefile @@ -39,10 +39,10 @@ SHELL = /bin/sh ######################################################################################## #auto values will be set by setup_code.py -FFTWROOT := auto -IMKLROOT := auto -ACMLROOT := auto -LAPACKROOT := auto +FFTWROOT := /usr/local +IMKLROOT := +ACMLROOT := +LAPACKROOT := /usr F90 ?= ifort COMPILERNAME ?= $(F90) @@ -96,11 +96,13 @@ endif LIBRARIES +=-lfftw3 LIB_DIRS +=-L$(FFTWROOT)/lib -#ifdef PETSC_DIR -#include ${PETSC_DIR}/conf/variables -#INCLUDE_DIRS +=${PETSC_FC_INCLUDES} -DPETSC -#LIBRARIES +=${PETSC_WITH_EXTERNAL_LIB} -#endif +ifeq "$(SOLVER)" "NEW" +ifdef PETSC_DIR +include ${PETSC_DIR}/conf/variables +INCLUDE_DIRS +=${PETSC_FC_INCLUDES} -DPETSC +LIBRARIES +=${PETSC_WITH_EXTERNAL_LIB} +endif +endif ifdef IMKLROOT LIB_DIRS +=-L$(IMKLROOT)/lib/intel64 @@ -273,10 +275,29 @@ COMPILE_MAXOPTI =$(OPENMP_FLAG_$(F90)) $(COMPILE_OPTIONS_$(F90)) $(STANDARD_CHEC COMPILED_FILES = prec.o DAMASK_spectral_interface.o IO.o numerics.o debug.o math.o \ FEsolving.o mesh.o material.o lattice.o \ constitutive_dislotwin.o constitutive_j2.o constitutive_phenopowerlaw.o \ - constitutive_titanmod.o constitutive_nonlocal.o constitutive_none.o constitutive.o \ - homogenization_RGC.o homogenization_isostrain.o homogenization.o CPFEM.o crystallite.o + constitutive_titanmod.o constitutive_nonlocal.o constitutive_none.o constitutive.o crystallite.o \ + homogenization_RGC.o homogenization_isostrain.o homogenization.o CPFEM.o -ifneq "$(SOLVER)" "AL" +ifeq "$(SOLVER)" "NEW" +COMPILED_FILES += DAMASK_spectral_Utilities.o DAMASK_spectral_SolverBasic.o DAMASK_spectral_SolverAL.o + +DAMASK_spectral.exe: DAMASK_spectral_Driver.o + $(PREFIX) $(COMPILERNAME) $(OPENMP_FLAG_$(F90)) $(OPTIMIZATION_$(MAXOPTI)_$(F90)) $(STANDARD_CHECK_$(F90)) \ + -o DAMASK_spectral.exe DAMASK_spectral_Driver.o \ + $(COMPILED_FILES) $(LIB_DIRS) $(LIBRARIES) $(SUFFIX) + +DAMASK_spectral_Driver.o: DAMASK_spectral_Driver.f90 DAMASK_spectral_SolverBasic.o DAMASK_spectral_SolverAL.o + $(PREFIX) $(COMPILERNAME) $(COMPILE_MAXOPTI) -c DAMASK_spectral_Driver.f90 $(SUFFIX) + +DAMASK_spectral_SolverAL.o: DAMASK_spectral_SolverAL.f90\ + DAMASK_spectral_Utilities.o + +DAMASK_spectral_SolverBasic.o: DAMASK_spectral_SolverBasic.f90\ + DAMASK_spectral_Utilities.o + +DAMASK_spectral_Utilities.o: DAMASK_spectral_Utilities.f90\ + CPFEM.o +else DAMASK_spectral.exe: DAMASK_spectral.o $(PREFIX) $(COMPILERNAME) $(OPENMP_FLAG_$(F90)) $(OPTIMIZATION_$(MAXOPTI)_$(F90)) $(STANDARD_CHECK_$(F90)) \ -o DAMASK_spectral.exe DAMASK_spectral.o \ @@ -284,15 +305,7 @@ DAMASK_spectral.exe: DAMASK_spectral.o DAMASK_spectral.o: DAMASK_spectral.f90 CPFEM.o $(PREFIX) $(COMPILERNAME) $(COMPILE_MAXOPTI) -c DAMASK_spectral.f90 $(SUFFIX) -else -DAMASK_spectral_AL.exe: DAMASK_spectral_AL.o - $(PREFIX) $(COMPILERNAME) $(OPENMP_FLAG_$(F90)) $(OPTIMIZATION_$(MAXOPTI)_$(F90)) $(STANDARD_CHECK_$(F90)) \ - -o DAMASK_spectral_AL.exe DAMASK_spectral_AL.o \ - $(COMPILED_FILES) $(LIB_DIRS) $(LIBRARIES) $(SUFFIX) - -DAMASK_spectral_AL.o: DAMASK_spectral_AL.f90 CPFEM.o - $(PREFIX) $(COMPILERNAME) $(COMPILE_MAXOPTI) -c DAMASK_spectral_AL.f90 $(SUFFIX) -endif +endif CPFEM.o: CPFEM.f90\ homogenization.o diff --git a/code/mesh.f90 b/code/mesh.f90 index 48bff3c06..b93d50d69 100644 --- a/code/mesh.f90 +++ b/code/mesh.f90 @@ -29,7 +29,7 @@ module mesh use prec, only: pReal, pInt - + implicit none private @@ -94,7 +94,14 @@ module mesh mesh_subNodeCoord !< coordinates of subnodes per element logical, private :: noPart !< for cases where the ABAQUS input file does not use part/assembly information - + +#ifdef Spectral + real(pReal), dimension(3), public :: geomdim, virt_dim ! physical dimension of volume element per direction + integer(pInt), dimension(3), public :: res + real(pReal), public :: wgt + integer(pInt), public :: res1_red, homog + integer(pInt), private :: i +#endif ! These definitions should actually reside in the FE-solver specific part (different for MARC/ABAQUS) ! Hence, I suggest to prefix with "FE_" @@ -272,14 +279,14 @@ module mesh mesh_build_ipVolumes, & mesh_build_ipCoordinates #ifdef Spectral - public :: mesh_spectral_getResolution, & - mesh_spectral_getDimension, & - mesh_spectral_getHomogenization, & - mesh_regrid + public :: mesh_regrid #endif private :: & #ifdef Spectral + mesh_spectral_getResolution, & + mesh_spectral_getDimension, & + mesh_spectral_getHomogenization, & mesh_spectral_count_nodesAndElements, & mesh_spectral_count_cpElements, & mesh_spectral_map_elements, & @@ -340,6 +347,8 @@ subroutine mesh_init(ip,element) #endif #ifdef Spectral IO_open_file + use numerics, only: & + divergence_correction #else IO_open_InputFile #endif @@ -375,8 +384,6 @@ subroutine mesh_init(ip,element) if (allocated(mesh_ipNeighborhood)) deallocate(mesh_ipNeighborhood) if (allocated(mesh_ipVolume)) deallocate(mesh_ipVolume) if (allocated(mesh_nodeTwins)) deallocate(mesh_nodeTwins) - if (allocated(calcMode)) deallocate(calcMode) - if (allocated(FEsolving_execIP)) deallocate(FEsolving_execIP) if (allocated(FE_nodesAtIP)) deallocate(FE_nodesAtIP) if (allocated(FE_ipNeighbor))deallocate(FE_ipNeighbor) if (allocated(FE_subNodeParent)) deallocate(FE_subNodeParent) @@ -384,12 +391,27 @@ subroutine mesh_init(ip,element) call mesh_build_FEdata ! get properties of the different types of elements #ifdef Spectral call IO_open_file(fileUnit,geometryFile) ! parse info from geometry file... - call mesh_spectral_count_nodesAndElements(fileUnit) + res = mesh_spectral_getResolution(fileUnit) + res1_red = res(1)/2_pInt + 1_pInt + wgt = 1.0/real(res(1)*res(2)*res(3),pReal) + geomdim = mesh_spectral_getDimension(fileUnit) + homog = mesh_spectral_getHomogenization(fileUnit) + if (divergence_correction) then + do i = 1_pInt, 3_pInt + if (i /= minloc(geomdim,1) .and. i /= maxloc(geomdim,1)) virt_dim = geomdim/geomdim(i) + enddo + else + virt_dim = geomdim + endif + write(6,'(a,3(i12 ))') ' resolution a b c:', res + write(6,'(a,3(f12.5))') ' dimension x y z:', geomdim + write(6,'(a,i5,/)') ' homogenization: ', homog + call mesh_spectral_count_nodesAndElements call mesh_spectral_count_cpElements call mesh_spectral_map_elements call mesh_spectral_map_nodes call mesh_spectral_count_cpSizes - call mesh_spectral_build_nodes(fileUnit) + call mesh_spectral_build_nodes call mesh_spectral_build_elements(fileUnit) #endif #ifdef Marc @@ -864,10 +886,10 @@ function mesh_regrid(adaptive,resNewInput,minRes) integer(pInt), dimension(3), optional, intent(in) :: minRes integer(pInt), dimension(3) :: mesh_regrid, ratio integer(pInt), dimension(3,2) :: possibleResNew - integer(pInt):: maxsize, i, j, k, ielem, Npoints, NpointsNew, spatialDim - integer(pInt), dimension(3) :: res, resNew + integer(pInt):: maxsize, i, j, k, ielem, NpointsNew, spatialDim + integer(pInt), dimension(3) :: resNew integer(pInt), dimension(:), allocatable :: indices - real(pReal), dimension(3) :: geomdim,geomdimNew + real(pReal), dimension(3) :: geomdimNew real(pReal), dimension(3,3) :: Favg, Favg_LastInc real(pReal), dimension(:,:), allocatable :: & coordinatesNew, & @@ -878,10 +900,10 @@ function mesh_regrid(adaptive,resNewInput,minRes) real(pReal), dimension(:,:,:,:,:), allocatable :: & F, FNew, & - F_lastInc, F_lastIncNew, & Fp, FpNew, & Lp, LpNew, & - dcsdE, dcsdENew + dcsdE, dcsdENew, & + F_lastIncNew real(pReal), dimension (:,:,:,:,:,:,:), allocatable :: & dPdF, dPdFNew @@ -892,10 +914,6 @@ function mesh_regrid(adaptive,resNewInput,minRes) stateConst integer(pInt), dimension(:,:), allocatable :: & sizeStateConst, sizeStateHomog - - res = mesh_spectral_getResolution() - geomdim = mesh_spectral_getDimension() - Npoints = res(1)*res(2)*res(3) if (adaptive) then if (present(resNewInput)) then @@ -911,16 +929,17 @@ function mesh_regrid(adaptive,resNewInput,minRes) read (777,rec=1) F close (777) -! ----Calculate deformed configuration and average-------- - do i= 1_pInt,3_pInt; do j = 1_pInt,3_pInt - Favg(i,j) = sum(F(1:res(1),1:res(2),1:res(3),i,j)) / real(Npoints,pReal) - enddo; enddo +! ----read in average deformation-------- + call IO_read_jobBinaryFile(777,'F_aim',trim(getSolverJobName()),size(Favg)) + read (777,rec=1) Favg + close (777) + allocate(coordinates(res(1),res(2),res(3),3)) call deformed_fft(res,geomdim,Favg,1.0_pReal,F,coordinates) deallocate(F) ! ----Store coordinates into a linear list-------- - allocate(coordinatesLinear(3,Npoints)) + allocate(coordinatesLinear(3,mesh_NcpElems)) ielem = 0_pInt do k=1_pInt,res(3); do j=1_pInt, res(2); do i=1_pInt, res(1) ielem = ielem + 1_pInt @@ -954,7 +973,6 @@ function mesh_regrid(adaptive,resNewInput,minRes) ratio = floor(real(resNewInput,pReal) * (geomdimNew/geomdim), pInt) possibleResNew = 1_pInt - do i = 1_pInt, spatialDim if (mod(ratio(i),2) == 0_pInt) then possibleResNew(i,1:2) = [ratio(i),ratio(i) + 2_pInt] @@ -1016,7 +1034,7 @@ function mesh_regrid(adaptive,resNewInput,minRes) !----- Nearest neighbour search ------------------------------------ allocate(indices(NpointsNew)) - call math_nearestNeighborSearch(spatialDim, Favg, geomdim, NpointsNew, Npoints, & + call math_nearestNeighborSearch(spatialDim, Favg, geomdim, NpointsNew, mesh_NcpElems, & coordinatesNew, coordinatesLinear, indices) deallocate(coordinatesNew) @@ -1045,7 +1063,7 @@ function mesh_regrid(adaptive,resNewInput,minRes) call IO_write_jobFile(777,'idx2') ! make it a general open-write file write(777, '(A)') '1 header' - write(777, '(A)') 'Numbered indices as per the large set' + write(777, '(A)') 'Numbered indices as per the small set' do i = 1_pInt, NpointsNew write(777,trim(formatString),advance='no') indices(i), ' ' if(mod(i,resNew(1)) == 0_pInt) write(777,'(A)') '' @@ -1053,7 +1071,7 @@ function mesh_regrid(adaptive,resNewInput,minRes) close(777) !------Adjusting the point-to-grain association--------------------- - write(N_Digits, '(I16.16)') 1_pInt+int(log10(real(maxval(mesh_element(4,1:Npoints)),pReal)),pInt) + write(N_Digits, '(I16.16)') 1_pInt+int(log10(real(maxval(mesh_element(4,1:mesh_NcpElems)),pReal)),pInt) N_Digits = adjustl(N_Digits) formatString = '(I'//trim(N_Digits)//'.'//trim(N_Digits)//',a)' @@ -1069,25 +1087,24 @@ function mesh_regrid(adaptive,resNewInput,minRes) enddo close(777) -!--------------------------------------------------------- - allocate(F_lastInc(res(1),res(2),res(3),3,3)) - allocate(F_lastIncNew(resNew(1),resNew(2),resNew(3),3,3)) - call IO_read_jobBinaryFile(777,'convergedSpectralDefgrad_lastInc',trim(getSolverJobName()),size(F_lastInc)) - read (777,rec=1) F_lastInc +!---set F_lastInc to homogeneous deformation------------------------------------------------------ + + call IO_read_jobBinaryFile(777,'F_aim_lastInc',trim(getSolverJobName()),size(Favg_LastInc)) + read (777,rec=1) Favg_LastInc close (777) - do i= 1_pInt,3_pInt; do j = 1_pInt,3_pInt - Favg_LastInc(i,j) = sum(F_lastInc(1:res(1),1:res(2),1:res(3),i,j)) / real(Npoints,pReal) - enddo; enddo + + allocate(F_lastIncNew(resNew(1),resNew(2),resNew(3),3,3)) do k=1_pInt,resNew(3); do j=1_pInt, resNew(2); do i=1_pInt, resNew(1) - F_lastIncNew(i,j,k,1:3,1:3) = Favg + F_lastIncNew(i,j,k,1:3,1:3) = Favg_LastInc enddo; enddo; enddo + call IO_write_jobBinaryFile(777,'convergedSpectralDefgrad_lastInc',size(F_lastIncNew)) write (777,rec=1) F_lastIncNew close (777) - deallocate(F_lastInc) deallocate(F_lastIncNew) + !------------------------------------------------------------------- - allocate(material_phase (1,1, Npoints)) + allocate(material_phase (1,1, mesh_NcpElems)) allocate(material_phaseNew (1,1, NpointsNew)) call IO_read_jobBinaryFile(777,'recordedPhase',trim(getSolverJobName()),size(material_phase)) read (777,rec=1) material_phase @@ -1103,7 +1120,7 @@ function mesh_regrid(adaptive,resNewInput,minRes) deallocate(material_phaseNew) !--------------------------------------------------------------------------- - allocate(F (3,3,1,1, Npoints)) + allocate(F (3,3,1,1, mesh_NcpElems)) allocate(FNew (3,3,1,1, NpointsNew)) call IO_read_jobBinaryFile(777,'convergedF',trim(getSolverJobName()),size(F)) read (777,rec=1) F @@ -1118,7 +1135,7 @@ function mesh_regrid(adaptive,resNewInput,minRes) deallocate(F) deallocate(FNew) !--------------------------------------------------------------------- - allocate(Fp (3,3,1,1,Npoints)) + allocate(Fp (3,3,1,1,mesh_NcpElems)) allocate(FpNew (3,3,1,1,NpointsNew)) call IO_read_jobBinaryFile(777,'convergedFp',trim(getSolverJobName()),size(Fp)) read (777,rec=1) Fp @@ -1134,7 +1151,7 @@ function mesh_regrid(adaptive,resNewInput,minRes) deallocate(FpNew) !------------------------------------------------------------------------ - allocate(Lp (3,3,1,1,Npoints)) + allocate(Lp (3,3,1,1,mesh_NcpElems)) allocate(LpNew (3,3,1,1,NpointsNew)) call IO_read_jobBinaryFile(777,'convergedLp',trim(getSolverJobName()),size(Lp)) read (777,rec=1) Lp @@ -1149,7 +1166,7 @@ function mesh_regrid(adaptive,resNewInput,minRes) deallocate(LpNew) !---------------------------------------------------------------------------- - allocate(dcsdE (6,6,1,1,Npoints)) + allocate(dcsdE (6,6,1,1,mesh_NcpElems)) allocate(dcsdENew (6,6,1,1,NpointsNew)) call IO_read_jobBinaryFile(777,'convergeddcsdE',trim(getSolverJobName()),size(dcsdE)) read (777,rec=1) dcsdE @@ -1164,7 +1181,7 @@ function mesh_regrid(adaptive,resNewInput,minRes) deallocate(dcsdENew) !--------------------------------------------------------------------------- - allocate(dPdF (3,3,3,3,1,1,Npoints)) + allocate(dPdF (3,3,3,3,1,1,mesh_NcpElems)) allocate(dPdFNew (3,3,3,3,1,1,NpointsNew)) call IO_read_jobBinaryFile(777,'convergeddPdF',trim(getSolverJobName()),size(dPdF)) read (777,rec=1) dPdF @@ -1179,7 +1196,7 @@ function mesh_regrid(adaptive,resNewInput,minRes) deallocate(dPdFNew) !--------------------------------------------------------------------------- - allocate(Tstar (6,1,1,Npoints)) + allocate(Tstar (6,1,1,mesh_NcpElems)) allocate(TstarNew (6,1,1,NpointsNew)) call IO_read_jobBinaryFile(777,'convergedTstar',trim(getSolverJobName()),size(Tstar)) read (777,rec=1) Tstar @@ -1194,16 +1211,16 @@ function mesh_regrid(adaptive,resNewInput,minRes) deallocate(TstarNew) !---------------------------------------------------------------------------- - allocate(sizeStateConst(1,Npoints)) + allocate(sizeStateConst(1,mesh_NcpElems)) call IO_read_jobBinaryFile(777,'sizeStateConst',trim(getSolverJobName()),size(sizeStateConst)) read (777,rec=1) sizeStateConst close (777) - maxsize = maxval(sizeStateConst(1,1:Npoints)) - allocate(StateConst (1,1,Npoints,maxsize)) + maxsize = maxval(sizeStateConst(1,1:mesh_NcpElems)) + allocate(StateConst (1,1,mesh_NcpElems,maxsize)) call IO_read_jobBinaryFile(777,'convergedStateConst',trim(getSolverJobName())) k = 0_pInt - do i =1, Npoints + do i =1, mesh_NcpElems do j = 1,sizeStateConst(1,i) k = k+1_pInt read(777,rec=k) StateConst(1,1,i,j) @@ -1223,16 +1240,16 @@ function mesh_regrid(adaptive,resNewInput,minRes) deallocate(StateConst) !---------------------------------------------------------------------------- - allocate(sizeStateHomog(1,Npoints)) + allocate(sizeStateHomog(1,mesh_NcpElems)) call IO_read_jobBinaryFile(777,'sizeStateHomog',trim(getSolverJobName()),size(sizeStateHomog)) read (777,rec=1) sizeStateHomog close (777) - maxsize = maxval(sizeStateHomog(1,1:Npoints)) - allocate(stateHomog (1,1,Npoints,maxsize)) + maxsize = maxval(sizeStateHomog(1,1:mesh_NcpElems)) + allocate(stateHomog (1,1,mesh_NcpElems,maxsize)) call IO_read_jobBinaryFile(777,'convergedStateHomog',trim(getSolverJobName())) k = 0_pInt - do i =1, Npoints + do i =1, mesh_NcpElems do j = 1,sizeStateHomog(1,i) k = k+1_pInt read(777,rec=k) stateHomog(1,1,i,j) @@ -1261,13 +1278,9 @@ end function mesh_regrid !> @brief Count overall number of nodes and elements in mesh and stores them in !! 'mesh_Nelems' and 'mesh_Nnodes' !-------------------------------------------------------------------------------------------------- -subroutine mesh_spectral_count_nodesAndElements(myUnit) +subroutine mesh_spectral_count_nodesAndElements() implicit none - integer(pInt), intent(in) :: myUnit - integer(pInt), dimension(3) :: res - - res = mesh_spectral_getResolution(myUnit) mesh_Nelems = res(1)*res(2)*res(3) mesh_Nnodes = (1_pInt + res(1))*(1_pInt + res(2))*(1_pInt + res(3)) @@ -1344,22 +1357,16 @@ end subroutine mesh_spectral_count_cpSizes !> @brief Store x,y,z coordinates of all nodes in mesh. !! Allocates global arrays 'mesh_node0' and 'mesh_node' !-------------------------------------------------------------------------------------------------- -subroutine mesh_spectral_build_nodes(myUnit) +subroutine mesh_spectral_build_nodes() use IO, only: & IO_error implicit none - integer(pInt), intent(in) :: myUnit integer(pInt) :: n - integer(pInt), dimension(3) :: res = 1_pInt - real(pReal), dimension(3) :: geomdim = 1.0_pReal allocate ( mesh_node0 (3,mesh_Nnodes) ); mesh_node0 = 0.0_pReal allocate ( mesh_node (3,mesh_Nnodes) ); mesh_node = 0.0_pReal - - res = mesh_spectral_getResolution(myUnit) - geomdim = mesh_spectral_getDimension(myUnit) forall (n = 0_pInt:mesh_Nnodes-1_pInt) mesh_node0(1,n+1_pInt) = & @@ -1399,16 +1406,12 @@ subroutine mesh_spectral_build_elements(myUnit) integer(pInt), parameter :: maxNchunks = 7_pInt integer(pInt), dimension (1_pInt+2_pInt*maxNchunks) :: myPos - integer(pInt), dimension(3) :: res - integer(pInt) :: e, i, homog = 0_pInt, headerLength = 0_pInt, maxIntCount + integer(pInt) :: e, i, headerLength = 0_pInt, maxIntCount integer(pInt), dimension(:), allocatable :: microstructures integer(pInt), dimension(1,1) :: dummySet = 0_pInt character(len=65536) :: line,keyword character(len=64), dimension(1) :: dummyName = '' - res = mesh_spectral_getResolution(myUnit) - homog = mesh_spectral_getHomogenization(myUnit) - call IO_checkAndRewind(myUnit) read(myUnit,'(a65536)') line diff --git a/code/numerics.f90 b/code/numerics.f90 index df6c44885..85a45019a 100644 --- a/code/numerics.f90 +++ b/code/numerics.f90 @@ -419,7 +419,7 @@ subroutine numerics_init .not. memory_efficient) call IO_error(error_ID = 847_pInt) #endif if (fixedSeed <= 0_pInt) then - write(6,'(a)') ' Random is random!' + write(6,'(a,/)') ' Random is random!' endif end subroutine numerics_init