diff --git a/VERSION b/VERSION index 6a2d0648d..a1ebc5acf 100644 --- a/VERSION +++ b/VERSION @@ -1 +1 @@ -v2.0.1-249-g54c0ba5 +v2.0.1-257-g68006f0 diff --git a/code/DAMASK_spectral.f90 b/code/DAMASK_spectral.f90 index 67e3f4042..6f1794e35 100644 --- a/code/DAMASK_spectral.f90 +++ b/code/DAMASK_spectral.f90 @@ -192,7 +192,7 @@ program DAMASK_spectral if ((N_def /= N_n) .or. (N_n /= N_t) .or. N_n < 1_pInt) & ! sanity check call IO_error(error_ID=837_pInt,ext_msg = trim(loadCaseFile)) ! error message for incomplete loadcase allocate (loadCases(N_n)) ! array of load cases - loadCases%P%myType='p' + loadCases%stress%myType='stress' do i = 1, size(loadCases) allocate(loadCases(i)%ID(nActiveFields)) @@ -244,10 +244,10 @@ program DAMASK_spectral temp_maskVector(j) = IO_stringValue(line,chunkPos,i+j) /= '*' ! true if not an asterisk if (temp_maskVector(j)) temp_valueVector(j) = IO_floatValue(line,chunkPos,i+j) ! read value where applicable enddo - loadCases(currentLoadCase)%P%maskLogical = transpose(reshape(temp_maskVector,[ 3,3])) - loadCases(currentLoadCase)%P%maskFloat = merge(ones,zeros,& - loadCases(currentLoadCase)%P%maskLogical) - loadCases(currentLoadCase)%P%values = math_plain9to33(temp_valueVector) + loadCases(currentLoadCase)%stress%maskLogical = transpose(reshape(temp_maskVector,[ 3,3])) + loadCases(currentLoadCase)%stress%maskFloat = merge(ones,zeros,& + loadCases(currentLoadCase)%stress%maskLogical) + loadCases(currentLoadCase)%stress%values = math_plain9to33(temp_valueVector) case('t','time','delta') ! increment time loadCases(currentLoadCase)%time = IO_floatValue(line,chunkPos,i+1_pInt) case('n','incs','increments','steps') ! number of increments @@ -318,16 +318,16 @@ program DAMASK_spectral endif enddo; write(6,'(/)',advance='no') enddo - if (any(loadCases(currentLoadCase)%P%maskLogical .eqv. & + if (any(loadCases(currentLoadCase)%stress%maskLogical .eqv. & loadCases(currentLoadCase)%deformation%maskLogical)) errorID = 831_pInt ! exclusive or masking only - if (any(loadCases(currentLoadCase)%P%maskLogical .and. & - transpose(loadCases(currentLoadCase)%P%maskLogical) .and. & + if (any(loadCases(currentLoadCase)%stress%maskLogical .and. & + transpose(loadCases(currentLoadCase)%stress%maskLogical) .and. & reshape([ .false.,.true.,.true.,.true.,.false.,.true.,.true.,.true.,.false.],[ 3,3]))) & errorID = 838_pInt ! no rotation is allowed by stress BC write(6,'(2x,a)') 'stress / GPa:' do i = 1_pInt, 3_pInt; do j = 1_pInt, 3_pInt - if(loadCases(currentLoadCase)%P%maskLogical(i,j)) then - write(6,'(2x,f12.7)',advance='no') loadCases(currentLoadCase)%P%values(i,j)*1e-9_pReal + if(loadCases(currentLoadCase)%stress%maskLogical(i,j)) then + write(6,'(2x,f12.7)',advance='no') loadCases(currentLoadCase)%stress%values(i,j)*1e-9_pReal else write(6,'(2x,12a)',advance='no') ' * ' endif @@ -524,30 +524,25 @@ program DAMASK_spectral case (DAMASK_spectral_SolverBasicPETSc_label) call BasicPETSc_forward (& guess,timeinc,timeIncOld,remainingLoadCaseTime, & - F_BC = loadCases(currentLoadCase)%deformation, & - P_BC = loadCases(currentLoadCase)%P, & + deformation_BC = loadCases(currentLoadCase)%deformation, & + stress_BC = loadCases(currentLoadCase)%stress, & rotation_BC = loadCases(currentLoadCase)%rotation) case (DAMASK_spectral_SolverAL_label) call AL_forward (& guess,timeinc,timeIncOld,remainingLoadCaseTime, & - F_BC = loadCases(currentLoadCase)%deformation, & - P_BC = loadCases(currentLoadCase)%P, & + deformation_BC = loadCases(currentLoadCase)%deformation, & + stress_BC = loadCases(currentLoadCase)%stress, & rotation_BC = loadCases(currentLoadCase)%rotation) case (DAMASK_spectral_SolverPolarisation_label) call Polarisation_forward (& guess,timeinc,timeIncOld,remainingLoadCaseTime, & - F_BC = loadCases(currentLoadCase)%deformation, & - P_BC = loadCases(currentLoadCase)%P, & + deformation_BC = loadCases(currentLoadCase)%deformation, & + stress_BC = loadCases(currentLoadCase)%stress, & rotation_BC = loadCases(currentLoadCase)%rotation) end select - case(FIELD_THERMAL_ID) - call spectral_thermal_forward (& - guess,timeinc,timeIncOld,remainingLoadCaseTime) - - case(FIELD_DAMAGE_ID) - call spectral_damage_forward (& - guess,timeinc,timeIncOld,remainingLoadCaseTime) + case(FIELD_THERMAL_ID); call spectral_thermal_forward() + case(FIELD_DAMAGE_ID); call spectral_damage_forward() end select enddo @@ -562,34 +557,29 @@ program DAMASK_spectral select case (spectral_solver) case (DAMASK_spectral_SolverBasicPETSc_label) solres(field) = BasicPETSC_solution (& - incInfo,guess,timeinc,timeIncOld,remainingLoadCaseTime, & - P_BC = loadCases(currentLoadCase)%P, & - F_BC = loadCases(currentLoadCase)%deformation, & + incInfo,timeinc,timeIncOld, & + stress_BC = loadCases(currentLoadCase)%stress, & rotation_BC = loadCases(currentLoadCase)%rotation) case (DAMASK_spectral_SolverAL_label) solres(field) = AL_solution (& - incInfo,guess,timeinc,timeIncOld,remainingLoadCaseTime, & - P_BC = loadCases(currentLoadCase)%P, & - F_BC = loadCases(currentLoadCase)%deformation, & + incInfo,timeinc,timeIncOld, & + stress_BC = loadCases(currentLoadCase)%stress, & rotation_BC = loadCases(currentLoadCase)%rotation) case (DAMASK_spectral_SolverPolarisation_label) solres(field) = Polarisation_solution (& - incInfo,guess,timeinc,timeIncOld,remainingLoadCaseTime, & - P_BC = loadCases(currentLoadCase)%P, & - F_BC = loadCases(currentLoadCase)%deformation, & + incInfo,timeinc,timeIncOld, & + stress_BC = loadCases(currentLoadCase)%stress, & rotation_BC = loadCases(currentLoadCase)%rotation) end select case(FIELD_THERMAL_ID) - solres(field) = spectral_thermal_solution (& - guess,timeinc,timeIncOld,remainingLoadCaseTime) + solres(field) = spectral_thermal_solution(timeinc,timeIncOld,remainingLoadCaseTime) case(FIELD_DAMAGE_ID) - solres(field) = spectral_damage_solution (& - guess,timeinc,timeIncOld,remainingLoadCaseTime) + solres(field) = spectral_damage_solution(timeinc,timeIncOld,remainingLoadCaseTime) end select if (.not. solres(field)%converged) exit ! no solution found diff --git a/code/Makefile b/code/Makefile index 5a43040aa..a3d50d6ef 100644 --- a/code/Makefile +++ b/code/Makefile @@ -22,7 +22,7 @@ DAMASKVERSION :=$(shell cat ../VERSION) include ${PETSC_DIR}/lib/petsc/conf/variables include ${PETSC_DIR}/lib/petsc/conf/rules -INCLUDE_DIRS := $(PETSC_FC_INCLUDES) -DPETSc -I../lib +INCLUDE_DIRS := $(PETSC_FC_INCLUDES) -DPETSc LIBRARIES := $(PETSC_WITH_EXTERNAL_LIB) FCOMPILERNAME ?= $(FC) CCOMPILERNAME ?= $(CC) diff --git a/code/kinematics_hydrogen_strain.f90 b/code/kinematics_hydrogen_strain.f90 index 154b97e79..7a33d1a5f 100644 --- a/code/kinematics_hydrogen_strain.f90 +++ b/code/kinematics_hydrogen_strain.f90 @@ -67,8 +67,6 @@ subroutine kinematics_hydrogen_strain_init(fileUnit) KINEMATICS_hydrogen_strain_ID, & material_Nphase, & MATERIAL_partPhase - use numerics,only: & - worldrank implicit none integer(pInt), intent(in) :: fileUnit @@ -79,11 +77,9 @@ subroutine kinematics_hydrogen_strain_init(fileUnit) tag = '', & line = '' - mainProcess: if (worldrank == 0) then - write(6,'(/,a)') ' <<<+- kinematics_'//KINEMATICS_hydrogen_strain_LABEL//' init -+>>>' - write(6,'(a15,a)') ' Current time: ',IO_timeStamp() + write(6,'(/,a)') ' <<<+- kinematics_'//KINEMATICS_hydrogen_strain_LABEL//' init -+>>>' + write(6,'(a15,a)') ' Current time: ',IO_timeStamp() #include "compilation_info.f90" - endif mainProcess maxNinstance = int(count(phase_kinematics == KINEMATICS_hydrogen_strain_ID),pInt) if (maxNinstance == 0_pInt) return diff --git a/code/mesh.f90 b/code/mesh.f90 index 32f94d66b..7da3d4968 100644 --- a/code/mesh.f90 +++ b/code/mesh.f90 @@ -549,13 +549,13 @@ subroutine mesh_init(ip,el) call IO_open_file(FILEUNIT,geometryFile) ! parse info from geometry file... if (myDebug) write(6,'(a)') ' Opened geometry file'; flush(6) grid = mesh_spectral_getGrid(fileUnit) - call MPI_comm_size(MPI_COMM_WORLD, worldsize, ierr) + call MPI_comm_size(PETSC_COMM_WORLD, worldsize, ierr) if(ierr /=0_pInt) call IO_error(894_pInt, ext_msg='MPI_comm_size') if(worldsize>grid(3)) call IO_error(894_pInt, ext_msg='number of processes exceeds grid(3)') geomSize = mesh_spectral_getSize(fileUnit) devNull = fftw_mpi_local_size_3d(int(grid(3),C_INTPTR_T),int(grid(2),C_INTPTR_T),& - int(grid(1),C_INTPTR_T)/2+1,MPI_COMM_WORLD,local_K,local_K_offset) + int(grid(1),C_INTPTR_T)/2+1,PETSC_COMM_WORLD,local_K,local_K_offset) grid3 = int(local_K,pInt) grid3Offset = int(local_K_offset,pInt) size3 = geomSize(3)*real(grid3,pReal) /real(grid(3),pReal) diff --git a/code/porosity_none.f90 b/code/porosity_none.f90 index 92d4d9fe5..1e6ea9dc9 100644 --- a/code/porosity_none.f90 +++ b/code/porosity_none.f90 @@ -23,21 +23,17 @@ subroutine porosity_none_init() use IO, only: & IO_timeStamp use material - use numerics, only: & - worldrank implicit none integer(pInt) :: & homog, & NofMyHomog - mainProcess: if (worldrank == 0) then - write(6,'(/,a)') ' <<<+- porosity_'//POROSITY_none_label//' init -+>>>' - write(6,'(a15,a)') ' Current time: ',IO_timeStamp() + write(6,'(/,a)') ' <<<+- porosity_'//POROSITY_none_label//' init -+>>>' + write(6,'(a15,a)') ' Current time: ',IO_timeStamp() #include "compilation_info.f90" - endif mainProcess - initializeInstances: do homog = 1_pInt, material_Nhomogenization + initializeInstances: do homog = 1_pInt, material_Nhomogenization myhomog: if (porosity_type(homog) == POROSITY_none_ID) then NofMyHomog = count(material_homog == homog) diff --git a/code/prec.f90 b/code/prec.f90 index 2b7ef8ab9..c6ae656b8 100644 --- a/code/prec.f90 +++ b/code/prec.f90 @@ -20,12 +20,7 @@ module prec private #if (FLOAT==8) integer, parameter, public :: pReal = 8 !< floating point double precision (was selected_real_kind(15,300), number with 15 significant digits, up to 1e+-300) -#ifdef __INTEL_COMPILER - real(pReal), parameter, public :: DAMASK_NaN = Z'7FF8000000000000' !< quiet NaN for double precision (from http://www.hpc.unimelb.edu.au/doc/f90lrm/dfum_035.html) -#endif -#ifdef __GFORTRAN__ real(pReal), parameter, public :: DAMASK_NaN = real(Z'7FF8000000000000',pReal) !< quiet NaN for double precision (from http://www.hpc.unimelb.edu.au/doc/f90lrm/dfum_035.html) -#endif #else NO SUITABLE PRECISION FOR REAL SELECTED, STOPPING COMPILATION #endif diff --git a/code/spectral_damage.f90 b/code/spectral_damage.f90 index 6260fe43a..9c259a2fb 100644 --- a/code/spectral_damage.f90 +++ b/code/spectral_damage.f90 @@ -81,7 +81,6 @@ subroutine spectral_damage_init() DM :: damage_grid Vec :: uBound, lBound PetscErrorCode :: ierr - PetscObject :: dummy character(len=100) :: snes_type external :: & @@ -99,11 +98,9 @@ subroutine spectral_damage_init() DMRestoreGlobalVector, & SNESVISetVariableBounds - mainProcess: if (worldrank == 0_pInt) then - write(6,'(/,a)') ' <<<+- spectral_damage init -+>>>' - write(6,'(a15,a)') ' Current time: ',IO_timeStamp() + write(6,'(/,a)') ' <<<+- spectral_damage init -+>>>' + write(6,'(a15,a)') ' Current time: ',IO_timeStamp() #include "compilation_info.f90" - endif mainProcess !-------------------------------------------------------------------------------------------------- ! initialize solver specific parts of PETSc @@ -124,7 +121,8 @@ subroutine spectral_damage_init() CHKERRQ(ierr) call SNESSetDM(damage_snes,damage_grid,ierr); CHKERRQ(ierr) !< connect snes to da call DMCreateGlobalVector(damage_grid,solution,ierr); CHKERRQ(ierr) !< global solution vector (grid x 1, i.e. every def grad tensor) - call DMDASNESSetFunctionLocal(damage_grid,INSERT_VALUES,spectral_damage_formResidual,dummy,ierr) !< residual vector of same shape as solution vector + call DMDASNESSetFunctionLocal(damage_grid,INSERT_VALUES,spectral_damage_formResidual,& + PETSC_NULL_OBJECT,ierr) !< residual vector of same shape as solution vector CHKERRQ(ierr) call SNESSetFromOptions(damage_snes,ierr); CHKERRQ(ierr) !< pull it all together with additional cli arguments call SNESGetType(damage_snes,snes_type,ierr); CHKERRQ(ierr) @@ -171,7 +169,7 @@ end subroutine spectral_damage_init !-------------------------------------------------------------------------------------------------- !> @brief solution for the spectral damage scheme with internal iterations !-------------------------------------------------------------------------------------------------- -type(tSolutionState) function spectral_damage_solution(guess,timeinc,timeinc_old,loadCaseTime) +type(tSolutionState) function spectral_damage_solution(timeinc,timeinc_old,loadCaseTime) use numerics, only: & itmax, & err_damage_tolAbs, & @@ -190,7 +188,6 @@ type(tSolutionState) function spectral_damage_solution(guess,timeinc,timeinc_old timeinc, & !< increment in time for current solution timeinc_old, & !< increment in time of last increment loadCaseTime !< remaining time of current load case - logical, intent(in) :: guess integer(pInt) :: i, j, k, cell PetscInt ::position PetscReal :: minDamage, maxDamage, stagNorm, solnNorm @@ -283,10 +280,10 @@ subroutine spectral_damage_formResidual(in,x_scal,f_scal,dummy,ierr) DMDALocalInfo, dimension(DMDA_LOCAL_INFO_SIZE) :: & in PetscScalar, dimension( & - XG_RANGE,YG_RANGE,ZG_RANGE) :: & + XG_RANGE,YG_RANGE,ZG_RANGE), intent(in) :: & x_scal PetscScalar, dimension( & - X_RANGE,Y_RANGE,Z_RANGE) :: & + X_RANGE,Y_RANGE,Z_RANGE), intent(out) :: & f_scal PetscObject :: dummy PetscErrorCode :: ierr @@ -341,7 +338,7 @@ end subroutine spectral_damage_formResidual !-------------------------------------------------------------------------------------------------- !> @brief spectral damage forwarding routine !-------------------------------------------------------------------------------------------------- -subroutine spectral_damage_forward(guess,timeinc,timeinc_old,loadCaseTime) +subroutine spectral_damage_forward() use mesh, only: & grid, & grid3 @@ -354,11 +351,6 @@ subroutine spectral_damage_forward(guess,timeinc,timeinc_old,loadCaseTime) damage_nonlocal_getMobility implicit none - real(pReal), intent(in) :: & - timeinc_old, & - timeinc, & - loadCaseTime !< remaining time of current load case - logical, intent(in) :: guess integer(pInt) :: i, j, k, cell DM :: dm_local PetscScalar, dimension(:,:,:), pointer :: x_scal diff --git a/code/spectral_interface.f90 b/code/spectral_interface.f90 index 7050d7b45..8d9a41619 100644 --- a/code/spectral_interface.f90 +++ b/code/spectral_interface.f90 @@ -89,7 +89,7 @@ subroutine DAMASK_interface_init() call PetscInitialize(PETSC_NULL_CHARACTER,ierr) ! according to PETSc manual, that should be the first line in the code CHKERRQ(ierr) ! this is a macro definition, it is case sensitive call MPI_Comm_rank(PETSC_COMM_WORLD,worldrank,ierr);CHKERRQ(ierr) - call MPI_Comm_size(MPI_COMM_WORLD, worldsize, ierr);CHKERRQ(ierr) + call MPI_Comm_size(PETSC_COMM_WORLD,worldsize,ierr);CHKERRQ(ierr) mainProcess: if (worldrank == 0) then if (output_unit /= 6) then write(output_unit,'(a)') ' STDOUT != 6' diff --git a/code/spectral_mech_AL.f90 b/code/spectral_mech_AL.f90 index 0077ba2c3..951ab2521 100644 --- a/code/spectral_mech_AL.f90 +++ b/code/spectral_mech_AL.f90 @@ -49,7 +49,7 @@ module spectral_mech_AL F_av = 0.0_pReal, & !< average incompatible def grad field P_av = 0.0_pReal, & !< average 1st Piola--Kirchhoff stress P_avLastEval = 0.0_pReal !< average 1st Piola--Kirchhoff stress last call of CPFEM_general - character(len=1024), private :: incInfo !< time and increment information + character(len=1024), private :: incInfo !< time and increment information real(pReal), private, dimension(3,3,3,3) :: & C_volAvg = 0.0_pReal, & !< current volume average stiffness C_volAvgLastInc = 0.0_pReal, & !< previous volume average stiffness @@ -57,7 +57,7 @@ module spectral_mech_AL S = 0.0_pReal, & !< current compliance (filled up with zeros) C_scale = 0.0_pReal, & S_scale = 0.0_pReal - + real(pReal), private :: & err_BC, & !< deviation from stress BC err_curl, & !< RMS of curl of F @@ -116,7 +116,6 @@ subroutine AL_init temp33_Real = 0.0_pReal PetscErrorCode :: ierr - PetscObject :: dummy PetscScalar, pointer, dimension(:,:,:,:) :: xx_psc, F, F_lambda integer(pInt), dimension(:), allocatable :: localK integer(pInt) :: proc @@ -133,11 +132,9 @@ subroutine AL_init SNESSetConvergenceTest, & SNESSetFromOptions - mainProcess: if (worldrank == 0_pInt) then - write(6,'(/,a)') ' <<<+- DAMASK_spectral_solverAL init -+>>>' - write(6,'(a15,a)') ' Current time: ',IO_timeStamp() + write(6,'(/,a)') ' <<<+- DAMASK_spectral_solverAL init -+>>>' + write(6,'(a15,a)') ' Current time: ',IO_timeStamp() #include "compilation_info.f90" - endif mainProcess !-------------------------------------------------------------------------------------------------- ! allocate global fields @@ -165,9 +162,9 @@ subroutine AL_init CHKERRQ(ierr) call SNESSetDM(snes,da,ierr); CHKERRQ(ierr) call DMCreateGlobalVector(da,solution_vec,ierr); CHKERRQ(ierr) - call DMDASNESSetFunctionLocal(da,INSERT_VALUES,AL_formResidual,dummy,ierr) + call DMDASNESSetFunctionLocal(da,INSERT_VALUES,AL_formResidual,PETSC_NULL_OBJECT,ierr) CHKERRQ(ierr) - call SNESSetConvergenceTest(snes,AL_converged,dummy,PETSC_NULL_FUNCTION,ierr) + call SNESSetConvergenceTest(snes,AL_converged,PETSC_NULL_OBJECT,PETSC_NULL_FUNCTION,ierr) CHKERRQ(ierr) call SNESSetFromOptions(snes,ierr); CHKERRQ(ierr) @@ -245,7 +242,7 @@ end subroutine AL_init !> @brief solution for the AL scheme with internal iterations !-------------------------------------------------------------------------------------------------- type(tSolutionState) function & - AL_solution(incInfoIn,guess,timeinc,timeinc_old,loadCaseTime,P_BC,F_BC,rotation_BC) + AL_solution(incInfoIn,timeinc,timeinc_old,stress_BC,rotation_BC) use IO, only: & IO_error use numerics, only: & @@ -266,13 +263,9 @@ type(tSolutionState) function & ! input data for solution real(pReal), intent(in) :: & timeinc, & !< increment in time for current solution - timeinc_old, & !< increment in time of last increment - loadCaseTime !< remaining time of current load case - logical, intent(in) :: & - guess + timeinc_old !< increment in time of last increment type(tBoundaryCondition), intent(in) :: & - P_BC, & - F_BC + stress_BC character(len=*), intent(in) :: & incInfoIn real(pReal), dimension(3,3), intent(in) :: rotation_BC @@ -290,7 +283,7 @@ type(tSolutionState) function & !-------------------------------------------------------------------------------------------------- ! update stiffness (and gamma operator) - S = Utilities_maskedCompliance(rotation_BC,P_BC%maskLogical,C_volAvg) + S = Utilities_maskedCompliance(rotation_BC,stress_BC%maskLogical,C_volAvg) if (update_gamma) then call Utilities_updateGamma(C_minMaxAvg,restartWrite) C_scale = C_minMaxAvg @@ -299,11 +292,11 @@ type(tSolutionState) function & !-------------------------------------------------------------------------------------------------- ! set module wide availabe data - mask_stress = P_BC%maskFloat - params%P_BC = P_BC%values + mask_stress = stress_BC%maskFloat + params%stress_BC = stress_BC%values params%rotation_BC = rotation_BC - params%timeinc = timeinc - params%timeincOld = timeinc_old + params%timeinc = timeinc + params%timeincOld = timeinc_old !-------------------------------------------------------------------------------------------------- ! solve BVP @@ -331,8 +324,7 @@ subroutine AL_formResidual(in,x_scal,f_scal,dummy,ierr) itmax, & itmin, & polarAlpha, & - polarBeta, & - worldrank + polarBeta use mesh, only: & grid3, & grid @@ -369,10 +361,10 @@ subroutine AL_formResidual(in,x_scal,f_scal,dummy,ierr) DMDA_LOCAL_INFO_SIZE) :: & in PetscScalar, target, dimension(3,3,2, & - XG_RANGE,YG_RANGE,ZG_RANGE) :: & + XG_RANGE,YG_RANGE,ZG_RANGE), intent(in) :: & x_scal PetscScalar, target, dimension(3,3,2, & - X_RANGE,Y_RANGE,Z_RANGE) :: & + X_RANGE,Y_RANGE,Z_RANGE), intent(out) :: & f_scal PetscScalar, pointer, dimension(:,:,:,:,:) :: & F, & @@ -411,16 +403,14 @@ subroutine AL_formResidual(in,x_scal,f_scal,dummy,ierr) !-------------------------------------------------------------------------------------------------- ! report begin of new iteration totalIter = totalIter + 1_pInt - if (worldrank == 0_pInt) then - write(6,'(1x,a,3(a,'//IO_intOut(itmax)//'))') trim(incInfo), & - ' @ Iteration ', itmin, '≤',totalIter, '≤', itmax - if (iand(debug_level(debug_spectral),debug_spectralRotation) /= 0) & - write(6,'(/,a,/,3(3(f12.7,1x)/))',advance='no') ' deformation gradient aim (lab) =', & - math_transpose33(math_rotate_backward33(F_aim,params%rotation_BC)) - write(6,'(/,a,/,3(3(f12.7,1x)/))',advance='no') ' deformation gradient aim =', & - math_transpose33(F_aim) - flush(6) - endif + write(6,'(1x,a,3(a,'//IO_intOut(itmax)//'))') trim(incInfo), & + ' @ Iteration ', itmin, '≤',totalIter, '≤', itmax + if (iand(debug_level(debug_spectral),debug_spectralRotation) /= 0) & + write(6,'(/,a,/,3(3(f12.7,1x)/))',advance='no') ' deformation gradient aim (lab) =', & + math_transpose33(math_rotate_backward33(F_aim,params%rotation_BC)) + write(6,'(/,a,/,3(3(f12.7,1x)/))',advance='no') ' deformation gradient aim =', & + math_transpose33(F_aim) + flush(6) endif newIteration !-------------------------------------------------------------------------------------------------- @@ -495,8 +485,7 @@ subroutine AL_converged(snes_local,PETScIter,xnorm,snorm,fnorm,reason,dummy,ierr err_curl_tolRel, & err_curl_tolAbs, & err_stress_tolAbs, & - err_stress_tolRel, & - worldrank + err_stress_tolRel use math, only: & math_mul3333xx33 use FEsolving, only: & @@ -519,9 +508,9 @@ subroutine AL_converged(snes_local,PETScIter,xnorm,snorm,fnorm,reason,dummy,ierr !-------------------------------------------------------------------------------------------------- ! stress BC handling - F_aim = F_aim - math_mul3333xx33(S, ((P_av - params%P_BC))) ! S = 0.0 for no bc + F_aim = F_aim - math_mul3333xx33(S, ((P_av - params%stress_BC))) ! S = 0.0 for no bc err_BC = maxval(abs((-mask_stress+1.0_pReal)*math_mul3333xx33(C_scale,F_aim-F_av) + & - mask_stress *(P_av - params%P_BC))) ! mask = 0.0 for no bc + mask_stress *(P_av - params%stress_BC))) ! mask = 0.0 for no bc !-------------------------------------------------------------------------------------------------- ! error calculation @@ -543,24 +532,22 @@ subroutine AL_converged(snes_local,PETScIter,xnorm,snorm,fnorm,reason,dummy,ierr !-------------------------------------------------------------------------------------------------- ! report - if (worldrank == 0_pInt) then - write(6,'(1/,a)') ' ... reporting .............................................................' - write(6,'(/,a,f12.2,a,es8.2,a,es9.2,a)') ' error curl = ', & - err_curl/curlTol,' (',err_curl,' -, tol =',curlTol,')' - write(6,' (a,f12.2,a,es8.2,a,es9.2,a)') ' error divergence = ', & - err_div/divTol, ' (',err_div, ' / m, tol =',divTol,')' - write(6,' (a,f12.2,a,es8.2,a,es9.2,a)') ' error BC = ', & - err_BC/BC_tol, ' (',err_BC, ' Pa, tol =',BC_tol,')' - write(6,'(/,a)') ' ===========================================================================' - flush(6) - endif + write(6,'(1/,a)') ' ... reporting .............................................................' + write(6,'(/,a,f12.2,a,es8.2,a,es9.2,a)') ' error curl = ', & + err_curl/curlTol,' (',err_curl,' -, tol =',curlTol,')' + write(6,' (a,f12.2,a,es8.2,a,es9.2,a)') ' error divergence = ', & + err_div/divTol, ' (',err_div, ' / m, tol =',divTol,')' + write(6,' (a,f12.2,a,es8.2,a,es9.2,a)') ' error BC = ', & + err_BC/BC_tol, ' (',err_BC, ' Pa, tol =',BC_tol,')' + write(6,'(/,a)') ' ===========================================================================' + flush(6) end subroutine AL_converged !-------------------------------------------------------------------------------------------------- !> @brief forwarding routine !-------------------------------------------------------------------------------------------------- -subroutine AL_forward(guess,timeinc,timeinc_old,loadCaseTime,F_BC,P_BC,rotation_BC) +subroutine AL_forward(guess,timeinc,timeinc_old,loadCaseTime,deformation_BC,stress_BC,rotation_BC) use math, only: & math_mul33x33, & math_mul3333xx33, & @@ -588,8 +575,8 @@ subroutine AL_forward(guess,timeinc,timeinc_old,loadCaseTime,F_BC,P_BC,rotation_ timeinc, & loadCaseTime !< remaining time of current load case type(tBoundaryCondition), intent(in) :: & - P_BC, & - F_BC + stress_BC, & + deformation_BC real(pReal), dimension(3,3), intent(in) :: rotation_BC logical, intent(in) :: & guess @@ -605,21 +592,19 @@ subroutine AL_forward(guess,timeinc,timeinc_old,loadCaseTime,F_BC,P_BC,rotation_ F => xx_psc(0:8,:,:,:) F_lambda => xx_psc(9:17,:,:,:) if (restartWrite) then - if (worldrank == 0_pInt) then - write(6,'(/,a)') ' writing converged results for restart' - flush(6) - endif + write(6,'(/,a)') ' writing converged results for restart' + flush(6) write(rankStr,'(a1,i0)')'_',worldrank - call IO_write_jobRealFile(777,'F'//trim(rankStr),size(F)) ! writing deformation gradient field to file + call IO_write_jobRealFile(777,'F'//trim(rankStr),size(F)) ! writing deformation gradient field to file write (777,rec=1) F close (777) - call IO_write_jobRealFile(777,'F_lastInc'//trim(rankStr),size(F_lastInc)) ! writing F_lastInc field to file + call IO_write_jobRealFile(777,'F_lastInc'//trim(rankStr),size(F_lastInc)) ! writing F_lastInc field to file write (777,rec=1) F_lastInc close (777) - call IO_write_jobRealFile(777,'F_lambda'//trim(rankStr),size(F_lambda)) ! writing deformation gradient field to file + call IO_write_jobRealFile(777,'F_lambda'//trim(rankStr),size(F_lambda)) ! writing deformation gradient field to file write (777,rec=1) F_lambda close (777) - call IO_write_jobRealFile(777,'F_lambda_lastInc'//trim(rankStr),size(F_lambda_lastInc)) ! writing F_lastInc field to file + call IO_write_jobRealFile(777,'F_lambda_lastInc'//trim(rankStr),size(F_lambda_lastInc)) ! writing F_lastInc field to file write (777,rec=1) F_lambda_lastInc close (777) if (worldrank == 0_pInt) then @@ -653,14 +638,14 @@ subroutine AL_forward(guess,timeinc,timeinc_old,loadCaseTime,F_BC,P_BC,rotation_ C_volAvgLastInc = C_volAvg !-------------------------------------------------------------------------------------------------- ! calculate rate for aim - if (F_BC%myType=='l') then ! calculate f_aimDot from given L and current F - f_aimDot = F_BC%maskFloat * math_mul33x33(F_BC%values, F_aim) - elseif(F_BC%myType=='fdot') then ! f_aimDot is prescribed - f_aimDot = F_BC%maskFloat * F_BC%values - elseif(F_BC%myType=='f') then ! aim at end of load case is prescribed - f_aimDot = F_BC%maskFloat * (F_BC%values -F_aim)/loadCaseTime + if (deformation_BC%myType=='l') then ! calculate f_aimDot from given L and current F + f_aimDot = deformation_BC%maskFloat * math_mul33x33(deformation_BC%values, F_aim) + elseif(deformation_BC%myType=='fdot') then ! f_aimDot is prescribed + f_aimDot = deformation_BC%maskFloat * deformation_BC%values + elseif(deformation_BC%myType=='f') then ! aim at end of load case is prescribed + f_aimDot = deformation_BC%maskFloat * (deformation_BC%values -F_aim)/loadCaseTime endif - if (guess) f_aimDot = f_aimDot + P_BC%maskFloat * (F_aim - F_aim_lastInc)/timeinc_old + if (guess) f_aimDot = f_aimDot + stress_BC%maskFloat * (F_aim - F_aim_lastInc)/timeinc_old F_aim_lastInc = F_aim !-------------------------------------------------------------------------------------------------- diff --git a/code/spectral_mech_Basic.f90 b/code/spectral_mech_Basic.f90 index 445fbdf10..e20ed6761 100644 --- a/code/spectral_mech_Basic.f90 +++ b/code/spectral_mech_Basic.f90 @@ -105,7 +105,6 @@ subroutine basicPETSc_init temp33_Real = 0.0_pReal PetscErrorCode :: ierr - PetscObject :: dummy PetscScalar, pointer, dimension(:,:,:,:) :: F integer(pInt), dimension(:), allocatable :: localK @@ -123,11 +122,9 @@ subroutine basicPETSc_init SNESSetConvergenceTest, & SNESSetFromOptions - mainProcess: if (worldrank == 0_pInt) then - write(6,'(/,a)') ' <<<+- DAMASK_spectral_solverBasicPETSc init -+>>>' - write(6,'(a15,a)') ' Current time: ',IO_timeStamp() + write(6,'(/,a)') ' <<<+- DAMASK_spectral_solverBasicPETSc init -+>>>' + write(6,'(a15,a)') ' Current time: ',IO_timeStamp() #include "compilation_info.f90" - endif mainProcess !-------------------------------------------------------------------------------------------------- ! allocate global fields @@ -153,10 +150,10 @@ subroutine basicPETSc_init CHKERRQ(ierr) call SNESSetDM(snes,da,ierr); CHKERRQ(ierr) call DMCreateGlobalVector(da,solution_vec,ierr); CHKERRQ(ierr) ! global solution vector (grid x 9, i.e. every def grad tensor) - call DMDASNESSetFunctionLocal(da,INSERT_VALUES,BasicPETSC_formResidual,dummy,ierr) ! residual vector of same shape as solution vector + call DMDASNESSetFunctionLocal(da,INSERT_VALUES,BasicPETSC_formResidual,PETSC_NULL_OBJECT,ierr) ! residual vector of same shape as solution vector CHKERRQ(ierr) call SNESSetDM(snes,da,ierr); CHKERRQ(ierr) ! connect snes to da - call SNESSetConvergenceTest(snes,BasicPETSC_converged,dummy,PETSC_NULL_FUNCTION,ierr) ! specify custom convergence check function "_converged" + call SNESSetConvergenceTest(snes,BasicPETSC_converged,PETSC_NULL_OBJECT,PETSC_NULL_FUNCTION,ierr) ! specify custom convergence check function "_converged" CHKERRQ(ierr) call SNESSetFromOptions(snes,ierr); CHKERRQ(ierr) ! pull it all together with additional cli arguments @@ -165,7 +162,7 @@ subroutine basicPETSc_init call DMDAVecGetArrayF90(da,solution_vec,F,ierr); CHKERRQ(ierr) ! get the data out of PETSc to work with restart: if (restartInc > 1_pInt) then - if (iand(debug_level(debug_spectral),debug_spectralRestart)/= 0 .and. worldrank == 0_pInt) & + if (iand(debug_level(debug_spectral),debug_spectralRestart)/= 0) & write(6,'(/,a,'//IO_intOut(restartInc-1_pInt)//',a)') & 'reading values of increment ', restartInc - 1_pInt, ' from file' flush(6) @@ -220,7 +217,7 @@ end subroutine basicPETSc_init !> @brief solution for the Basic PETSC scheme with internal iterations !-------------------------------------------------------------------------------------------------- type(tSolutionState) function & - basicPETSc_solution(incInfoIn,guess,timeinc,timeinc_old,loadCaseTime,P_BC,F_BC,rotation_BC) + basicPETSc_solution(incInfoIn,timeinc,timeinc_old,stress_BC,rotation_BC) use IO, only: & IO_error use numerics, only: & @@ -239,13 +236,9 @@ type(tSolutionState) function & ! input data for solution real(pReal), intent(in) :: & timeinc, & !< increment in time for current solution - timeinc_old, & !< increment in time of last increment - loadCaseTime !< remaining time of current load case - logical, intent(in) :: & - guess + timeinc_old !< increment in time of last increment type(tBoundaryCondition), intent(in) :: & - P_BC, & - F_BC + stress_BC character(len=*), intent(in) :: & incInfoIn real(pReal), dimension(3,3), intent(in) :: rotation_BC @@ -263,17 +256,17 @@ type(tSolutionState) function & !-------------------------------------------------------------------------------------------------- ! update stiffness (and gamma operator) - S = Utilities_maskedCompliance(rotation_BC,P_BC%maskLogical,C_volAvg) + S = Utilities_maskedCompliance(rotation_BC,stress_BC%maskLogical,C_volAvg) if (update_gamma) call Utilities_updateGamma(C_minmaxAvg,restartWrite) !-------------------------------------------------------------------------------------------------- ! set module wide availabe data - mask_stress = P_BC%maskFloat - params%P_BC = P_BC%values + mask_stress = stress_BC%maskFloat + params%stress_BC = stress_BC%values params%rotation_BC = rotation_BC - params%timeinc = timeinc - params%timeincOld = timeinc_old + params%timeinc = timeinc + params%timeincOld = timeinc_old !-------------------------------------------------------------------------------------------------- ! solve BVP @@ -301,8 +294,6 @@ subroutine BasicPETSC_formResidual(in,x_scal,f_scal,dummy,ierr) use numerics, only: & itmax, & itmin - use numerics, only: & - worldrank use mesh, only: & grid, & grid3 @@ -331,10 +322,10 @@ subroutine BasicPETSC_formResidual(in,x_scal,f_scal,dummy,ierr) DMDALocalInfo, dimension(DMDA_LOCAL_INFO_SIZE) :: & in PetscScalar, dimension(3,3, & - XG_RANGE,YG_RANGE,ZG_RANGE) :: & + XG_RANGE,YG_RANGE,ZG_RANGE), intent(in) :: & x_scal PetscScalar, dimension(3,3, & - X_RANGE,Y_RANGE,Z_RANGE) :: & + X_RANGE,Y_RANGE,Z_RANGE), intent(out) :: & f_scal PetscInt :: & PETScIter, & @@ -354,16 +345,14 @@ subroutine BasicPETSC_formResidual(in,x_scal,f_scal,dummy,ierr) !-------------------------------------------------------------------------------------------------- ! report begin of new iteration totalIter = totalIter + 1_pInt - if (worldrank == 0_pInt) then - write(6,'(1x,a,3(a,'//IO_intOut(itmax)//'))') trim(incInfo), & - ' @ Iteration ', itmin, '≤',totalIter, '≤', itmax - if (iand(debug_level(debug_spectral),debug_spectralRotation) /= 0) & - write(6,'(/,a,/,3(3(f12.7,1x)/))',advance='no') ' deformation gradient aim (lab) =', & - math_transpose33(math_rotate_backward33(F_aim,params%rotation_BC)) - write(6,'(/,a,/,3(3(f12.7,1x)/))',advance='no') ' deformation gradient aim =', & - math_transpose33(F_aim) - flush(6) - endif + write(6,'(1x,a,3(a,'//IO_intOut(itmax)//'))') trim(incInfo), & + ' @ Iteration ', itmin, '≤',totalIter, '≤', itmax + if (iand(debug_level(debug_spectral),debug_spectralRotation) /= 0) & + write(6,'(/,a,/,3(3(f12.7,1x)/))',advance='no') ' deformation gradient aim (lab) =', & + math_transpose33(math_rotate_backward33(F_aim,params%rotation_BC)) + write(6,'(/,a,/,3(3(f12.7,1x)/))',advance='no') ' deformation gradient aim =', & + math_transpose33(F_aim) + flush(6) endif newIteration !-------------------------------------------------------------------------------------------------- @@ -376,8 +365,8 @@ subroutine BasicPETSC_formResidual(in,x_scal,f_scal,dummy,ierr) !-------------------------------------------------------------------------------------------------- ! stress BC handling F_aim_lastIter = F_aim - F_aim = F_aim - math_mul3333xx33(S, ((P_av - params%P_BC))) ! S = 0.0 for no bc - err_stress = maxval(abs(mask_stress * (P_av - params%P_BC))) ! mask = 0.0 for no bc + F_aim = F_aim - math_mul3333xx33(S, ((P_av - params%stress_BC))) ! S = 0.0 for no bc + err_stress = maxval(abs(mask_stress * (P_av - params%stress_BC))) ! mask = 0.0 for no bc !-------------------------------------------------------------------------------------------------- ! updated deformation gradient using fix point algorithm of basic scheme @@ -405,8 +394,7 @@ subroutine BasicPETSc_converged(snes_local,PETScIter,xnorm,snorm,fnorm,reason,du err_div_tolRel, & err_div_tolAbs, & err_stress_tolRel, & - err_stress_tolAbs, & - worldrank + err_stress_tolAbs use FEsolving, only: & terminallyIll @@ -440,22 +428,20 @@ subroutine BasicPETSc_converged(snes_local,PETScIter,xnorm,snorm,fnorm,reason,du !-------------------------------------------------------------------------------------------------- ! report - if (worldrank == 0_pInt) then - write(6,'(1/,a)') ' ... reporting .............................................................' - write(6,'(1/,a,f12.2,a,es8.2,a,es9.2,a)') ' error divergence = ', & - err_div/divTol, ' (',err_div,' / m, tol =',divTol,')' - write(6,'(a,f12.2,a,es8.2,a,es9.2,a)') ' error stress BC = ', & - err_stress/stressTol, ' (',err_stress, ' Pa, tol =',stressTol,')' - write(6,'(/,a)') ' ===========================================================================' - flush(6) - endif + write(6,'(1/,a)') ' ... reporting .............................................................' + write(6,'(1/,a,f12.2,a,es8.2,a,es9.2,a)') ' error divergence = ', & + err_div/divTol, ' (',err_div,' / m, tol =',divTol,')' + write(6,'(a,f12.2,a,es8.2,a,es9.2,a)') ' error stress BC = ', & + err_stress/stressTol, ' (',err_stress, ' Pa, tol =',stressTol,')' + write(6,'(/,a)') ' ===========================================================================' +flush(6) end subroutine BasicPETSc_converged !-------------------------------------------------------------------------------------------------- !> @brief forwarding routine !-------------------------------------------------------------------------------------------------- -subroutine BasicPETSc_forward(guess,timeinc,timeinc_old,loadCaseTime,F_BC,P_BC,rotation_BC) +subroutine BasicPETSc_forward(guess,timeinc,timeinc_old,loadCaseTime,deformation_BC,stress_BC,rotation_BC) use math, only: & math_mul33x33 ,& math_rotate_backward33 @@ -481,8 +467,8 @@ subroutine BasicPETSc_forward(guess,timeinc,timeinc_old,loadCaseTime,F_BC,P_BC,r timeinc, & loadCaseTime !< remaining time of current load case type(tBoundaryCondition), intent(in) :: & - P_BC, & - F_BC + stress_BC, & + deformation_BC real(pReal), dimension(3,3), intent(in) :: rotation_BC logical, intent(in) :: & guess @@ -495,10 +481,8 @@ subroutine BasicPETSc_forward(guess,timeinc,timeinc_old,loadCaseTime,F_BC,P_BC,r !-------------------------------------------------------------------------------------------------- ! restart information for spectral solver if (restartWrite) then - if (worldrank == 0_pInt) then - write(6,'(/,a)') ' writing converged results for restart' - flush(6) - endif + write(6,'(/,a)') ' writing converged results for restart' + flush(6) write(rankStr,'(a1,i0)')'_',worldrank call IO_write_jobRealFile(777,'F'//trim(rankStr),size(F)) ! writing deformation gradient field to file write (777,rec=1) F @@ -530,14 +514,14 @@ subroutine BasicPETSc_forward(guess,timeinc,timeinc_old,loadCaseTime,F_BC,P_BC,r C_volAvgLastInc = C_volAvg !-------------------------------------------------------------------------------------------------- ! calculate rate for aim - if (F_BC%myType=='l') then ! calculate f_aimDot from given L and current F - f_aimDot = F_BC%maskFloat * math_mul33x33(F_BC%values, F_aim) - elseif(F_BC%myType=='fdot') then ! f_aimDot is prescribed - f_aimDot = F_BC%maskFloat * F_BC%values - elseif(F_BC%myType=='f') then ! aim at end of load case is prescribed - f_aimDot = F_BC%maskFloat * (F_BC%values -F_aim)/loadCaseTime + if (deformation_BC%myType=='l') then ! calculate f_aimDot from given L and current F + f_aimDot = deformation_BC%maskFloat * math_mul33x33(deformation_BC%values, F_aim) + elseif(deformation_BC%myType=='fdot') then ! f_aimDot is prescribed + f_aimDot = deformation_BC%maskFloat * deformation_BC%values + elseif(deformation_BC%myType=='f') then ! aim at end of load case is prescribed + f_aimDot = deformation_BC%maskFloat * (deformation_BC%values -F_aim)/loadCaseTime endif - if (guess) f_aimDot = f_aimDot + P_BC%maskFloat * (F_aim - F_aim_lastInc)/timeinc_old + if (guess) f_aimDot = f_aimDot + stress_BC%maskFloat * (F_aim - F_aim_lastInc)/timeinc_old F_aim_lastInc = F_aim !-------------------------------------------------------------------------------------------------- @@ -552,8 +536,8 @@ subroutine BasicPETSc_forward(guess,timeinc,timeinc_old,loadCaseTime,F_BC,P_BC,r !-------------------------------------------------------------------------------------------------- ! update local deformation gradient - F = reshape(Utilities_forwardField(timeinc,F_lastInc,Fdot, & ! ensure that it matches rotated F_aim - math_rotate_backward33(F_aim,rotation_BC)),[9,grid(1),grid(2),grid3]) + F = reshape(Utilities_forwardField(timeinc,F_lastInc,Fdot, & ! ensure that it matches rotated F_aim + math_rotate_backward33(F_aim,rotation_BC)),[9,grid(1),grid(2),grid3]) call DMDAVecRestoreArrayF90(da,solution_vec,F,ierr); CHKERRQ(ierr) end subroutine BasicPETSc_forward diff --git a/code/spectral_mech_Polarisation.f90 b/code/spectral_mech_Polarisation.f90 index fc9fae7d7..ed44793bb 100644 --- a/code/spectral_mech_Polarisation.f90 +++ b/code/spectral_mech_Polarisation.f90 @@ -116,7 +116,6 @@ subroutine Polarisation_init temp33_Real = 0.0_pReal PetscErrorCode :: ierr - PetscObject :: dummy PetscScalar, pointer, dimension(:,:,:,:) :: xx_psc, F, F_tau integer(pInt), dimension(:), allocatable :: localK integer(pInt) :: proc @@ -133,11 +132,9 @@ subroutine Polarisation_init SNESSetConvergenceTest, & SNESSetFromOptions - mainProcess: if (worldrank == 0_pInt) then - write(6,'(/,a)') ' <<<+- DAMASK_spectral_solverPolarisation init -+>>>' - write(6,'(a15,a)') ' Current time: ',IO_timeStamp() + write(6,'(/,a)') ' <<<+- DAMASK_spectral_solverPolarisation init -+>>>' + write(6,'(a15,a)') ' Current time: ',IO_timeStamp() #include "compilation_info.f90" - endif mainProcess !-------------------------------------------------------------------------------------------------- ! allocate global fields @@ -165,9 +162,9 @@ subroutine Polarisation_init CHKERRQ(ierr) call SNESSetDM(snes,da,ierr); CHKERRQ(ierr) call DMCreateGlobalVector(da,solution_vec,ierr); CHKERRQ(ierr) - call DMDASNESSetFunctionLocal(da,INSERT_VALUES,Polarisation_formResidual,dummy,ierr) + call DMDASNESSetFunctionLocal(da,INSERT_VALUES,Polarisation_formResidual,PETSC_NULL_OBJECT,ierr) CHKERRQ(ierr) - call SNESSetConvergenceTest(snes,Polarisation_converged,dummy,PETSC_NULL_FUNCTION,ierr) + call SNESSetConvergenceTest(snes,Polarisation_converged,PETSC_NULL_OBJECT,PETSC_NULL_FUNCTION,ierr) CHKERRQ(ierr) call SNESSetFromOptions(snes,ierr); CHKERRQ(ierr) @@ -177,7 +174,7 @@ subroutine Polarisation_init F => xx_psc(0:8,:,:,:) F_tau => xx_psc(9:17,:,:,:) restart: if (restartInc > 1_pInt) then - if (iand(debug_level(debug_spectral),debug_spectralRestart)/= 0 .and. worldrank == 0_pInt) & + if (iand(debug_level(debug_spectral),debug_spectralRestart)/= 0) & write(6,'(/,a,'//IO_intOut(restartInc-1_pInt)//',a)') & 'reading values of increment ', restartInc - 1_pInt, ' from file' flush(6) @@ -245,7 +242,7 @@ end subroutine Polarisation_init !> @brief solution for the Polarisation scheme with internal iterations !-------------------------------------------------------------------------------------------------- type(tSolutionState) function & - Polarisation_solution(incInfoIn,guess,timeinc,timeinc_old,loadCaseTime,P_BC,F_BC,rotation_BC) + Polarisation_solution(incInfoIn,timeinc,timeinc_old,stress_BC,rotation_BC) use IO, only: & IO_error use numerics, only: & @@ -266,13 +263,9 @@ type(tSolutionState) function & ! input data for solution real(pReal), intent(in) :: & timeinc, & !< increment in time for current solution - timeinc_old, & !< increment in time of last increment - loadCaseTime !< remaining time of current load case - logical, intent(in) :: & - guess + timeinc_old !< increment in time of last increment type(tBoundaryCondition), intent(in) :: & - P_BC, & - F_BC + stress_BC character(len=*), intent(in) :: & incInfoIn real(pReal), dimension(3,3), intent(in) :: rotation_BC @@ -290,7 +283,7 @@ type(tSolutionState) function & !-------------------------------------------------------------------------------------------------- ! update stiffness (and gamma operator) - S = Utilities_maskedCompliance(rotation_BC,P_BC%maskLogical,C_volAvg) + S = Utilities_maskedCompliance(rotation_BC,stress_BC%maskLogical,C_volAvg) if (update_gamma) then call Utilities_updateGamma(C_minMaxAvg,restartWrite) C_scale = C_minMaxAvg @@ -299,11 +292,11 @@ type(tSolutionState) function & !-------------------------------------------------------------------------------------------------- ! set module wide availabe data - mask_stress = P_BC%maskFloat - params%P_BC = P_BC%values + mask_stress = stress_BC%maskFloat + params%stress_BC = stress_BC%values params%rotation_BC = rotation_BC - params%timeinc = timeinc - params%timeincOld = timeinc_old + params%timeinc = timeinc + params%timeincOld = timeinc_old !-------------------------------------------------------------------------------------------------- ! solve BVP @@ -331,8 +324,7 @@ subroutine Polarisation_formResidual(in,x_scal,f_scal,dummy,ierr) itmax, & itmin, & polarAlpha, & - polarBeta, & - worldrank + polarBeta use mesh, only: & grid3, & grid @@ -369,10 +361,10 @@ subroutine Polarisation_formResidual(in,x_scal,f_scal,dummy,ierr) DMDA_LOCAL_INFO_SIZE) :: & in PetscScalar, target, dimension(3,3,2, & - XG_RANGE,YG_RANGE,ZG_RANGE) :: & + XG_RANGE,YG_RANGE,ZG_RANGE), intent(in) :: & x_scal PetscScalar, target, dimension(3,3,2, & - X_RANGE,Y_RANGE,Z_RANGE) :: & + X_RANGE,Y_RANGE,Z_RANGE), intent(out) :: & f_scal PetscScalar, pointer, dimension(:,:,:,:,:) :: & F, & @@ -411,16 +403,14 @@ subroutine Polarisation_formResidual(in,x_scal,f_scal,dummy,ierr) !-------------------------------------------------------------------------------------------------- ! report begin of new iteration totalIter = totalIter + 1_pInt - if (worldrank == 0_pInt) then - write(6,'(1x,a,3(a,'//IO_intOut(itmax)//'))') trim(incInfo), & - ' @ Iteration ', itmin, '≤',totalIter, '≤', itmax - if (iand(debug_level(debug_spectral),debug_spectralRotation) /= 0) & - write(6,'(/,a,/,3(3(f12.7,1x)/))',advance='no') ' deformation gradient aim (lab) =', & - math_transpose33(math_rotate_backward33(F_aim,params%rotation_BC)) - write(6,'(/,a,/,3(3(f12.7,1x)/))',advance='no') ' deformation gradient aim =', & - math_transpose33(F_aim) - flush(6) - endif + write(6,'(1x,a,3(a,'//IO_intOut(itmax)//'))') trim(incInfo), & + ' @ Iteration ', itmin, '≤',totalIter, '≤', itmax + if (iand(debug_level(debug_spectral),debug_spectralRotation) /= 0) & + write(6,'(/,a,/,3(3(f12.7,1x)/))',advance='no') ' deformation gradient aim (lab) =', & + math_transpose33(math_rotate_backward33(F_aim,params%rotation_BC)) + write(6,'(/,a,/,3(3(f12.7,1x)/))',advance='no') ' deformation gradient aim =', & + math_transpose33(F_aim) + flush(6) endif newIteration !-------------------------------------------------------------------------------------------------- @@ -494,8 +484,7 @@ subroutine Polarisation_converged(snes_local,PETScIter,xnorm,snorm,fnorm,reason, err_curl_tolRel, & err_curl_tolAbs, & err_stress_tolAbs, & - err_stress_tolRel, & - worldrank + err_stress_tolRel use math, only: & math_mul3333xx33 use FEsolving, only: & @@ -518,9 +507,9 @@ subroutine Polarisation_converged(snes_local,PETScIter,xnorm,snorm,fnorm,reason, !-------------------------------------------------------------------------------------------------- ! stress BC handling - F_aim = F_aim - math_mul3333xx33(S, ((P_av - params%P_BC))) ! S = 0.0 for no bc + F_aim = F_aim - math_mul3333xx33(S, ((P_av - params%stress_BC))) ! S = 0.0 for no bc err_BC = maxval(abs((-mask_stress+1.0_pReal)*math_mul3333xx33(C_scale,F_aim-F_av) + & - mask_stress *(P_av - params%P_BC))) ! mask = 0.0 for no bc + mask_stress *(P_av - params%stress_BC))) ! mask = 0.0 for no bc !-------------------------------------------------------------------------------------------------- ! error calculation @@ -542,31 +531,29 @@ subroutine Polarisation_converged(snes_local,PETScIter,xnorm,snorm,fnorm,reason, !-------------------------------------------------------------------------------------------------- ! report -if (worldrank == 0_pInt) then - write(6,'(1/,a)') ' ... reporting .............................................................' - write(6,'(/,a,f12.2,a,es8.2,a,es9.2,a)') ' error curl = ', & - err_curl/curlTol,' (',err_curl,' -, tol =',curlTol,')' - write(6,' (a,f12.2,a,es8.2,a,es9.2,a)') ' error divergence = ', & - err_div/divTol, ' (',err_div, ' / m, tol =',divTol,')' - write(6,' (a,f12.2,a,es8.2,a,es9.2,a)') ' error BC = ', & - err_BC/BC_tol, ' (',err_BC, ' Pa, tol =',BC_tol,')' - write(6,'(/,a)') ' ===========================================================================' - flush(6) -endif + write(6,'(1/,a)') ' ... reporting .............................................................' + write(6,'(/,a,f12.2,a,es8.2,a,es9.2,a)') ' error curl = ', & + err_curl/curlTol,' (',err_curl,' -, tol =',curlTol,')' + write(6,' (a,f12.2,a,es8.2,a,es9.2,a)') ' error divergence = ', & + err_div/divTol, ' (',err_div, ' / m, tol =',divTol,')' + write(6,' (a,f12.2,a,es8.2,a,es9.2,a)') ' error BC = ', & + err_BC/BC_tol, ' (',err_BC, ' Pa, tol =',BC_tol,')' + write(6,'(/,a)') ' ===========================================================================' + flush(6) end subroutine Polarisation_converged !-------------------------------------------------------------------------------------------------- !> @brief forwarding routine !-------------------------------------------------------------------------------------------------- -subroutine Polarisation_forward(guess,timeinc,timeinc_old,loadCaseTime,F_BC,P_BC,rotation_BC) +subroutine Polarisation_forward(guess,timeinc,timeinc_old,loadCaseTime,deformation_BC,stress_BC,rotation_BC) use math, only: & math_mul33x33, & math_mul3333xx33, & math_transpose33, & math_rotate_backward33 use numerics, only: & - worldrank + worldrank use mesh, only: & grid3, & grid @@ -587,8 +574,8 @@ subroutine Polarisation_forward(guess,timeinc,timeinc_old,loadCaseTime,F_BC,P_BC timeinc, & loadCaseTime !< remaining time of current load case type(tBoundaryCondition), intent(in) :: & - P_BC, & - F_BC + stress_BC, & + deformation_BC real(pReal), dimension(3,3), intent(in) :: rotation_BC logical, intent(in) :: & guess @@ -604,19 +591,19 @@ subroutine Polarisation_forward(guess,timeinc,timeinc_old,loadCaseTime,F_BC,P_BC F => xx_psc(0:8,:,:,:) F_tau => xx_psc(9:17,:,:,:) if (restartWrite) then - if (worldrank == 0_pInt) write(6,'(/,a)') ' writing converged results for restart' + write(6,'(/,a)') ' writing converged results for restart' flush(6) write(rankStr,'(a1,i0)')'_',worldrank - call IO_write_jobRealFile(777,'F'//trim(rankStr),size(F)) ! writing deformation gradient field to file + call IO_write_jobRealFile(777,'F'//trim(rankStr),size(F)) ! writing deformation gradient field to file write (777,rec=1) F close (777) - call IO_write_jobRealFile(777,'F_lastInc'//trim(rankStr),size(F_lastInc)) ! writing F_lastInc field to file + call IO_write_jobRealFile(777,'F_lastInc'//trim(rankStr),size(F_lastInc)) ! writing F_lastInc field to file write (777,rec=1) F_lastInc close (777) - call IO_write_jobRealFile(777,'F_tau'//trim(rankStr),size(F_tau)) ! writing deformation gradient field to file + call IO_write_jobRealFile(777,'F_tau'//trim(rankStr),size(F_tau)) ! writing deformation gradient field to file write (777,rec=1) F_tau close (777) - call IO_write_jobRealFile(777,'F_tau_lastInc'//trim(rankStr),size(F_tau_lastInc)) ! writing F_lastInc field to file + call IO_write_jobRealFile(777,'F_tau_lastInc'//trim(rankStr),size(F_tau_lastInc)) ! writing F_lastInc field to file write (777,rec=1) F_tau_lastInc close (777) if (worldrank == 0_pInt) then @@ -650,14 +637,14 @@ subroutine Polarisation_forward(guess,timeinc,timeinc_old,loadCaseTime,F_BC,P_BC C_volAvgLastInc = C_volAvg !-------------------------------------------------------------------------------------------------- ! calculate rate for aim - if (F_BC%myType=='l') then ! calculate f_aimDot from given L and current F - f_aimDot = F_BC%maskFloat * math_mul33x33(F_BC%values, F_aim) - elseif(F_BC%myType=='fdot') then ! f_aimDot is prescribed - f_aimDot = F_BC%maskFloat * F_BC%values - elseif(F_BC%myType=='f') then ! aim at end of load case is prescribed - f_aimDot = F_BC%maskFloat * (F_BC%values -F_aim)/loadCaseTime + if (deformation_BC%myType=='l') then ! calculate f_aimDot from given L and current F + f_aimDot = deformation_BC%maskFloat * math_mul33x33(deformation_BC%values, F_aim) + elseif(deformation_BC%myType=='fdot') then ! f_aimDot is prescribed + f_aimDot = deformation_BC%maskFloat * deformation_BC%values + elseif(deformation_BC%myType=='f') then ! aim at end of load case is prescribed + f_aimDot = deformation_BC%maskFloat * (deformation_BC%values -F_aim)/loadCaseTime endif - if (guess) f_aimDot = f_aimDot + P_BC%maskFloat * (F_aim - F_aim_lastInc)/timeinc_old + if (guess) f_aimDot = f_aimDot + stress_BC%maskFloat * (F_aim - F_aim_lastInc)/timeinc_old F_aim_lastInc = F_aim !-------------------------------------------------------------------------------------------------- diff --git a/code/spectral_thermal.f90 b/code/spectral_thermal.f90 index bb5747574..490325ab7 100644 --- a/code/spectral_thermal.f90 +++ b/code/spectral_thermal.f90 @@ -86,7 +86,6 @@ subroutine spectral_thermal_init DM :: thermal_grid PetscScalar, dimension(:,:,:), pointer :: x_scal PetscErrorCode :: ierr - PetscObject :: dummy external :: & SNESCreate, & @@ -123,7 +122,8 @@ subroutine spectral_thermal_init CHKERRQ(ierr) call SNESSetDM(thermal_snes,thermal_grid,ierr); CHKERRQ(ierr) ! connect snes to da call DMCreateGlobalVector(thermal_grid,solution ,ierr); CHKERRQ(ierr) ! global solution vector (grid x 1, i.e. every def grad tensor) - call DMDASNESSetFunctionLocal(thermal_grid,INSERT_VALUES,spectral_thermal_formResidual,dummy,ierr) ! residual vector of same shape as solution vector + call DMDASNESSetFunctionLocal(thermal_grid,INSERT_VALUES,spectral_thermal_formResidual,& + PETSC_NULL_OBJECT,ierr) ! residual vector of same shape as solution vector CHKERRQ(ierr) call SNESSetFromOptions(thermal_snes,ierr); CHKERRQ(ierr) ! pull it all together with additional cli arguments @@ -170,7 +170,7 @@ end subroutine spectral_thermal_init !-------------------------------------------------------------------------------------------------- !> @brief solution for the spectral thermal scheme with internal iterations !-------------------------------------------------------------------------------------------------- -type(tSolutionState) function spectral_thermal_solution(guess,timeinc,timeinc_old,loadCaseTime) +type(tSolutionState) function spectral_thermal_solution(timeinc,timeinc_old,loadCaseTime) use numerics, only: & itmax, & err_thermal_tolAbs, & @@ -189,7 +189,6 @@ type(tSolutionState) function spectral_thermal_solution(guess,timeinc,timeinc_ol timeinc, & !< increment in time for current solution timeinc_old, & !< increment in time of last increment loadCaseTime !< remaining time of current load case - logical, intent(in) :: guess integer(pInt) :: i, j, k, cell PetscInt :: position PetscReal :: minTemperature, maxTemperature, stagNorm, solnNorm @@ -283,10 +282,10 @@ subroutine spectral_thermal_formResidual(in,x_scal,f_scal,dummy,ierr) DMDALocalInfo, dimension(DMDA_LOCAL_INFO_SIZE) :: & in PetscScalar, dimension( & - XG_RANGE,YG_RANGE,ZG_RANGE) :: & + XG_RANGE,YG_RANGE,ZG_RANGE), intent(in) :: & x_scal PetscScalar, dimension( & - X_RANGE,Y_RANGE,Z_RANGE) :: & + X_RANGE,Y_RANGE,Z_RANGE), intent(out) :: & f_scal PetscObject :: dummy PetscErrorCode :: ierr @@ -337,7 +336,7 @@ end subroutine spectral_thermal_formResidual !-------------------------------------------------------------------------------------------------- !> @brief forwarding routine !-------------------------------------------------------------------------------------------------- -subroutine spectral_thermal_forward(guess,timeinc,timeinc_old,loadCaseTime) +subroutine spectral_thermal_forward() use mesh, only: & grid, & grid3 @@ -351,11 +350,6 @@ subroutine spectral_thermal_forward(guess,timeinc,timeinc_old,loadCaseTime) thermal_conduction_getSpecificHeat implicit none - real(pReal), intent(in) :: & - timeinc_old, & - timeinc, & - loadCaseTime !< remaining time of current load case - logical, intent(in) :: guess integer(pInt) :: i, j, k, cell DM :: dm_local PetscScalar, dimension(:,:,:), pointer :: x_scal diff --git a/code/spectral_utilities.f90 b/code/spectral_utilities.f90 index 1a86b7648..34eb0eab0 100644 --- a/code/spectral_utilities.f90 +++ b/code/spectral_utilities.f90 @@ -1,3 +1,4 @@ +!-------------------------------------------------------------------------------------------------- !> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH !> @author Philip Eisenlohr, Max-Planck-Institut für Eisenforschung GmbH !> @brief Utilities used by the different spectral solver variants @@ -84,7 +85,7 @@ module spectral_utilities type, public :: tLoadCase real(pReal), dimension (3,3) :: rotation = math_I3 !< rotation of BC - type(tBoundaryCondition) :: P, & !< stress BC + type(tBoundaryCondition) :: stress, & !< stress BC deformation !< deformation BC (Fdot or L) real(pReal) :: time = 0.0_pReal !< length of increment integer(pInt) :: incs = 0_pInt, & !< number of increments @@ -96,7 +97,7 @@ module spectral_utilities end type tLoadCase type, public :: tSolutionParams !< @todo use here the type definition for a full loadcase including mask - real(pReal), dimension(3,3) :: P_BC, rotation_BC + real(pReal), dimension(3,3) :: stress_BC, rotation_BC real(pReal) :: timeinc real(pReal) :: timeincOld real(pReal) :: density @@ -227,9 +228,12 @@ subroutine utilities_init() ' add more using the PETSc_Options keyword in numerics.config '; flush(6) call PetscOptionsClear(ierr); CHKERRQ(ierr) - if(debugPETSc) call PetscOptionsInsertString(trim(PETSCDEBUG),ierr); CHKERRQ(ierr) - call PetscOptionsInsertString(trim(petsc_defaultOptions),ierr); CHKERRQ(ierr) - call PetscOptionsInsertString(trim(petsc_options),ierr); CHKERRQ(ierr) + if(debugPETSc) call PetscOptionsInsertString(trim(PETSCDEBUG),ierr) + CHKERRQ(ierr) + call PetscOptionsInsertString(trim(petsc_defaultOptions),ierr) + CHKERRQ(ierr) + call PetscOptionsInsertString(trim(petsc_options),ierr) + CHKERRQ(ierr) grid1Red = grid(1)/2_pInt + 1_pInt wgt = 1.0/real(product(grid),pReal) @@ -238,11 +242,11 @@ subroutine utilities_init() write(6,'(a,3(es12.5))') ' size x y z: ', geomSize select case (spectral_derivative) - case ('continuous') ! default, no weighting + case ('continuous') spectral_derivative_ID = DERIVATIVE_CONTINUOUS_ID - case ('central_difference') ! cosine curve with 1 for avg and zero for highest freq + case ('central_difference') spectral_derivative_ID = DERIVATIVE_CENTRAL_DIFF_ID - case ('fwbw_difference') ! gradient, might need grid scaling as for cosine filter + case ('fwbw_difference') spectral_derivative_ID = DERIVATIVE_FWBW_DIFF_ID case default call IO_error(892_pInt,ext_msg=trim(spectral_derivative)) @@ -271,9 +275,9 @@ subroutine utilities_init() ! MPI allocation gridFFTW = int(grid,C_INTPTR_T) alloc_local = fftw_mpi_local_size_3d(gridFFTW(3), gridFFTW(2), gridFFTW(1)/2 +1, & - MPI_COMM_WORLD, local_K, local_K_offset) - allocate (xi1st (3,grid1Red,grid(2),grid3),source = cmplx(0.0_pReal,0.0_pReal,pReal)) ! frequencies, only half the size for first dimension - allocate (xi2nd (3,grid1Red,grid(2),grid3),source = cmplx(0.0_pReal,0.0_pReal,pReal)) ! frequencies, only half the size for first dimension + PETSC_COMM_WORLD, local_K, local_K_offset) + allocate (xi1st (3,grid1Red,grid(2),grid3),source = cmplx(0.0_pReal,0.0_pReal,pReal)) ! frequencies for first derivatives, only half the size for first dimension + allocate (xi2nd (3,grid1Red,grid(2),grid3),source = cmplx(0.0_pReal,0.0_pReal,pReal)) ! frequencies for second derivatives, only half the size for first dimension tensorField = fftw_alloc_complex(tensorSize*alloc_local) call c_f_pointer(tensorField, tensorField_real, [3_C_INTPTR_T,3_C_INTPTR_T, & @@ -298,12 +302,12 @@ subroutine utilities_init() planTensorForth = fftw_mpi_plan_many_dft_r2c(3, [gridFFTW(3),gridFFTW(2),gridFFTW(1)], & ! dimension, logical length in each dimension in reversed order tensorSize, FFTW_MPI_DEFAULT_BLOCK, FFTW_MPI_DEFAULT_BLOCK, &! no. of transforms, default iblock and oblock tensorField_real, tensorField_fourier, & ! input data, output data - MPI_COMM_WORLD, fftw_planner_flag) ! use all processors, planer precision + PETSC_COMM_WORLD, fftw_planner_flag) ! use all processors, planer precision if (.not. C_ASSOCIATED(planTensorForth)) call IO_error(810, ext_msg='planTensorForth') planTensorBack = fftw_mpi_plan_many_dft_c2r(3, [gridFFTW(3),gridFFTW(2),gridFFTW(1)], & ! dimension, logical length in each dimension in reversed order tensorSize, FFTW_MPI_DEFAULT_BLOCK, FFTW_MPI_DEFAULT_BLOCK, &! no. of transforms, default iblock and oblock tensorField_fourier,tensorField_real, & ! input data, output data - MPI_COMM_WORLD, fftw_planner_flag) ! all processors, planer precision + PETSC_COMM_WORLD, fftw_planner_flag) ! all processors, planer precision if (.not. C_ASSOCIATED(planTensorBack)) call IO_error(810, ext_msg='planTensorBack') !-------------------------------------------------------------------------------------------------- @@ -311,12 +315,12 @@ subroutine utilities_init() planVectorForth = fftw_mpi_plan_many_dft_r2c(3, [gridFFTW(3),gridFFTW(2),gridFFTW(1)], & ! dimension, logical length in each dimension in reversed order vecSize, FFTW_MPI_DEFAULT_BLOCK, FFTW_MPI_DEFAULT_BLOCK, &! no. of transforms, default iblock and oblock vectorField_real, vectorField_fourier, & ! input data, output data - MPI_COMM_WORLD, fftw_planner_flag) ! use all processors, planer precision + PETSC_COMM_WORLD, fftw_planner_flag) ! use all processors, planer precision if (.not. C_ASSOCIATED(planVectorForth)) call IO_error(810, ext_msg='planVectorForth') planVectorBack = fftw_mpi_plan_many_dft_c2r(3, [gridFFTW(3),gridFFTW(2),gridFFTW(1)], & ! dimension, logical length in each dimension in reversed order vecSize, FFTW_MPI_DEFAULT_BLOCK, FFTW_MPI_DEFAULT_BLOCK, & ! no. of transforms, default iblock and oblock vectorField_fourier,vectorField_real, & ! input data, output data - MPI_COMM_WORLD, fftw_planner_flag) ! all processors, planer precision + PETSC_COMM_WORLD, fftw_planner_flag) ! all processors, planer precision if (.not. C_ASSOCIATED(planVectorBack)) call IO_error(810, ext_msg='planVectorBack') !-------------------------------------------------------------------------------------------------- @@ -324,12 +328,12 @@ subroutine utilities_init() planScalarForth = fftw_mpi_plan_many_dft_r2c(3, [gridFFTW(3),gridFFTW(2),gridFFTW(1)], & ! dimension, logical length in each dimension in reversed order scalarSize, FFTW_MPI_DEFAULT_BLOCK, FFTW_MPI_DEFAULT_BLOCK, & ! no. of transforms, default iblock and oblock scalarField_real, scalarField_fourier, & ! input data, output data - MPI_COMM_WORLD, fftw_planner_flag) ! use all processors, planer precision + PETSC_COMM_WORLD, fftw_planner_flag) ! use all processors, planer precision if (.not. C_ASSOCIATED(planScalarForth)) call IO_error(810, ext_msg='planScalarForth') - planScalarBack = fftw_mpi_plan_many_dft_c2r(3, [gridFFTW(3),gridFFTW(2),gridFFTW(1)], & ! dimension, logical length in each dimension in reversed order, no. of transforms + planScalarBack = fftw_mpi_plan_many_dft_c2r(3, [gridFFTW(3),gridFFTW(2),gridFFTW(1)], & ! dimension, logical length in each dimension in reversed order, no. of transforms scalarSize, FFTW_MPI_DEFAULT_BLOCK, FFTW_MPI_DEFAULT_BLOCK, & ! no. of transforms, default iblock and oblock - scalarField_fourier,scalarField_real, & ! input data, output data - MPI_COMM_WORLD, fftw_planner_flag) ! use all processors, planer precision + scalarField_fourier,scalarField_real, & ! input data, output data + PETSC_COMM_WORLD, fftw_planner_flag) ! use all processors, planer precision if (.not. C_ASSOCIATED(planScalarBack)) call IO_error(810, ext_msg='planScalarBack') !-------------------------------------------------------------------------------------------------- @@ -699,8 +703,8 @@ real(pReal) function utilities_curlRMS() curl_fourier(l,3) = (+tensorField_fourier(l,2,i,j,k)*xi1st(1,i,j,k)*rescaledGeom(1) & -tensorField_fourier(l,1,i,j,k)*xi1st(2,i,j,k)*rescaledGeom(2)) enddo - utilities_curlRMS = utilities_curlRMS + & - 2.0_pReal*sum(real(curl_fourier)**2.0_pReal + aimag(curl_fourier)**2.0_pReal)! Has somewhere a conj. complex counterpart. Therefore count it twice. + utilities_curlRMS = utilities_curlRMS & + +2.0_pReal*sum(real(curl_fourier)**2.0_pReal+aimag(curl_fourier)**2.0_pReal) ! Has somewhere a conj. complex counterpart. Therefore count it twice. enddo do l = 1_pInt, 3_pInt curl_fourier = (+tensorField_fourier(l,3,1,j,k)*xi1st(2,1,j,k)*rescaledGeom(2) & @@ -710,8 +714,8 @@ real(pReal) function utilities_curlRMS() curl_fourier = (+tensorField_fourier(l,2,1,j,k)*xi1st(1,1,j,k)*rescaledGeom(1) & -tensorField_fourier(l,1,1,j,k)*xi1st(2,1,j,k)*rescaledGeom(2)) enddo - utilities_curlRMS = utilities_curlRMS + & - sum(real(curl_fourier)**2.0_pReal + aimag(curl_fourier)**2.0_pReal)! this layer (DC) does not have a conjugate complex counterpart (if grid(1) /= 1) + utilities_curlRMS = utilities_curlRMS & + + sum(real(curl_fourier)**2.0_pReal + aimag(curl_fourier)**2.0_pReal) ! this layer (DC) does not have a conjugate complex counterpart (if grid(1) /= 1) do l = 1_pInt, 3_pInt curl_fourier = (+tensorField_fourier(l,3,grid1Red,j,k)*xi1st(2,grid1Red,j,k)*rescaledGeom(2) & -tensorField_fourier(l,2,grid1Red,j,k)*xi1st(3,grid1Red,j,k)*rescaledGeom(3)) @@ -720,14 +724,14 @@ real(pReal) function utilities_curlRMS() curl_fourier = (+tensorField_fourier(l,2,grid1Red,j,k)*xi1st(1,grid1Red,j,k)*rescaledGeom(1) & -tensorField_fourier(l,1,grid1Red,j,k)*xi1st(2,grid1Red,j,k)*rescaledGeom(2)) enddo - utilities_curlRMS = utilities_curlRMS + & - sum(real(curl_fourier)**2.0_pReal + aimag(curl_fourier)**2.0_pReal)! this layer (Nyquist) does not have a conjugate complex counterpart (if grid(1) /= 1) + utilities_curlRMS = utilities_curlRMS & + + sum(real(curl_fourier)**2.0_pReal + aimag(curl_fourier)**2.0_pReal) ! this layer (Nyquist) does not have a conjugate complex counterpart (if grid(1) /= 1) enddo; enddo call MPI_Allreduce(MPI_IN_PLACE,utilities_curlRMS,1,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD,ierr) if(ierr /=0_pInt) call IO_error(894_pInt, ext_msg='utilities_curlRMS') utilities_curlRMS = sqrt(utilities_curlRMS) * wgt - if(grid(1) == 1_pInt) utilities_curlRMS = utilities_curlRMS * 0.5_pReal ! counted twice in case of grid(1) == 1 + if(grid(1) == 1_pInt) utilities_curlRMS = utilities_curlRMS * 0.5_pReal ! counted twice in case of grid(1) == 1 end function utilities_curlRMS diff --git a/code/thermal_adiabatic.f90 b/code/thermal_adiabatic.f90 index f4f01f3c3..a0626af62 100644 --- a/code/thermal_adiabatic.f90 +++ b/code/thermal_adiabatic.f90 @@ -74,8 +74,6 @@ subroutine thermal_adiabatic_init(fileUnit) temperature, & temperatureRate, & material_partHomogenization - use numerics,only: & - worldrank implicit none integer(pInt), intent(in) :: fileUnit @@ -88,11 +86,9 @@ subroutine thermal_adiabatic_init(fileUnit) tag = '', & line = '' - mainProcess: if (worldrank == 0) then - write(6,'(/,a)') ' <<<+- thermal_'//THERMAL_ADIABATIC_label//' init -+>>>' - write(6,'(a15,a)') ' Current time: ',IO_timeStamp() + write(6,'(/,a)') ' <<<+- thermal_'//THERMAL_ADIABATIC_label//' init -+>>>' + write(6,'(a15,a)') ' Current time: ',IO_timeStamp() #include "compilation_info.f90" - endif mainProcess maxNinstance = int(count(thermal_type == THERMAL_adiabatic_ID),pInt) if (maxNinstance == 0_pInt) return diff --git a/code/thermal_conduction.f90 b/code/thermal_conduction.f90 index 38e642be4..973ae2d03 100644 --- a/code/thermal_conduction.f90 +++ b/code/thermal_conduction.f90 @@ -75,8 +75,6 @@ subroutine thermal_conduction_init(fileUnit) temperature, & temperatureRate, & material_partHomogenization - use numerics,only: & - worldrank implicit none integer(pInt), intent(in) :: fileUnit @@ -89,11 +87,9 @@ subroutine thermal_conduction_init(fileUnit) tag = '', & line = '' - mainProcess: if (worldrank == 0) then - write(6,'(/,a)') ' <<<+- thermal_'//THERMAL_CONDUCTION_label//' init -+>>>' - write(6,'(a15,a)') ' Current time: ',IO_timeStamp() + write(6,'(/,a)') ' <<<+- thermal_'//THERMAL_CONDUCTION_label//' init -+>>>' + write(6,'(a15,a)') ' Current time: ',IO_timeStamp() #include "compilation_info.f90" - endif mainProcess maxNinstance = int(count(thermal_type == THERMAL_conduction_ID),pInt) if (maxNinstance == 0_pInt) return diff --git a/code/thermal_isothermal.f90 b/code/thermal_isothermal.f90 index a3ac8f685..30ca7562a 100644 --- a/code/thermal_isothermal.f90 +++ b/code/thermal_isothermal.f90 @@ -23,8 +23,6 @@ subroutine thermal_isothermal_init() use IO, only: & IO_timeStamp use material - use numerics, only: & - worldrank implicit none integer(pInt) :: & @@ -32,13 +30,11 @@ subroutine thermal_isothermal_init() NofMyHomog, & sizeState - mainProcess: if (worldrank == 0) then - write(6,'(/,a)') ' <<<+- thermal_'//THERMAL_isothermal_label//' init -+>>>' - write(6,'(a15,a)') ' Current time: ',IO_timeStamp() + write(6,'(/,a)') ' <<<+- thermal_'//THERMAL_isothermal_label//' init -+>>>' + write(6,'(a15,a)') ' Current time: ',IO_timeStamp() #include "compilation_info.f90" - endif mainProcess - initializeInstances: do homog = 1_pInt, material_Nhomogenization + initializeInstances: do homog = 1_pInt, material_Nhomogenization myhomog: if (thermal_type(homog) == THERMAL_isothermal_ID) then NofMyHomog = count(material_homog == homog) diff --git a/code/vacancyflux_cahnhilliard.f90 b/code/vacancyflux_cahnhilliard.f90 index 5abc923dc..f73f66631 100644 --- a/code/vacancyflux_cahnhilliard.f90 +++ b/code/vacancyflux_cahnhilliard.f90 @@ -90,8 +90,6 @@ subroutine vacancyflux_cahnhilliard_init(fileUnit) vacancyflux_initialCv, & material_partHomogenization, & material_partPhase - use numerics,only: & - worldrank implicit none integer(pInt), intent(in) :: fileUnit @@ -104,11 +102,9 @@ subroutine vacancyflux_cahnhilliard_init(fileUnit) tag = '', & line = '' - mainProcess: if (worldrank == 0) then - write(6,'(/,a)') ' <<<+- vacancyflux_'//VACANCYFLUX_cahnhilliard_label//' init -+>>>' - write(6,'(a15,a)') ' Current time: ',IO_timeStamp() + write(6,'(/,a)') ' <<<+- vacancyflux_'//VACANCYFLUX_cahnhilliard_label//' init -+>>>' + write(6,'(a15,a)') ' Current time: ',IO_timeStamp() #include "compilation_info.f90" - endif mainProcess maxNinstance = int(count(vacancyflux_type == VACANCYFLUX_cahnhilliard_ID),pInt) if (maxNinstance == 0_pInt) return diff --git a/code/vacancyflux_isochempot.f90 b/code/vacancyflux_isochempot.f90 index b5b060a18..642d5a2e0 100644 --- a/code/vacancyflux_isochempot.f90 +++ b/code/vacancyflux_isochempot.f90 @@ -72,8 +72,6 @@ subroutine vacancyflux_isochempot_init(fileUnit) vacancyConcRate, & vacancyflux_initialCv, & material_partHomogenization - use numerics,only: & - worldrank implicit none integer(pInt), intent(in) :: fileUnit @@ -86,11 +84,9 @@ subroutine vacancyflux_isochempot_init(fileUnit) tag = '', & line = '' - mainProcess: if (worldrank == 0) then - write(6,'(/,a)') ' <<<+- vacancyflux_'//VACANCYFLUX_isochempot_label//' init -+>>>' - write(6,'(a15,a)') ' Current time: ',IO_timeStamp() + write(6,'(/,a)') ' <<<+- vacancyflux_'//VACANCYFLUX_isochempot_label//' init -+>>>' + write(6,'(a15,a)') ' Current time: ',IO_timeStamp() #include "compilation_info.f90" - endif mainProcess maxNinstance = int(count(vacancyflux_type == VACANCYFLUX_isochempot_ID),pInt) if (maxNinstance == 0_pInt) return diff --git a/code/vacancyflux_isoconc.f90 b/code/vacancyflux_isoconc.f90 index c32cb648d..e4c20b246 100644 --- a/code/vacancyflux_isoconc.f90 +++ b/code/vacancyflux_isoconc.f90 @@ -23,21 +23,17 @@ subroutine vacancyflux_isoconc_init() use IO, only: & IO_timeStamp use material - use numerics, only: & - worldrank implicit none integer(pInt) :: & homog, & NofMyHomog - mainProcess: if (worldrank == 0) then - write(6,'(/,a)') ' <<<+- vacancyflux_'//VACANCYFLUX_isoconc_label//' init -+>>>' - write(6,'(a15,a)') ' Current time: ',IO_timeStamp() + write(6,'(/,a)') ' <<<+- vacancyflux_'//VACANCYFLUX_isoconc_label//' init -+>>>' + write(6,'(a15,a)') ' Current time: ',IO_timeStamp() #include "compilation_info.f90" - endif mainProcess - initializeInstances: do homog = 1_pInt, material_Nhomogenization + initializeInstances: do homog = 1_pInt, material_Nhomogenization myhomog: if (vacancyflux_type(homog) == VACANCYFLUX_isoconc_ID) then NofMyHomog = count(material_homog == homog)