From 3d51dd36faa1480b7b13b148f4c150cd46f83484 Mon Sep 17 00:00:00 2001 From: Christoph Kords Date: Tue, 29 Mar 2011 07:27:19 +0000 Subject: [PATCH] * Introduced preprocessor directives in order to suppress compilation of most write statements when using openmp. This tremendously improves efficiency of parallelization. * Also added some more openmp directives to increase percentage of parallelized code. * "implicit none" was missing in two subroutines of homogenization and constitutive. --- code/CPFEM.f90 | 36 ++- code/constitutive.f90 | 570 +++++++++++++++++---------------- code/constitutive_nonlocal.f90 | 129 ++++---- code/crystallite.f90 | 440 +++++++++++++------------ code/homogenization.f90 | 368 ++++++++++----------- code/mpie_cpfem_abaqus_exp.f | 40 +-- code/mpie_cpfem_abaqus_std.f | 40 +-- code/mpie_cpfem_marc.f90 | 40 +-- 8 files changed, 842 insertions(+), 821 deletions(-) diff --git a/code/CPFEM.f90 b/code/CPFEM.f90 index 9184c5938..3329c60bd 100644 --- a/code/CPFEM.f90 +++ b/code/CPFEM.f90 @@ -373,12 +373,14 @@ subroutine CPFEM_general(mode, ffn, ffn1, Temperature, dt, element, IP, cauchySt endif !$OMP END CRITICAL (write2out) endif - do k = 1,mesh_NcpElems - do j = 1,mesh_maxNips - if (homogenization_sizeState(j,k) > 0_pInt) & - homogenization_state0(j,k)%p = homogenization_state(j,k)%p ! internal state of homogenization scheme + !$OMP PARALLEL DO + do k = 1,mesh_NcpElems + do j = 1,mesh_maxNips + if (homogenization_sizeState(j,k) > 0_pInt) & + homogenization_state0(j,k)%p = homogenization_state(j,k)%p ! internal state of homogenization scheme + enddo enddo - enddo + !$OMP END PARALLEL DO ! *** dump the last converged values of each essential variable to a binary file @@ -492,16 +494,18 @@ subroutine CPFEM_general(mode, ffn, ffn1, Temperature, dt, element, IP, cauchySt endif call materialpoint_stressAndItsTangent(updateJaco, dt) ! calculate stress and its tangent (parallel execution inside) call materialpoint_postResults(dt) ! post results - do e = FEsolving_execElem(1),FEsolving_execElem(2) ! loop over all parallely processed elements - if (microstructure_elemhomo(mesh_element(4,e))) then ! dealing with homogeneous element? - forall (i = 2:FE_Nips(mesh_element(2,e))) ! copy results of first IP to all others - materialpoint_P(1:3,1:3,i,e) = materialpoint_P(1:3,1:3,1,e) - materialpoint_F(1:3,1:3,i,e) = materialpoint_F(1:3,1:3,1,e) - materialpoint_dPdF(1:3,1:3,1:3,1:3,i,e) = materialpoint_dPdF(1:3,1:3,1:3,1:3,1,e) - materialpoint_results(1:materialpoint_sizeResults,i,e) = materialpoint_results(1:materialpoint_sizeResults,1,e) - end forall - endif - enddo + !$OMP PARALLEL DO + do e = FEsolving_execElem(1),FEsolving_execElem(2) ! loop over all parallely processed elements + if (microstructure_elemhomo(mesh_element(4,e))) then ! dealing with homogeneous element? + forall (i = 2:FE_Nips(mesh_element(2,e))) ! copy results of first IP to all others + materialpoint_P(1:3,1:3,i,e) = materialpoint_P(1:3,1:3,1,e) + materialpoint_F(1:3,1:3,i,e) = materialpoint_F(1:3,1:3,1,e) + materialpoint_dPdF(1:3,1:3,1:3,1:3,i,e) = materialpoint_dPdF(1:3,1:3,1:3,1:3,1,e) + materialpoint_results(1:materialpoint_sizeResults,i,e) = materialpoint_results(1:materialpoint_sizeResults,1,e) + end forall + endif + enddo + !$OMP END PARALLEL DO CPFEM_calc_done = .true. endif @@ -511,7 +515,7 @@ subroutine CPFEM_general(mode, ffn, ffn1, Temperature, dt, element, IP, cauchySt CPFEM_cs(1:6,IP,cp_en) = rnd * CPFEM_odd_stress CPFEM_dcsde(1:6,1:6,IP,cp_en) = CPFEM_odd_jacobian * math_identity2nd(6) else - ! translate from P to CS + ! translate from P to CS Kirchhoff = math_mul33x33(materialpoint_P(1:3,1:3,IP, cp_en), math_transpose3x3(materialpoint_F(1:3,1:3,IP,cp_en))) J_inverse = 1.0_pReal / math_det3x3(materialpoint_F(1:3,1:3,IP,cp_en)) CPFEM_cs(1:6,IP,cp_en) = math_Mandel33to6(J_inverse * Kirchhoff) diff --git a/code/constitutive.f90 b/code/constitutive.f90 index 967bf44c8..0d0f866df 100644 --- a/code/constitutive.f90 +++ b/code/constitutive.f90 @@ -10,27 +10,27 @@ MODULE constitutive !*** Include other modules *** - use prec - implicit none +use prec +implicit none - type(p_vec), dimension(:,:,:), allocatable :: constitutive_state0, & ! pointer array to microstructure at start of FE inc - constitutive_partionedState0, & ! pointer array to microstructure at start of homogenization inc - constitutive_subState0, & ! pointer array to microstructure at start of crystallite inc - constitutive_state, & ! pointer array to current microstructure (end of converged time step) - constitutive_state_backup, & ! pointer array to backed up microstructure (end of converged time step) - constitutive_dotState, & ! pointer array to evolution of current microstructure - constitutive_previousDotState,& ! pointer array to previous evolution of current microstructure - constitutive_previousDotState2,&! pointer array to 2nd previous evolution of current microstructure - constitutive_dotState_backup, & ! pointer array to backed up evolution of current microstructure - constitutive_RK4dotState, & ! pointer array to evolution of microstructure defined by classical Runge-Kutta method - constitutive_aTolState ! pointer array to absolute state tolerance - type(p_vec), dimension(:,:,:,:), allocatable :: constitutive_RKCK45dotState ! pointer array to evolution of microstructure used by Cash-Karp Runge-Kutta method - integer(pInt), dimension(:,:,:), allocatable :: constitutive_sizeDotState, & ! size of dotState array - constitutive_sizeState, & ! size of state array per grain - constitutive_sizePostResults ! size of postResults array per grain - integer(pInt) constitutive_maxSizeDotState, & - constitutive_maxSizeState, & - constitutive_maxSizePostResults +type(p_vec), dimension(:,:,:), allocatable :: constitutive_state0, & ! pointer array to microstructure at start of FE inc + constitutive_partionedState0, & ! pointer array to microstructure at start of homogenization inc + constitutive_subState0, & ! pointer array to microstructure at start of crystallite inc + constitutive_state, & ! pointer array to current microstructure (end of converged time step) + constitutive_state_backup, & ! pointer array to backed up microstructure (end of converged time step) + constitutive_dotState, & ! pointer array to evolution of current microstructure + constitutive_previousDotState,& ! pointer array to previous evolution of current microstructure + constitutive_previousDotState2,&! pointer array to 2nd previous evolution of current microstructure + constitutive_dotState_backup, & ! pointer array to backed up evolution of current microstructure + constitutive_RK4dotState, & ! pointer array to evolution of microstructure defined by classical Runge-Kutta method + constitutive_aTolState ! pointer array to absolute state tolerance +type(p_vec), dimension(:,:,:,:), allocatable :: constitutive_RKCK45dotState ! pointer array to evolution of microstructure used by Cash-Karp Runge-Kutta method +integer(pInt), dimension(:,:,:), allocatable :: constitutive_sizeDotState, & ! size of dotState array + constitutive_sizeState, & ! size of state array per grain + constitutive_sizePostResults ! size of postResults array per grain +integer(pInt) constitutive_maxSizeDotState, & + constitutive_maxSizeState, & + constitutive_maxSizePostResults CONTAINS !**************************************** @@ -44,273 +44,299 @@ CONTAINS !* - constitutive_postResults !**************************************** -subroutine constitutive_init() + !************************************** !* Module initialization * !************************************** - use prec, only: pReal,pInt - use debug, only: debug_verbosity, debug_selectiveDebugger, debug_e, debug_i, debug_g - use numerics, only: numerics_integrator - use IO, only: IO_error, IO_open_file, IO_open_jobFile - use mesh, only: mesh_maxNips,mesh_NcpElems,mesh_element,FE_Nips - use material - use constitutive_j2 - use constitutive_phenopowerlaw - use constitutive_titanmod - use constitutive_dislotwin - use constitutive_nonlocal +subroutine constitutive_init() +use prec, only: pReal,pInt +use debug, only: debug_verbosity, debug_selectiveDebugger, debug_e, debug_i, debug_g +use numerics, only: numerics_integrator +use IO, only: IO_error, IO_open_file, IO_open_jobFile +use mesh, only: mesh_maxNips,mesh_NcpElems,mesh_element,FE_Nips +use material +use constitutive_j2 +use constitutive_phenopowerlaw +use constitutive_titanmod +use constitutive_dislotwin +use constitutive_nonlocal - integer(pInt), parameter :: fileunit = 200 - integer(pInt) e,i,g,p,s,myInstance,myNgrains - integer(pInt), dimension(:,:), pointer :: thisSize - character(len=64), dimension(:,:), pointer :: thisOutput - logical knownConstitution - - if(.not. IO_open_file(fileunit,material_configFile)) call IO_error (100) ! corrupt config file +implicit none - call constitutive_j2_init(fileunit) ! parse all phases of this constitution - call constitutive_phenopowerlaw_init(fileunit) - call constitutive_titanmod_init(fileunit) - call constitutive_dislotwin_init(fileunit) - call constitutive_nonlocal_init(fileunit) +integer(pInt), parameter :: fileunit = 200 +integer(pInt) g, & ! grain number + i, & ! integration point number + e, & ! element number + gMax, & ! maximum number of grains + iMax, & ! maximum number of integration points + eMax, & ! maximum number of elements + p, & + s, & + myInstance,& + myNgrains +integer(pInt), dimension(:,:), pointer :: thisSize +character(len=64), dimension(:,:), pointer :: thisOutput +logical knownConstitution - close(fileunit) -! write description file for constitutive phase output +! --- PARSE CONSTITUTIONS FROM CONFIG FILE --- - if(.not. IO_open_jobFile(fileunit,'outputConstitutive')) call IO_error (50) ! problems in writing file - - do p = 1,material_Nphase - i = phase_constitutionInstance(p) ! which instance of a constitution is present phase - knownConstitution = .true. ! assume valid - select case(phase_constitution(p)) ! split per constitiution - case (constitutive_j2_label) - thisOutput => constitutive_j2_output - thisSize => constitutive_j2_sizePostResult - case (constitutive_phenopowerlaw_label) - thisOutput => constitutive_phenopowerlaw_output - thisSize => constitutive_phenopowerlaw_sizePostResult - case (constitutive_titanmod_label) - thisOutput => constitutive_titanmod_output - thisSize => constitutive_titanmod_sizePostResult - case (constitutive_dislotwin_label) - thisOutput => constitutive_dislotwin_output - thisSize => constitutive_dislotwin_sizePostResult - case (constitutive_nonlocal_label) - thisOutput => constitutive_nonlocal_output - thisSize => constitutive_nonlocal_sizePostResult - case default - knownConstitution = .false. - end select +if(.not. IO_open_file(fileunit,material_configFile)) call IO_error (100) ! corrupt config file +call constitutive_j2_init(fileunit) +call constitutive_phenopowerlaw_init(fileunit) +call constitutive_titanmod_init(fileunit) +call constitutive_dislotwin_init(fileunit) +call constitutive_nonlocal_init(fileunit) +close(fileunit) - write(fileunit,*) - write(fileunit,'(a)') '['//trim(phase_name(p))//']' - write(fileunit,*) - if (knownConstitution) then - write(fileunit,'(a)') '(constitution)'//char(9)//trim(phase_constitution(p)) - do e = 1,phase_Noutput(p) - write(fileunit,'(a,i4)') trim(thisOutput(e,i))//char(9),thisSize(e,i) - enddo - endif - enddo - close(fileunit) +! --- WRITE DESCRIPTION FILE FOR CONSTITUTIVE PHASE OUTPUT --- -! allocate memory for state management +if(.not. IO_open_jobFile(fileunit,'outputConstitutive')) then ! problems in writing file + call IO_error (50) +endif +do p = 1,material_Nphase + i = phase_constitutionInstance(p) ! which instance of a constitution is present phase + knownConstitution = .true. ! assume valid + select case(phase_constitution(p)) ! split per constitiution + case (constitutive_j2_label) + thisOutput => constitutive_j2_output + thisSize => constitutive_j2_sizePostResult + case (constitutive_phenopowerlaw_label) + thisOutput => constitutive_phenopowerlaw_output + thisSize => constitutive_phenopowerlaw_sizePostResult + case (constitutive_titanmod_label) + thisOutput => constitutive_titanmod_output + thisSize => constitutive_titanmod_sizePostResult + case (constitutive_dislotwin_label) + thisOutput => constitutive_dislotwin_output + thisSize => constitutive_dislotwin_sizePostResult + case (constitutive_nonlocal_label) + thisOutput => constitutive_nonlocal_output + thisSize => constitutive_nonlocal_sizePostResult + case default + knownConstitution = .false. + end select + write(fileunit,*) + write(fileunit,'(a)') '['//trim(phase_name(p))//']' + write(fileunit,*) + if (knownConstitution) then + write(fileunit,'(a)') '(constitution)'//char(9)//trim(phase_constitution(p)) + do e = 1,phase_Noutput(p) + write(fileunit,'(a,i4)') trim(thisOutput(e,i))//char(9),thisSize(e,i) + enddo + endif +enddo +close(fileunit) - allocate(constitutive_state0(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems)) - allocate(constitutive_partionedState0(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems)) - allocate(constitutive_subState0(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems)) - allocate(constitutive_state(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems)) - allocate(constitutive_state_backup(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems)) - allocate(constitutive_dotState(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems)) - allocate(constitutive_dotState_backup(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems)) - allocate(constitutive_aTolState(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems)) - allocate(constitutive_sizeDotState(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems)) ; constitutive_sizeDotState = 0_pInt - allocate(constitutive_sizeState(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems)) ; constitutive_sizeState = 0_pInt - allocate(constitutive_sizePostResults(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems)); constitutive_sizePostResults = 0_pInt - if (any(numerics_integrator == 1)) then - allocate(constitutive_previousDotState(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems)) - allocate(constitutive_previousDotState2(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems)) - endif - if (any(numerics_integrator == 4)) & - allocate(constitutive_RK4dotState(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems)) - if (any(numerics_integrator == 5)) & - allocate(constitutive_RKCK45dotState(6,homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems)) - - do e = 1,mesh_NcpElems ! loop over elements - myNgrains = homogenization_Ngrains(mesh_element(3,e)) - do i = 1,FE_Nips(mesh_element(2,e)) ! loop over IPs - do g = 1,myNgrains ! loop over grains - myInstance = phase_constitutionInstance(material_phase(g,i,e)) - select case(phase_constitution(material_phase(g,i,e))) - - case (constitutive_j2_label) - allocate(constitutive_state0(g,i,e)%p(constitutive_j2_sizeState(myInstance))) - allocate(constitutive_partionedState0(g,i,e)%p(constitutive_j2_sizeState(myInstance))) - allocate(constitutive_subState0(g,i,e)%p(constitutive_j2_sizeState(myInstance))) - allocate(constitutive_state(g,i,e)%p(constitutive_j2_sizeState(myInstance))) - allocate(constitutive_state_backup(g,i,e)%p(constitutive_j2_sizeState(myInstance))) - allocate(constitutive_aTolState(g,i,e)%p(constitutive_j2_sizeState(myInstance))) - allocate(constitutive_dotState(g,i,e)%p(constitutive_j2_sizeDotState(myInstance))) - allocate(constitutive_dotState_backup(g,i,e)%p(constitutive_j2_sizeDotState(myInstance))) - if (any(numerics_integrator == 1)) then - allocate(constitutive_previousDotState(g,i,e)%p(constitutive_j2_sizeDotState(myInstance))) - allocate(constitutive_previousDotState2(g,i,e)%p(constitutive_j2_sizeDotState(myInstance))) - endif - if (any(numerics_integrator == 4)) & - allocate(constitutive_RK4dotState(g,i,e)%p(constitutive_j2_sizeDotState(myInstance))) - if (any(numerics_integrator == 5)) then - do s = 1,6 - allocate(constitutive_RKCK45dotState(s,g,i,e)%p(constitutive_j2_sizeDotState(myInstance))) - enddo - endif - constitutive_state0(g,i,e)%p = constitutive_j2_stateInit(myInstance) - constitutive_aTolState(g,i,e)%p = constitutive_j2_aTolState(myInstance) - constitutive_sizeState(g,i,e) = constitutive_j2_sizeState(myInstance) - constitutive_sizeDotState(g,i,e) = constitutive_j2_sizeDotState(myInstance) - constitutive_sizePostResults(g,i,e) = constitutive_j2_sizePostResults(myInstance) + +! --- ALLOCATION OF STATES --- + +gMax = homogenization_maxNgrains +iMax = mesh_maxNips +eMax = mesh_NcpElems + +allocate(constitutive_state0(gMax,iMax,eMax)) +allocate(constitutive_partionedState0(gMax,iMax,eMax)) +allocate(constitutive_subState0(gMax,iMax,eMax)) +allocate(constitutive_state(gMax,iMax,eMax)) +allocate(constitutive_state_backup(gMax,iMax,eMax)) +allocate(constitutive_dotState(gMax,iMax,eMax)) +allocate(constitutive_dotState_backup(gMax,iMax,eMax)) +allocate(constitutive_aTolState(gMax,iMax,eMax)) +allocate(constitutive_sizeDotState(gMax,iMax,eMax)) ; constitutive_sizeDotState = 0_pInt +allocate(constitutive_sizeState(gMax,iMax,eMax)) ; constitutive_sizeState = 0_pInt +allocate(constitutive_sizePostResults(gMax,iMax,eMax)); constitutive_sizePostResults = 0_pInt +if (any(numerics_integrator == 1)) then + allocate(constitutive_previousDotState(gMax,iMax,eMax)) + allocate(constitutive_previousDotState2(gMax,iMax,eMax)) +endif +if (any(numerics_integrator == 4)) then + allocate(constitutive_RK4dotState(gMax,iMax,eMax)) +endif +if (any(numerics_integrator == 5)) then + allocate(constitutive_RKCK45dotState(6,gMax,iMax,eMax)) +endif + +!$OMP PARALLEL DO PRIVATE(myNgrains,myInstance) + do e = 1,mesh_NcpElems ! loop over elements + myNgrains = homogenization_Ngrains(mesh_element(3,e)) + do i = 1,FE_Nips(mesh_element(2,e)) ! loop over IPs + do g = 1,myNgrains ! loop over grains + myInstance = phase_constitutionInstance(material_phase(g,i,e)) + select case(phase_constitution(material_phase(g,i,e))) + + case (constitutive_j2_label) + allocate(constitutive_state0(g,i,e)%p(constitutive_j2_sizeState(myInstance))) + allocate(constitutive_partionedState0(g,i,e)%p(constitutive_j2_sizeState(myInstance))) + allocate(constitutive_subState0(g,i,e)%p(constitutive_j2_sizeState(myInstance))) + allocate(constitutive_state(g,i,e)%p(constitutive_j2_sizeState(myInstance))) + allocate(constitutive_state_backup(g,i,e)%p(constitutive_j2_sizeState(myInstance))) + allocate(constitutive_aTolState(g,i,e)%p(constitutive_j2_sizeState(myInstance))) + allocate(constitutive_dotState(g,i,e)%p(constitutive_j2_sizeDotState(myInstance))) + allocate(constitutive_dotState_backup(g,i,e)%p(constitutive_j2_sizeDotState(myInstance))) + if (any(numerics_integrator == 1)) then + allocate(constitutive_previousDotState(g,i,e)%p(constitutive_j2_sizeDotState(myInstance))) + allocate(constitutive_previousDotState2(g,i,e)%p(constitutive_j2_sizeDotState(myInstance))) + endif + if (any(numerics_integrator == 4)) then + allocate(constitutive_RK4dotState(g,i,e)%p(constitutive_j2_sizeDotState(myInstance))) + endif + if (any(numerics_integrator == 5)) then + do s = 1,6 + allocate(constitutive_RKCK45dotState(s,g,i,e)%p(constitutive_j2_sizeDotState(myInstance))) + enddo + endif + constitutive_state0(g,i,e)%p = constitutive_j2_stateInit(myInstance) + constitutive_aTolState(g,i,e)%p = constitutive_j2_aTolState(myInstance) + constitutive_sizeState(g,i,e) = constitutive_j2_sizeState(myInstance) + constitutive_sizeDotState(g,i,e) = constitutive_j2_sizeDotState(myInstance) + constitutive_sizePostResults(g,i,e) = constitutive_j2_sizePostResults(myInstance) - case (constitutive_phenopowerlaw_label) - allocate(constitutive_state0(g,i,e)%p(constitutive_phenopowerlaw_sizeState(myInstance))) - allocate(constitutive_partionedState0(g,i,e)%p(constitutive_phenopowerlaw_sizeState(myInstance))) - allocate(constitutive_subState0(g,i,e)%p(constitutive_phenopowerlaw_sizeState(myInstance))) - allocate(constitutive_state(g,i,e)%p(constitutive_phenopowerlaw_sizeState(myInstance))) - allocate(constitutive_state_backup(g,i,e)%p(constitutive_phenopowerlaw_sizeState(myInstance))) - allocate(constitutive_aTolState(g,i,e)%p(constitutive_phenopowerlaw_sizeState(myInstance))) - allocate(constitutive_dotState(g,i,e)%p(constitutive_phenopowerlaw_sizeDotState(myInstance))) - allocate(constitutive_dotState_backup(g,i,e)%p(constitutive_phenopowerlaw_sizeDotState(myInstance))) - if (any(numerics_integrator == 1)) then - allocate(constitutive_previousDotState(g,i,e)%p(constitutive_phenopowerlaw_sizeDotState(myInstance))) - allocate(constitutive_previousDotState2(g,i,e)%p(constitutive_phenopowerlaw_sizeDotState(myInstance))) - endif - if (any(numerics_integrator == 4)) & - allocate(constitutive_RK4dotState(g,i,e)%p(constitutive_phenopowerlaw_sizeDotState(myInstance))) - if (any(numerics_integrator == 5)) then - do s = 1,6 - allocate(constitutive_RKCK45dotState(s,g,i,e)%p(constitutive_phenopowerlaw_sizeDotState(myInstance))) - enddo - endif - constitutive_state0(g,i,e)%p = constitutive_phenopowerlaw_stateInit(myInstance) - constitutive_aTolState(g,i,e)%p = constitutive_phenopowerlaw_aTolState(myInstance) - constitutive_sizeState(g,i,e) = constitutive_phenopowerlaw_sizeState(myInstance) - constitutive_sizeDotState(g,i,e) = constitutive_phenopowerlaw_sizeDotState(myInstance) - constitutive_sizePostResults(g,i,e) = constitutive_phenopowerlaw_sizePostResults(myInstance) + case (constitutive_phenopowerlaw_label) + allocate(constitutive_state0(g,i,e)%p(constitutive_phenopowerlaw_sizeState(myInstance))) + allocate(constitutive_partionedState0(g,i,e)%p(constitutive_phenopowerlaw_sizeState(myInstance))) + allocate(constitutive_subState0(g,i,e)%p(constitutive_phenopowerlaw_sizeState(myInstance))) + allocate(constitutive_state(g,i,e)%p(constitutive_phenopowerlaw_sizeState(myInstance))) + allocate(constitutive_state_backup(g,i,e)%p(constitutive_phenopowerlaw_sizeState(myInstance))) + allocate(constitutive_aTolState(g,i,e)%p(constitutive_phenopowerlaw_sizeState(myInstance))) + allocate(constitutive_dotState(g,i,e)%p(constitutive_phenopowerlaw_sizeDotState(myInstance))) + allocate(constitutive_dotState_backup(g,i,e)%p(constitutive_phenopowerlaw_sizeDotState(myInstance))) + if (any(numerics_integrator == 1)) then + allocate(constitutive_previousDotState(g,i,e)%p(constitutive_phenopowerlaw_sizeDotState(myInstance))) + allocate(constitutive_previousDotState2(g,i,e)%p(constitutive_phenopowerlaw_sizeDotState(myInstance))) + endif + if (any(numerics_integrator == 4)) then + allocate(constitutive_RK4dotState(g,i,e)%p(constitutive_phenopowerlaw_sizeDotState(myInstance))) + endif + if (any(numerics_integrator == 5)) then + do s = 1,6 + allocate(constitutive_RKCK45dotState(s,g,i,e)%p(constitutive_phenopowerlaw_sizeDotState(myInstance))) + enddo + endif + constitutive_state0(g,i,e)%p = constitutive_phenopowerlaw_stateInit(myInstance) + constitutive_aTolState(g,i,e)%p = constitutive_phenopowerlaw_aTolState(myInstance) + constitutive_sizeState(g,i,e) = constitutive_phenopowerlaw_sizeState(myInstance) + constitutive_sizeDotState(g,i,e) = constitutive_phenopowerlaw_sizeDotState(myInstance) + constitutive_sizePostResults(g,i,e) = constitutive_phenopowerlaw_sizePostResults(myInstance) + + case (constitutive_titanmod_label) + allocate(constitutive_state0(g,i,e)%p(constitutive_titanmod_sizeState(myInstance))) + allocate(constitutive_partionedState0(g,i,e)%p(constitutive_titanmod_sizeState(myInstance))) + allocate(constitutive_subState0(g,i,e)%p(constitutive_titanmod_sizeState(myInstance))) + allocate(constitutive_state(g,i,e)%p(constitutive_titanmod_sizeState(myInstance))) + allocate(constitutive_state_backup(g,i,e)%p(constitutive_titanmod_sizeState(myInstance))) + allocate(constitutive_aTolState(g,i,e)%p(constitutive_titanmod_sizeState(myInstance))) + allocate(constitutive_dotState(g,i,e)%p(constitutive_titanmod_sizeDotState(myInstance))) + allocate(constitutive_dotState_backup(g,i,e)%p(constitutive_titanmod_sizeDotState(myInstance))) + if (any(numerics_integrator == 1)) then + allocate(constitutive_previousDotState(g,i,e)%p(constitutive_titanmod_sizeDotState(myInstance))) + allocate(constitutive_previousDotState2(g,i,e)%p(constitutive_titanmod_sizeDotState(myInstance))) + endif + if (any(numerics_integrator == 4)) then + allocate(constitutive_RK4dotState(g,i,e)%p(constitutive_titanmod_sizeDotState(myInstance))) + endif + if (any(numerics_integrator == 5)) then + do s = 1,6 + allocate(constitutive_RKCK45dotState(s,g,i,e)%p(constitutive_titanmod_sizeDotState(myInstance))) + enddo + endif + constitutive_state0(g,i,e)%p = constitutive_titanmod_stateInit(myInstance) + constitutive_aTolState(g,i,e)%p = constitutive_titanmod_aTolState(myInstance) + constitutive_sizeState(g,i,e) = constitutive_titanmod_sizeState(myInstance) + constitutive_sizeDotState(g,i,e) = constitutive_titanmod_sizeDotState(myInstance) + constitutive_sizePostResults(g,i,e) = constitutive_titanmod_sizePostResults(myInstance) + + case (constitutive_dislotwin_label) + allocate(constitutive_state0(g,i,e)%p(constitutive_dislotwin_sizeState(myInstance))) + allocate(constitutive_partionedState0(g,i,e)%p(constitutive_dislotwin_sizeState(myInstance))) + allocate(constitutive_subState0(g,i,e)%p(constitutive_dislotwin_sizeState(myInstance))) + allocate(constitutive_state(g,i,e)%p(constitutive_dislotwin_sizeState(myInstance))) + allocate(constitutive_state_backup(g,i,e)%p(constitutive_dislotwin_sizeState(myInstance))) + allocate(constitutive_aTolState(g,i,e)%p(constitutive_dislotwin_sizeState(myInstance))) + allocate(constitutive_dotState(g,i,e)%p(constitutive_dislotwin_sizeDotState(myInstance))) + allocate(constitutive_dotState_backup(g,i,e)%p(constitutive_dislotwin_sizeDotState(myInstance))) + if (any(numerics_integrator == 1)) then + allocate(constitutive_previousDotState(g,i,e)%p(constitutive_dislotwin_sizeDotState(myInstance))) + allocate(constitutive_previousDotState2(g,i,e)%p(constitutive_dislotwin_sizeDotState(myInstance))) + endif + if (any(numerics_integrator == 4)) then + allocate(constitutive_RK4dotState(g,i,e)%p(constitutive_dislotwin_sizeDotState(myInstance))) + endif + if (any(numerics_integrator == 5)) then + do s = 1,6 + allocate(constitutive_RKCK45dotState(s,g,i,e)%p(constitutive_dislotwin_sizeDotState(myInstance))) + enddo + endif + constitutive_state0(g,i,e)%p = constitutive_dislotwin_stateInit(myInstance) + constitutive_aTolState(g,i,e)%p = constitutive_dislotwin_aTolState(myInstance) + constitutive_sizeState(g,i,e) = constitutive_dislotwin_sizeState(myInstance) + constitutive_sizeDotState(g,i,e) = constitutive_dislotwin_sizeDotState(myInstance) + constitutive_sizePostResults(g,i,e) = constitutive_dislotwin_sizePostResults(myInstance) + + case (constitutive_nonlocal_label) + allocate(constitutive_state0(g,i,e)%p(constitutive_nonlocal_sizeState(myInstance))) + allocate(constitutive_partionedState0(g,i,e)%p(constitutive_nonlocal_sizeState(myInstance))) + allocate(constitutive_subState0(g,i,e)%p(constitutive_nonlocal_sizeState(myInstance))) + allocate(constitutive_state(g,i,e)%p(constitutive_nonlocal_sizeState(myInstance))) + allocate(constitutive_state_backup(g,i,e)%p(constitutive_nonlocal_sizeState(myInstance))) + allocate(constitutive_aTolState(g,i,e)%p(constitutive_nonlocal_sizeState(myInstance))) + allocate(constitutive_dotState(g,i,e)%p(constitutive_nonlocal_sizeDotState(myInstance))) + allocate(constitutive_dotState_backup(g,i,e)%p(constitutive_nonlocal_sizeDotState(myInstance))) + if (any(numerics_integrator == 1)) then + allocate(constitutive_previousDotState(g,i,e)%p(constitutive_nonlocal_sizeDotState(myInstance))) + allocate(constitutive_previousDotState2(g,i,e)%p(constitutive_nonlocal_sizeDotState(myInstance))) + endif + if (any(numerics_integrator == 4)) then + allocate(constitutive_RK4dotState(g,i,e)%p(constitutive_nonlocal_sizeDotState(myInstance))) + endif + if (any(numerics_integrator == 5)) then + do s = 1,6 + allocate(constitutive_RKCK45dotState(s,g,i,e)%p(constitutive_nonlocal_sizeDotState(myInstance))) + enddo + endif + constitutive_state0(g,i,e)%p = constitutive_nonlocal_stateInit(myInstance) + constitutive_aTolState(g,i,e)%p = constitutive_nonlocal_aTolState(myInstance) + constitutive_sizeState(g,i,e) = constitutive_nonlocal_sizeState(myInstance) + constitutive_sizeDotState(g,i,e) = constitutive_nonlocal_sizeDotState(myInstance) + constitutive_sizePostResults(g,i,e) = constitutive_nonlocal_sizePostResults(myInstance) + + case default + call IO_error(200,material_phase(g,i,e)) ! unknown constitution + + end select + constitutive_partionedState0(g,i,e)%p = constitutive_state0(g,i,e)%p + enddo + enddo + enddo +!$OMP END PARALLEL DO - case (constitutive_titanmod_label) - allocate(constitutive_state0(g,i,e)%p(constitutive_titanmod_sizeState(myInstance))) - allocate(constitutive_partionedState0(g,i,e)%p(constitutive_titanmod_sizeState(myInstance))) - allocate(constitutive_subState0(g,i,e)%p(constitutive_titanmod_sizeState(myInstance))) - allocate(constitutive_state(g,i,e)%p(constitutive_titanmod_sizeState(myInstance))) - allocate(constitutive_state_backup(g,i,e)%p(constitutive_titanmod_sizeState(myInstance))) - allocate(constitutive_aTolState(g,i,e)%p(constitutive_titanmod_sizeState(myInstance))) - allocate(constitutive_dotState(g,i,e)%p(constitutive_titanmod_sizeDotState(myInstance))) - allocate(constitutive_dotState_backup(g,i,e)%p(constitutive_titanmod_sizeDotState(myInstance))) - if (any(numerics_integrator == 1)) then - allocate(constitutive_previousDotState(g,i,e)%p(constitutive_titanmod_sizeDotState(myInstance))) - allocate(constitutive_previousDotState2(g,i,e)%p(constitutive_titanmod_sizeDotState(myInstance))) - endif - if (any(numerics_integrator == 4)) & - allocate(constitutive_RK4dotState(g,i,e)%p(constitutive_titanmod_sizeDotState(myInstance))) - if (any(numerics_integrator == 5)) then - do s = 1,6 - allocate(constitutive_RKCK45dotState(s,g,i,e)%p(constitutive_titanmod_sizeDotState(myInstance))) - enddo - endif - constitutive_state0(g,i,e)%p = constitutive_titanmod_stateInit(myInstance) - constitutive_aTolState(g,i,e)%p = constitutive_titanmod_aTolState(myInstance) - constitutive_sizeState(g,i,e) = constitutive_titanmod_sizeState(myInstance) - constitutive_sizeDotState(g,i,e) = constitutive_titanmod_sizeDotState(myInstance) - constitutive_sizePostResults(g,i,e) = constitutive_titanmod_sizePostResults(myInstance) - - case (constitutive_dislotwin_label) - allocate(constitutive_state0(g,i,e)%p(constitutive_dislotwin_sizeState(myInstance))) - allocate(constitutive_partionedState0(g,i,e)%p(constitutive_dislotwin_sizeState(myInstance))) - allocate(constitutive_subState0(g,i,e)%p(constitutive_dislotwin_sizeState(myInstance))) - allocate(constitutive_state(g,i,e)%p(constitutive_dislotwin_sizeState(myInstance))) - allocate(constitutive_state_backup(g,i,e)%p(constitutive_dislotwin_sizeState(myInstance))) - allocate(constitutive_aTolState(g,i,e)%p(constitutive_dislotwin_sizeState(myInstance))) - allocate(constitutive_dotState(g,i,e)%p(constitutive_dislotwin_sizeDotState(myInstance))) - allocate(constitutive_dotState_backup(g,i,e)%p(constitutive_dislotwin_sizeDotState(myInstance))) - if (any(numerics_integrator == 1)) then - allocate(constitutive_previousDotState(g,i,e)%p(constitutive_dislotwin_sizeDotState(myInstance))) - allocate(constitutive_previousDotState2(g,i,e)%p(constitutive_dislotwin_sizeDotState(myInstance))) - endif - if (any(numerics_integrator == 4)) & - allocate(constitutive_RK4dotState(g,i,e)%p(constitutive_dislotwin_sizeDotState(myInstance))) - if (any(numerics_integrator == 5)) then - do s = 1,6 - allocate(constitutive_RKCK45dotState(s,g,i,e)%p(constitutive_dislotwin_sizeDotState(myInstance))) - enddo - endif - constitutive_state0(g,i,e)%p = constitutive_dislotwin_stateInit(myInstance) - constitutive_aTolState(g,i,e)%p = constitutive_dislotwin_aTolState(myInstance) - constitutive_sizeState(g,i,e) = constitutive_dislotwin_sizeState(myInstance) - constitutive_sizeDotState(g,i,e) = constitutive_dislotwin_sizeDotState(myInstance) - constitutive_sizePostResults(g,i,e) = constitutive_dislotwin_sizePostResults(myInstance) - - case (constitutive_nonlocal_label) - allocate(constitutive_state0(g,i,e)%p(constitutive_nonlocal_sizeState(myInstance))) - allocate(constitutive_partionedState0(g,i,e)%p(constitutive_nonlocal_sizeState(myInstance))) - allocate(constitutive_subState0(g,i,e)%p(constitutive_nonlocal_sizeState(myInstance))) - allocate(constitutive_state(g,i,e)%p(constitutive_nonlocal_sizeState(myInstance))) - allocate(constitutive_state_backup(g,i,e)%p(constitutive_nonlocal_sizeState(myInstance))) - allocate(constitutive_aTolState(g,i,e)%p(constitutive_nonlocal_sizeState(myInstance))) - allocate(constitutive_dotState(g,i,e)%p(constitutive_nonlocal_sizeDotState(myInstance))) - allocate(constitutive_dotState_backup(g,i,e)%p(constitutive_nonlocal_sizeDotState(myInstance))) - if (any(numerics_integrator == 1)) then - allocate(constitutive_previousDotState(g,i,e)%p(constitutive_nonlocal_sizeDotState(myInstance))) - allocate(constitutive_previousDotState2(g,i,e)%p(constitutive_nonlocal_sizeDotState(myInstance))) - endif - if (any(numerics_integrator == 4)) & - allocate(constitutive_RK4dotState(g,i,e)%p(constitutive_nonlocal_sizeDotState(myInstance))) - if (any(numerics_integrator == 5)) then - do s = 1,6 - allocate(constitutive_RKCK45dotState(s,g,i,e)%p(constitutive_nonlocal_sizeDotState(myInstance))) - enddo - endif - constitutive_state0(g,i,e)%p = constitutive_nonlocal_stateInit(myInstance) - constitutive_aTolState(g,i,e)%p = constitutive_nonlocal_aTolState(myInstance) - constitutive_sizeState(g,i,e) = constitutive_nonlocal_sizeState(myInstance) - constitutive_sizeDotState(g,i,e) = constitutive_nonlocal_sizeDotState(myInstance) - constitutive_sizePostResults(g,i,e) = constitutive_nonlocal_sizePostResults(myInstance) - - case default - call IO_error(200,material_phase(g,i,e)) ! unknown constitution - - end select - constitutive_partionedState0(g,i,e)%p = constitutive_state0(g,i,e)%p - enddo - enddo - enddo +constitutive_maxSizeState = maxval(constitutive_sizeState) +constitutive_maxSizeDotState = maxval(constitutive_sizeDotState) +constitutive_maxSizePostResults = maxval(constitutive_sizePostResults) - constitutive_maxSizeState = maxval(constitutive_sizeState) - constitutive_maxSizeDotState = maxval(constitutive_sizeDotState) - constitutive_maxSizePostResults = maxval(constitutive_sizePostResults) - - !$OMP CRITICAL (write2out) - write(6,*) - write(6,*) '<<<+- constitutive init -+>>>' - write(6,*) '$Id$' - write(6,*) - if (debug_verbosity > 0) then - write(6,'(a32,x,7(i5,x))') 'constitutive_state0: ', shape(constitutive_state0) - write(6,'(a32,x,7(i5,x))') 'constitutive_partionedState0: ', shape(constitutive_partionedState0) - write(6,'(a32,x,7(i5,x))') 'constitutive_subState0: ', shape(constitutive_subState0) - write(6,'(a32,x,7(i5,x))') 'constitutive_state: ', shape(constitutive_state) - write(6,'(a32,x,7(i5,x))') 'constitutive_aTolState: ', shape(constitutive_aTolState) - write(6,'(a32,x,7(i5,x))') 'constitutive_dotState: ', shape(constitutive_dotState) - write(6,'(a32,x,7(i5,x))') 'constitutive_sizeState: ', shape(constitutive_sizeState) - write(6,'(a32,x,7(i5,x))') 'constitutive_sizeDotState: ', shape(constitutive_sizeDotState) - write(6,'(a32,x,7(i5,x))') 'constitutive_sizePostResults: ', shape(constitutive_sizePostResults) - write(6,*) - write(6,'(a32,x,7(i5,x))') 'maxSizeState: ', constitutive_maxSizeState - write(6,'(a32,x,7(i5,x))') 'maxSizeDotState: ', constitutive_maxSizeDotState - write(6,'(a32,x,7(i5,x))') 'maxSizePostResults: ', constitutive_maxSizePostResults - endif - call flush(6) - !$OMP END CRITICAL (write2out) - return +!$OMP CRITICAL (write2out) + write(6,*) + write(6,*) '<<<+- constitutive init -+>>>' + write(6,*) '$Id$' + write(6,*) + if (debug_verbosity > 0) then + write(6,'(a32,x,7(i5,x))') 'constitutive_state0: ', shape(constitutive_state0) + write(6,'(a32,x,7(i5,x))') 'constitutive_partionedState0: ', shape(constitutive_partionedState0) + write(6,'(a32,x,7(i5,x))') 'constitutive_subState0: ', shape(constitutive_subState0) + write(6,'(a32,x,7(i5,x))') 'constitutive_state: ', shape(constitutive_state) + write(6,'(a32,x,7(i5,x))') 'constitutive_aTolState: ', shape(constitutive_aTolState) + write(6,'(a32,x,7(i5,x))') 'constitutive_dotState: ', shape(constitutive_dotState) + write(6,'(a32,x,7(i5,x))') 'constitutive_sizeState: ', shape(constitutive_sizeState) + write(6,'(a32,x,7(i5,x))') 'constitutive_sizeDotState: ', shape(constitutive_sizeDotState) + write(6,'(a32,x,7(i5,x))') 'constitutive_sizePostResults: ', shape(constitutive_sizePostResults) + write(6,*) + write(6,'(a32,x,7(i5,x))') 'maxSizeState: ', constitutive_maxSizeState + write(6,'(a32,x,7(i5,x))') 'maxSizeDotState: ', constitutive_maxSizeDotState + write(6,'(a32,x,7(i5,x))') 'maxSizePostResults: ', constitutive_maxSizePostResults + endif + call flush(6) +!$OMP END CRITICAL (write2out) endsubroutine diff --git a/code/constitutive_nonlocal.f90 b/code/constitutive_nonlocal.f90 index 8ad6f9a2e..75290b359 100644 --- a/code/constitutive_nonlocal.f90 +++ b/code/constitutive_nonlocal.f90 @@ -280,6 +280,7 @@ constitutive_nonlocal_interactionSlipSlip = 0.0_pReal constitutive_nonlocal_dLowerEdgePerSlipFamily = 0.0_pReal constitutive_nonlocal_dLowerScrewPerSlipFamily = 0.0_pReal + !*** readout data from material.config file rewind(file) @@ -384,7 +385,8 @@ enddo lattice_initializeStructure(constitutive_nonlocal_structureName(i), constitutive_nonlocal_CoverA(i)) ! our lattice structure is defined in the material.config file by the structureName (and the c/a ratio) myStructure = constitutive_nonlocal_structure(i) -!*** sanity checks + + !*** sanity checks if (myStructure < 1 .or. myStructure > 3) call IO_error(205) if (sum(constitutive_nonlocal_Nslip(:,i)) <= 0_pInt) call IO_error(235,ext_msg='Nslip') @@ -421,7 +423,7 @@ enddo .or. constitutive_nonlocal_surfaceTransmissivity(i) > 1.0_pReal) call IO_error(235,ext_msg='surfaceTransmissivity') -!*** determine total number of active slip systems + !*** determine total number of active slip systems constitutive_nonlocal_Nslip(1:lattice_maxNslipFamily,i) = min( lattice_NslipSystem(1:lattice_maxNslipFamily, myStructure), & constitutive_nonlocal_Nslip(1:lattice_maxNslipFamily,i) ) ! we can't use more slip systems per family than specified in lattice @@ -471,7 +473,8 @@ do i = 1,maxNinstance myStructure = constitutive_nonlocal_structure(i) ! lattice structure of this instance -!*** Inverse lookup of my slip system family and the slip system in lattice + + !*** Inverse lookup of my slip system family and the slip system in lattice l = 0_pInt do f = 1,lattice_maxNslipFamily @@ -482,7 +485,7 @@ do i = 1,maxNinstance enddo; enddo -!*** determine size of state array + !*** determine size of state array ns = constitutive_nonlocal_totalNslip(i) constitutive_nonlocal_sizeState(i) = size(constitutive_nonlocal_listBasicStates) * ns & @@ -490,7 +493,7 @@ do i = 1,maxNinstance constitutive_nonlocal_sizeDotState(i) = size(constitutive_nonlocal_listBasicStates) * ns -!*** determine size of postResults array + !*** determine size of postResults array do o = 1,maxval(phase_Noutput) select case(constitutive_nonlocal_output(o,i)) @@ -572,7 +575,7 @@ do i = 1,maxNinstance enddo -!*** elasticity matrix and shear modulus according to material.config + !*** elasticity matrix and shear modulus according to material.config select case (myStructure) case(1:2) ! cubic(s) @@ -605,12 +608,10 @@ do i = 1,maxNinstance / ( 4.0_pReal*constitutive_nonlocal_C11(i) + 6.0_pReal*constitutive_nonlocal_C12(i) & + 2.0_pReal*constitutive_nonlocal_C44(i) ) ! C12iso/(C11iso+C12iso) with C11iso=(3*C11+2*C12+4*C44)/5 and C12iso=(C11+4*C12-2*C44)/5 - - do s1 = 1,ns - + do s1 = 1,ns f = constitutive_nonlocal_slipFamily(s1,i) -!*** burgers vector, mean free path prefactor and minimum dipole distance for each slip system + !*** burgers vector, mean free path prefactor and minimum dipole distance for each slip system constitutive_nonlocal_burgersPerSlipSystem(s1,i) = constitutive_nonlocal_burgersPerSlipFamily(f,i) constitutive_nonlocal_lambda0PerSlipSystem(s1,i) = constitutive_nonlocal_lambda0PerSlipFamily(f,i) @@ -619,31 +620,33 @@ do i = 1,maxNinstance do s2 = 1,ns -!*** calculation of forest projections for edge and screw dislocations. s2 acts as forest to s1 + !*** calculation of forest projections for edge and screw dislocations. s2 acts as forest for s1 - constitutive_nonlocal_forestProjectionEdge(s1, s2, i) & - = abs(math_mul3x3(lattice_sn(1:3, constitutive_nonlocal_slipSystemLattice(s1,i), myStructure), & - lattice_st(1:3, constitutive_nonlocal_slipSystemLattice(s2,i), myStructure))) ! forest projection of edge dislocations is the projection of (t = b x n) onto the slip normal of the respective slip plane + constitutive_nonlocal_forestProjectionEdge(s1,s2,i) & + = abs(math_mul3x3(lattice_sn(1:3,constitutive_nonlocal_slipSystemLattice(s1,i),myStructure), & + lattice_st(1:3,constitutive_nonlocal_slipSystemLattice(s2,i),myStructure))) ! forest projection of edge dislocations is the projection of (t = b x n) onto the slip normal of the respective slip plane - constitutive_nonlocal_forestProjectionScrew(s1, s2, i) & - = abs(math_mul3x3(lattice_sn(1:3, constitutive_nonlocal_slipSystemLattice(s1,i), myStructure), & - lattice_sd(1:3, constitutive_nonlocal_slipSystemLattice(s2,i), myStructure))) ! forest projection of screw dislocations is the projection of b onto the slip normal of the respective splip plane + constitutive_nonlocal_forestProjectionScrew(s1,s2,i) & + = abs(math_mul3x3(lattice_sn(1:3,constitutive_nonlocal_slipSystemLattice(s1,i),myStructure), & + lattice_sd(1:3,constitutive_nonlocal_slipSystemLattice(s2,i),myStructure))) ! forest projection of screw dislocations is the projection of b onto the slip normal of the respective splip plane -!*** calculation of interaction matrices + !*** calculation of interaction matrices - constitutive_nonlocal_interactionMatrixSlipSlip(s1, s2, i) & - = constitutive_nonlocal_interactionSlipSlip( lattice_interactionSlipSlip(constitutive_nonlocal_slipSystemLattice(s1,i), & - constitutive_nonlocal_slipSystemLattice(s2,i), & - myStructure), & - i ) + constitutive_nonlocal_interactionMatrixSlipSlip(s1,s2,i) & + = constitutive_nonlocal_interactionSlipSlip(lattice_interactionSlipSlip(constitutive_nonlocal_slipSystemLattice(s1,i), & + constitutive_nonlocal_slipSystemLattice(s2,i), & + myStructure), & + i) - enddo; enddo + enddo + enddo -!*** calculation of prefactor for activation enthalpy for dislocation glide + + !*** calculation of prefactor for activation enthalpy for dislocation glide constitutive_nonlocal_Qeff0(1:ns,i) = constitutive_nonlocal_burgersPerSlipSystem(1:ns,i) ** 3.0_pReal & - * sqrt(0.5_pReal * constitutive_nonlocal_d0(i) ** 3.0_pReal & - * constitutive_nonlocal_Gmod(i) * constitutive_nonlocal_tauObs(i)) + * sqrt(0.5_pReal * constitutive_nonlocal_d0(i) ** 3.0_pReal & + * constitutive_nonlocal_Gmod(i) * constitutive_nonlocal_tauObs(i)) enddo @@ -1029,16 +1032,6 @@ if (.not. phase_localConstitution(myPhase)) then endif -if (debug_verbosity > 6 .and. ((debug_e == el .and. debug_i == ip .and. debug_g == g) .or. .not. debug_selectiveDebugger)) then - !$OMP CRITICAL (write2out) - write(6,*) - write(6,'(a,i5,x,i2,x,i1)') '<< CONST >> nonlocal_microstructure at el ip g',el,ip,g - write(6,*) - write(6,'(a,/,12(x),12(e10.3,x))') '<< CONST >> rhoForest', rhoForest - write(6,'(a,/,12(x),12(f10.5,x))') '<< CONST >> tauThreshold / MPa', tauThreshold/1e6 - write(6,'(a,/,3(12(x),3(f10.5,x),/))') '<< CONST >> Tdislocation / MPa', math_Mandel6to33(Tdislocation_v)/1e6 - !$OMP END CRITICAL (write2out) -endif !********************************************************************** !*** set states @@ -1049,6 +1042,18 @@ state(g,ip,el)%p(10*ns+1:11*ns) = rhoForest state(g,ip,el)%p(11*ns+1:12*ns) = tauThreshold state(g,ip,el)%p(12*ns+1:12*ns+6) = Tdislocation_v + +#ifndef _OPENMP + if (debug_verbosity > 6 .and. ((debug_e == el .and. debug_i == ip .and. debug_g == g) .or. .not. debug_selectiveDebugger)) then + write(6,*) + write(6,'(a,i5,x,i2,x,i1)') '<< CONST >> nonlocal_microstructure at el ip g',el,ip,g + write(6,*) + write(6,'(a,/,12(x),12(e10.3,x))') '<< CONST >> rhoForest', rhoForest + write(6,'(a,/,12(x),12(f10.5,x))') '<< CONST >> tauThreshold / MPa', tauThreshold/1e6 + write(6,'(a,/,3(12(x),3(f10.5,x),/))') '<< CONST >> Tdislocation / MPa', math_Mandel6to33(Tdislocation_v)/1e6 + endif +#endif + endsubroutine @@ -1166,15 +1171,15 @@ endif constitutive_nonlocal_v(1:ns,1:4,g,ip,el) = v !$OMP FLUSH(constitutive_nonlocal_v) -if (debug_verbosity > 6 .and. ((debug_e == el .and. debug_i == ip .and. debug_g == g) .or. .not. debug_selectiveDebugger)) then - !$OMP CRITICAL (write2out) +#ifndef _OPENMP + if (debug_verbosity > 6 .and. ((debug_e == el .and. debug_i == ip .and. debug_g == g) .or. .not. debug_selectiveDebugger)) then write(6,*) write(6,'(a,i5,x,i2,x,i1)') '<< CONST >> nonlocal_kinetics at el ip g',el,ip,g write(6,*) write(6,'(a,/,12(x),12(f12.5,x))') '<< CONST >> tau / MPa', tau/1e6_pReal write(6,'(a,/,4(12(x),12(f12.5,x),/))') '<< CONST >> v / 1e-3m/s', constitutive_nonlocal_v(:,:,g,ip,el)*1e3 - !$OMP END CRITICAL (write2out) -endif + endif +#endif endsubroutine @@ -1289,16 +1294,16 @@ do s = 1,ns enddo dLp_dTstar99 = math_Plain3333to99(dLp_dTstar3333) -if (debug_verbosity > 6 .and. ((debug_e == el .and. debug_i == ip .and. debug_g == g) .or. .not. debug_selectiveDebugger)) then - !$OMP CRITICAL (write2out) +#ifndef _OPENMP + if (debug_verbosity > 6 .and. ((debug_e == el .and. debug_i == ip .and. debug_g == g) .or. .not. debug_selectiveDebugger)) then write(6,*) write(6,'(a,i5,x,i2,x,i1)') '<< CONST >> nonlocal_LpandItsTangent at el ip g ',el,ip,g write(6,*) write(6,'(a,/,4(12(x),12(f12.5,x)),/)') '<< CONST >> gdot / 1e-3',gdot*1e3_pReal write(6,'(a,/,12(x),12(f12.5,x))') '<< CONST >> gdot total / 1e-3',gdotTotal*1e3_pReal write(6,'(a,/,3(12(x),3(f12.7,x),/))') '<< CONST >> Lp',Lp - !$OMP END CRITICAL (write2out) -endif + endif +#endif endsubroutine @@ -1441,21 +1446,19 @@ real(pReal) area, & ! area logical considerEnteringFlux, & considerLeavingFlux -if (debug_verbosity > 6 .and. ((debug_e == el .and. debug_i == ip .and. debug_g == g) .or. .not. debug_selectiveDebugger)) then - !$OMP CRITICAL (write2out) +#ifndef _OPENMP + if (debug_verbosity > 6 .and. ((debug_e == el .and. debug_i == ip .and. debug_g == g) .or. .not. debug_selectiveDebugger)) then write(6,*) write(6,'(a,i5,x,i2,x,i1)') '<< CONST >> nonlocal_dotState at el ip g ',el,ip,g write(6,*) - !$OMP END CRITICAL (write2out) -endif + endif +#endif select case(mesh_element(2,el)) case (1,6,7,8,9) ! all fine case default call IO_error(-1,el,ip,g,'element type not supported for nonlocal constitution') - - end select myInstance = phase_constitutionInstance(material_phase(g,ip,el)) @@ -1489,7 +1492,7 @@ endif !**************************************************************************** !*** Calculate shear rate -call constitutive_nonlocal_kinetics(Tstar_v, Temperature, state(g,ip,el), g, ip, el) ! get velocities +call constitutive_nonlocal_kinetics(Tstar_v, Temperature, state(g,ip,el), g, ip, el) ! velocities forall (t = 1:4) & gdot(1:ns,t) = rhoSgl(1:ns,t) * constitutive_nonlocal_burgersPerSlipSystem(1:ns,myInstance) & @@ -1498,12 +1501,12 @@ forall (s = 1:ns, t = 1:4, rhoSgl(s,t+4) * constitutive_nonlocal_v(s,t,g,ip,el) gdot(s,t) = gdot(s,t) + abs(rhoSgl(s,t+4)) * constitutive_nonlocal_burgersPerSlipSystem(s,myInstance) & * constitutive_nonlocal_v(s,t,g,ip,el) -if (debug_verbosity > 6 .and. ((debug_e == el .and. debug_i == ip .and. debug_g == g) .or. .not. debug_selectiveDebugger)) then - !$OMP CRITICAL (write2out) +#ifndef _OPENMP + if (debug_verbosity > 6 .and. ((debug_e == el .and. debug_i == ip .and. debug_g == g) .or. .not. debug_selectiveDebugger)) then write(6,'(a,/,10(12(x),12(e12.5,x),/))') '<< CONST >> rho / 1/m^2', rhoSgl, rhoDip write(6,'(a,/,4(12(x),12(e12.5,x),/))') '<< CONST >> gdot / 1/s',gdot - !$OMP END CRITICAL (write2out) -endif + endif +#endif @@ -1751,8 +1754,10 @@ forall (t = 1:10) & + rhoDotAthermalAnnihilation(1:ns,t) & + rhoDotThermalAnnihilation(1:ns,t) -if (debug_verbosity > 6 .and. ((debug_e == el .and. debug_i == ip .and. debug_g == g) .or. .not. debug_selectiveDebugger)) then - !$OMP CRITICAL (write2out) +dotState%p(1:10*ns) = dotState%p(1:10*ns) + reshape(rhoDot,(/10*ns/)) + +#ifndef _OPENMP + if (debug_verbosity > 6 .and. ((debug_e == el .and. debug_i == ip .and. debug_g == g) .or. .not. debug_selectiveDebugger)) then write(6,'(a,/,8(12(x),12(e12.5,x),/))') '<< CONST >> dislocation remobilization', rhoDotRemobilization(1:ns,1:8) * timestep write(6,'(a,/,4(12(x),12(e12.5,x),/))') '<< CONST >> dislocation multiplication', rhoDotMultiplication(1:ns,1:4) * timestep write(6,'(a,/,8(12(x),12(e12.5,x),/))') '<< CONST >> dislocation flux', rhoDotFlux(1:ns,1:8) * timestep @@ -1765,10 +1770,8 @@ if (debug_verbosity > 6 .and. ((debug_e == el .and. debug_i == ip .and. debug_g write(6,'(a,/,10(12(x),12(f12.7,x),/))') '<< CONST >> relative density change', & rhoDot(1:ns,1:8) * timestep / (abs(rhoSgl)+1.0e-10), & rhoDot(1:ns,9:10) * timestep / (rhoDip+1.0e-10) - !$OMP END CRITICAL (write2out) -endif - -dotState%p(1:10*ns) = dotState%p(1:10*ns) + reshape(rhoDot,(/10*ns/)) + endif +#endif endsubroutine @@ -2089,7 +2092,7 @@ forall (t = 1:4) & do s = 1,ns sLattice = constitutive_nonlocal_slipSystemLattice(s,myInstance) - tau(s) = math_mul6x6( Tstar_v + Tdislocation_v, lattice_Sslip_v(1:6,sLattice,myStructure) ) + tau(s) = math_mul6x6(Tstar_v + Tdislocation_v, lattice_Sslip_v(1:6,sLattice,myStructure)) enddo dLower(1:ns,1) = constitutive_nonlocal_dLowerEdgePerSlipSystem(1:ns,myInstance) diff --git a/code/crystallite.f90 b/code/crystallite.f90 index 71d3259d5..0447bd529 100644 --- a/code/crystallite.f90 +++ b/code/crystallite.f90 @@ -15,7 +15,7 @@ MODULE crystallite use prec, only: pReal, pInt implicit none -! + ! **************************************************************** ! *** General variables for the crystallite calculation *** ! **************************************************************** @@ -76,6 +76,7 @@ logical, dimension (:,:,:), allocatable :: & CONTAINS + !******************************************************************** ! allocate and initialize per grain variables !******************************************************************** @@ -213,7 +214,7 @@ allocate(crystallite_output(maxval(crystallite_Noutput), & material_Ncrystallite)) ; crystallite_output = '' allocate(crystallite_sizePostResults(material_Ncrystallite)) ; crystallite_sizePostResults = 0_pInt allocate(crystallite_sizePostResult(maxval(crystallite_Noutput), & - material_Ncrystallite)) ; crystallite_sizePostResult = 0_pInt + material_Ncrystallite)) ; crystallite_sizePostResult = 0_pInt if(.not. IO_open_file(file,material_configFile)) call IO_error (100) ! corrupt config file @@ -450,17 +451,13 @@ use math, only: math_inv3x3, & math_transpose3x3, & math_I3 use FEsolving, only: FEsolving_execElem, & - FEsolving_execIP, & - theInc, & - cycleCounter + FEsolving_execIP use mesh, only: mesh_element, & mesh_NcpElems, & mesh_maxNips use material, only: homogenization_Ngrains, & homogenization_maxNgrains -use constitutive, only: constitutive_maxSizeState, & - constitutive_maxSizeDotState, & - constitutive_sizeState, & +use constitutive, only: constitutive_sizeState, & constitutive_sizeDotState, & constitutive_state, & constitutive_state_backup, & @@ -468,10 +465,7 @@ use constitutive, only: constitutive_maxSizeState, constitutive_partionedState0, & constitutive_homogenizedC, & constitutive_dotState, & - constitutive_dotState_backup, & - constitutive_collectDotState, & - constitutive_dotTemperature, & - constitutive_microstructure + constitutive_dotState_backup implicit none @@ -517,21 +511,21 @@ logical, dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems) :: & ! --+>> INITIALIZE TO STARTING CONDITION <<+-- - if (debug_verbosity > 4 .and. debug_e > 0 .and. debug_e <= mesh_NcpElems & - .and. debug_i > 0 .and. debug_i <= mesh_maxNips & - .and. debug_g > 0 .and. debug_g <= homogenization_maxNgrains) then - !$OMP CRITICAL (write2out) - write (6,*) - write (6,'(a,i5,x,i2,x,i3)') '<< CRYST >> crystallite start at el ip g ', debug_e, debug_i, debug_g - write (6,'(a,/,12(x),f14.9)') '<< CRYST >> Temp0', crystallite_partionedTemperature0(debug_g,debug_i,debug_e) - write (6,'(a,/,3(12(x),3(f14.9,x)/))') '<< CRYST >> F0 ', & - math_transpose3x3(crystallite_partionedF0(1:3,1:3,debug_g,debug_i,debug_e)) - write (6,'(a,/,3(12(x),3(f14.9,x)/))') '<< CRYST >> Fp0', & - math_transpose3x3(crystallite_partionedFp0(1:3,1:3,debug_g,debug_i,debug_e)) - write (6,'(a,/,3(12(x),3(f14.9,x)/))') '<< CRYST >> Lp0', & - math_transpose3x3(crystallite_partionedLp0(1:3,1:3,debug_g,debug_i,debug_e)) - !$OMP END CRITICAL (write2out) - endif +if (debug_verbosity > 4 .and. debug_e > 0 .and. debug_e <= mesh_NcpElems & + .and. debug_i > 0 .and. debug_i <= mesh_maxNips & + .and. debug_g > 0 .and. debug_g <= homogenization_maxNgrains) then + !$OMP CRITICAL (write2out) + write (6,*) + write (6,'(a,i5,x,i2,x,i3)') '<< CRYST >> crystallite start at el ip g ', debug_e, debug_i, debug_g + write (6,'(a,/,12(x),f14.9)') '<< CRYST >> Temp0', crystallite_partionedTemperature0(debug_g,debug_i,debug_e) + write (6,'(a,/,3(12(x),3(f14.9,x)/))') '<< CRYST >> F0 ', & + math_transpose3x3(crystallite_partionedF0(1:3,1:3,debug_g,debug_i,debug_e)) + write (6,'(a,/,3(12(x),3(f14.9,x)/))') '<< CRYST >> Fp0', & + math_transpose3x3(crystallite_partionedFp0(1:3,1:3,debug_g,debug_i,debug_e)) + write (6,'(a,/,3(12(x),3(f14.9,x)/))') '<< CRYST >> Lp0', & + math_transpose3x3(crystallite_partionedLp0(1:3,1:3,debug_g,debug_i,debug_e)) + !$OMP END CRITICAL (write2out) +endif crystallite_subStep = 0.0_pReal @@ -575,15 +569,15 @@ do while (any(crystallite_subStep(:,:,FEsolving_execELem(1):FEsolving_execElem(2 ! --- wind forward --- if (crystallite_converged(g,i,e)) then +#ifndef _OPENMP if (debug_verbosity > 4 & .and. ((e == debug_e .and. i == debug_i .and. g == debug_g) .or. .not. debug_selectiveDebugger)) then - !$OMP CRITICAL (write2out) - write(6,'(a,f10.8,a,f10.8,a)') '<< CRYST >> winding forward from ', & - crystallite_subFrac(g,i,e),' to current crystallite_subfrac ', & - crystallite_subFrac(g,i,e)+crystallite_subStep(g,i,e),' in crystallite_stressAndItsTangent' - write(6,*) - !$OMP END CRITICAL (write2out) + write(6,'(a,f10.8,a,f10.8,a)') '<< CRYST >> winding forward from ', & + crystallite_subFrac(g,i,e),' to current crystallite_subfrac ', & + crystallite_subFrac(g,i,e)+crystallite_subStep(g,i,e),' in crystallite_stressAndItsTangent' + write(6,*) endif +#endif crystallite_subFrac(g,i,e) = crystallite_subFrac(g,i,e) + crystallite_subStep(g,i,e) formerSubStep = crystallite_subStep(g,i,e) !$OMP FLUSH(crystallite_subFrac) @@ -618,15 +612,15 @@ do while (any(crystallite_subStep(:,:,FEsolving_execELem(1):FEsolving_execElem(2 constitutive_state(g,i,e)%p = constitutive_subState0(g,i,e)%p ! ...microstructure crystallite_Tstar_v(1:6,g,i,e) = crystallite_subTstar0_v(1:6,g,i,e) ! ...2nd PK stress ! cant restore dotState here, since not yet calculated in first cutback after initialization - !$OMP FLUSH(crystallite_subStep,crystallite_invFp) + !$OMP FLUSH(crystallite_invFp) +#ifndef _OPENMP if (debug_verbosity > 4 & .and. ((e == debug_e .and. i == debug_i .and. g == debug_g) .or. .not. debug_selectiveDebugger)) then - !$OMP CRITICAL (write2out) - write(6,'(a,f10.8)') '<< CRYST >> cutback step in crystallite_stressAndItsTangent with new crystallite_subStep: ',& - crystallite_subStep(g,i,e) - write(6,*) - !$OMP END CRITICAL (write2out) + write(6,'(a,f10.8)') '<< CRYST >> cutback step in crystallite_stressAndItsTangent with new crystallite_subStep: ',& + crystallite_subStep(g,i,e) + write(6,*) endif +#endif endif ! --- prepare for integration --- @@ -684,18 +678,17 @@ enddo math_Mandel33to6( math_mul33x33(transpose(Fe_guess),Fe_guess) - math_I3 ) ) ) crystallite_P(1:3,1:3,g,i,e) = math_mul33x33(Fe_guess,math_mul33x33(Tstar,transpose(invFp))) endif +#ifndef _OPENMP if (debug_verbosity > 4 & .and. ((e == debug_e .and. i == debug_i .and. g == debug_g) .or. .not. debug_selectiveDebugger)) then - !$OMP FLUSH(crystallite_P) - !$OMP CRITICAL (write2out) - write (6,'(a,i5,x,i2,x,i3)') '<< CRYST >> central solution of cryst_StressAndTangent at el ip g ',e,i,g - write (6,*) - write (6,'(a,/,3(12(x),3(f12.4,x)/))') '<< CRYST >> P / MPa', math_transpose3x3(crystallite_P(1:3,1:3,g,i,e)) / 1e6 - write (6,'(a,/,3(12(x),3(f14.9,x)/))') '<< CRYST >> Fp', math_transpose3x3(crystallite_Fp(1:3,1:3,g,i,e)) - write (6,'(a,/,3(12(x),3(f14.9,x)/))') '<< CRYST >> Lp', math_transpose3x3(crystallite_Lp(1:3,1:3,g,i,e)) - write (6,*) - !$OMP END CRITICAL (write2out) + write (6,'(a,i5,x,i2,x,i3)') '<< CRYST >> central solution of cryst_StressAndTangent at el ip g ',e,i,g + write (6,*) + write (6,'(a,/,3(12(x),3(f12.4,x)/))') '<< CRYST >> P / MPa', math_transpose3x3(crystallite_P(1:3,1:3,g,i,e)) / 1e6 + write (6,'(a,/,3(12(x),3(f14.9,x)/))') '<< CRYST >> Fp', math_transpose3x3(crystallite_Fp(1:3,1:3,g,i,e)) + write (6,'(a,/,3(12(x),3(f14.9,x)/))') '<< CRYST >> Lp', math_transpose3x3(crystallite_Lp(1:3,1:3,g,i,e)) + write (6,*) endif +#endif enddo enddo enddo @@ -1007,18 +1000,18 @@ do n = 1,4 if (crystallite_todo(g,i,e)) then if (crystallite_integrateStress(g,i,e,timeStepFraction(n))) then ! fraction of original times step if (n == 4) then ! final integration step +#ifndef _OPENMP if (debug_verbosity > 5 & .and. ((e == debug_e .and. i == debug_i .and. g == debug_g) .or. .not. debug_selectiveDebugger)) then mySizeDotState = constitutive_sizeDotState(g,i,e) - !$OMP CRITICAL (write2out) - write(6,'(a,i5,x,i2,x,i3)') '<< CRYST >> updateState at el ip g ',e,i,g - write(6,*) - write(6,'(a,/,(12(x),12(e12.5,x)))') '<< CRYST >> dotState', constitutive_dotState(g,i,e)%p(1:mySizeDotState) - write(6,*) - write(6,'(a,/,(12(x),12(e12.5,x)))') '<< CRYST >> new state', constitutive_state(g,i,e)%p(1:mySizeDotState) - write(6,*) - !$OMP END CRITICAL (write2out) + write(6,'(a,i5,x,i2,x,i3)') '<< CRYST >> updateState at el ip g ',e,i,g + write(6,*) + write(6,'(a,/,(12(x),12(e12.5,x)))') '<< CRYST >> dotState', constitutive_dotState(g,i,e)%p(1:mySizeDotState) + write(6,*) + write(6,'(a,/,(12(x),12(e12.5,x)))') '<< CRYST >> new state', constitutive_state(g,i,e)%p(1:mySizeDotState) + write(6,*) endif +#endif crystallite_converged(g,i,e) = .true. ! ... converged per definition crystallite_todo(g,i,e) = .false. ! ... integration done if (debug_verbosity > 0) then @@ -1238,13 +1231,11 @@ endif ! --- FIRST RUNGE KUTTA STEP --- +#ifndef _OPENMP if (debug_verbosity > 5) then - !$OMP SINGLE - !$OMP CRITICAL (write2out) - write(6,'(a,x,i1)') '<< CRYST >> RUNGE KUTTA STEP',1 - !$OMP END CRITICAL (write2out) - !$OMP END SINGLE + write(6,'(a,x,i1)') '<< CRYST >> RUNGE KUTTA STEP',1 endif +#endif !$OMP DO do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains if (crystallite_todo(g,i,e)) then @@ -1253,7 +1244,7 @@ endif crystallite_dotTemperature(g,i,e) = constitutive_dotTemperature(crystallite_Tstar_v(1:6,g,i,e), & crystallite_Temperature(g,i,e),g,i,e) endif - enddo; enddo; enddo + enddo; enddo; enddo !$OMP ENDDO !$OMP DO do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains @@ -1376,13 +1367,11 @@ do n = 1,5 ! --- dot state and RK dot state--- +#ifndef _OPENMP if (debug_verbosity > 5) then - !$OMP SINGLE - !$OMP CRITICAL (write2out) - write(6,'(a,x,i1)') '<< CRYST >> RUNGE KUTTA STEP',n+1 - !$OMP END CRITICAL (write2out) - !$OMP END SINGLE + write(6,'(a,x,i1)') '<< CRYST >> RUNGE KUTTA STEP',n+1 endif +#endif !$OMP DO do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains if (crystallite_todo(g,i,e)) then @@ -1502,24 +1491,23 @@ relTemperatureResiduum = 0.0_pReal .or. abs(stateResiduum(1:mySizeDotState,g,i,e)) < constitutive_aTolState(g,i,e)%p(1:mySizeDotState) ) & .and. abs(relTemperatureResiduum(g,i,e)) < rTol_crystalliteTemperature ) +#ifndef _OPENMP if (debug_verbosity > 5 & .and. ((e == debug_e .and. i == debug_i .and. g == debug_g) .or. .not. debug_selectiveDebugger)) then - !$OMP CRITICAL (write2out) - write(6,'(a,i5,x,i3,x,i3)') '<< CRYST >> updateState at el ip g ',e,i,g - write(6,*) - write(6,'(a,/,(12(x),12(f12.1,x)))') '<< CRYST >> absolute residuum tolerance', & - stateResiduum(1:mySizeDotState,g,i,e) / constitutive_aTolState(g,i,e)%p(1:mySizeDotState) - write(6,*) - write(6,'(a,/,(12(x),12(f12.1,x)))') '<< CRYST >> relative residuum tolerance', & - relStateResiduum(1:mySizeDotState,g,i,e) / rTol_crystalliteState - write(6,*) - write(6,'(a,/,(12(x),12(e12.5,x)))') '<< CRYST >> dotState', constitutive_dotState(g,i,e)%p(1:mySizeDotState) - write(6,*) - write(6,'(a,/,(12(x),12(e12.5,x)))') '<< CRYST >> new state', constitutive_state(g,i,e)%p(1:mySizeDotState) - write(6,*) - !$OMP END CRITICAL (write2out) + write(6,'(a,i5,x,i3,x,i3)') '<< CRYST >> updateState at el ip g ',e,i,g + write(6,*) + write(6,'(a,/,(12(x),12(f12.1,x)))') '<< CRYST >> absolute residuum tolerance', & + stateResiduum(1:mySizeDotState,g,i,e) / constitutive_aTolState(g,i,e)%p(1:mySizeDotState) + write(6,*) + write(6,'(a,/,(12(x),12(f12.1,x)))') '<< CRYST >> relative residuum tolerance', & + relStateResiduum(1:mySizeDotState,g,i,e) / rTol_crystalliteState + write(6,*) + write(6,'(a,/,(12(x),12(e12.5,x)))') '<< CRYST >> dotState', constitutive_dotState(g,i,e)%p(1:mySizeDotState) + write(6,*) + write(6,'(a,/,(12(x),12(e12.5,x)))') '<< CRYST >> new state', constitutive_state(g,i,e)%p(1:mySizeDotState) + write(6,*) endif - +#endif endif enddo; enddo; enddo !$OMP ENDDO @@ -1565,11 +1553,11 @@ relTemperatureResiduum = 0.0_pReal ! --- nonlocal convergence check --- +#ifndef _OPENMP if (debug_verbosity > 5) then - !$OMP CRITICAL (write2out) - write(6,'(a,L)') '<< CRYST >> crystallite_converged',crystallite_converged - !$OMP END CRITICAL (write2out) + write(6,'(a,L)') '<< CRYST >> crystallite_converged',crystallite_converged endif +#endif if (.not. singleRun) then ! if not requesting Integration of just a single IP if ( any(.not. crystallite_converged .and. .not. crystallite_localConstitution)) then ! any non-local not yet converged (or broken)... crystallite_converged = crystallite_converged .and. crystallite_localConstitution ! ...restart all non-local as not converged @@ -1811,26 +1799,25 @@ relTemperatureResiduum = 0.0_pReal if (crystallite_Temperature(g,i,e) > 0) & relTemperatureResiduum(g,i,e) = temperatureResiduum(g,i,e) / crystallite_Temperature(g,i,e) !$OMP FLUSH(relStateResiduum,relTemperatureResiduum) - + +#ifndef _OPENMP if (debug_verbosity > 5 & .and. ((e == debug_e .and. i == debug_i .and. g == debug_g) .or. .not. debug_selectiveDebugger)) then - !$OMP CRITICAL (write2out) - write(6,'(a,i5,x,i2,x,i3)') '<< CRYST >> updateState at el ip g ',e,i,g - write(6,*) - write(6,'(a,/,(12(x),12(f12.1,x)))') '<< CRYST >> absolute residuum tolerance', & - stateResiduum(1:mySizeDotState,g,i,e) / constitutive_aTolState(g,i,e)%p(1:mySizeDotState) - write(6,*) - write(6,'(a,/,(12(x),12(f12.1,x)))') '<< CRYST >> relative residuum tolerance', & - relStateResiduum(1:mySizeDotState,g,i,e) / rTol_crystalliteState - write(6,*) - write(6,'(a,/,(12(x),12(e12.5,x)))') '<< CRYST >> dotState', constitutive_dotState(g,i,e)%p(1:mySizeDotState) & - - 2.0_pReal * stateResiduum(1:mySizeDotState,g,i,e) / crystallite_subdt(g,i,e) ! calculate former dotstate from higher order solution and state residuum - write(6,*) - write(6,'(a,/,(12(x),12(e12.5,x)))') '<< CRYST >> new state', constitutive_state(g,i,e)%p(1:mySizeDotState) - write(6,*) - !$OMP END CRITICAL (write2out) + write(6,'(a,i5,x,i2,x,i3)') '<< CRYST >> updateState at el ip g ',e,i,g + write(6,*) + write(6,'(a,/,(12(x),12(f12.1,x)))') '<< CRYST >> absolute residuum tolerance', & + stateResiduum(1:mySizeDotState,g,i,e) / constitutive_aTolState(g,i,e)%p(1:mySizeDotState) + write(6,*) + write(6,'(a,/,(12(x),12(f12.1,x)))') '<< CRYST >> relative residuum tolerance', & + relStateResiduum(1:mySizeDotState,g,i,e) / rTol_crystalliteState + write(6,*) + write(6,'(a,/,(12(x),12(e12.5,x)))') '<< CRYST >> dotState', constitutive_dotState(g,i,e)%p(1:mySizeDotState) & + - 2.0_pReal * stateResiduum(1:mySizeDotState,g,i,e) / crystallite_subdt(g,i,e) ! calculate former dotstate from higher order solution and state residuum + write(6,*) + write(6,'(a,/,(12(x),12(e12.5,x)))') '<< CRYST >> new state', constitutive_state(g,i,e)%p(1:mySizeDotState) + write(6,*) endif - +#endif ! --- converged ? --- @@ -1854,11 +1841,11 @@ relTemperatureResiduum = 0.0_pReal ! --- NONLOCAL CONVERGENCE CHECK --- +#ifndef _OPENMP if (debug_verbosity > 5) then - !$OMP CRITICAL (write2out) - write(6,'(a,L)') '<< CRYST >> crystallite_converged',crystallite_converged - !$OMP END CRITICAL (write2out) + write(6,'(a,L)') '<< CRYST >> crystallite_converged',crystallite_converged endif +#endif if (.not. singleRun) then ! if not requesting Integration of just a single IP if ( any(.not. crystallite_converged .and. .not. crystallite_localConstitution)) then ! any non-local not yet converged (or broken)... crystallite_converged = crystallite_converged .and. crystallite_localConstitution ! ...restart all non-local as not converged @@ -1988,16 +1975,16 @@ endif endif enddo; enddo; enddo !$OMP ENDDO +#ifndef _OPENMP if (debug_verbosity > 5 .and. ((e == debug_e .and. i == debug_i .and. g == debug_g) .or. .not. debug_selectiveDebugger)) then - !$OMP CRITICAL (write2out) - write(6,'(a,i5,x,i2,x,i3)') '<< CRYST >> updateState at el ip g ',e,i,g - write(6,*) - write(6,'(a,/,(12(x),12(e12.5,x)))') '<< CRYST >> dotState', constitutive_dotState(debug_g,debug_i,debug_e)%p(1:mySizeDotState) - write(6,*) - write(6,'(a,/,(12(x),12(e12.5,x)))') '<< CRYST >> new state', constitutive_state(debug_g,debug_i,debug_e)%p(1:mySizeDotState) - write(6,*) - !$OMP END CRITICAL (write2out) + write(6,'(a,i5,x,i2,x,i3)') '<< CRYST >> updateState at el ip g ',e,i,g + write(6,*) + write(6,'(a,/,(12(x),12(e12.5,x)))') '<< CRYST >> dotState', constitutive_dotState(debug_g,debug_i,debug_e)%p(1:mySizeDotState) + write(6,*) + write(6,'(a,/,(12(x),12(e12.5,x)))') '<< CRYST >> new state', constitutive_state(debug_g,debug_i,debug_e)%p(1:mySizeDotState) + write(6,*) endif +#endif ! --- UPDATE DEPENDENT STATES --- @@ -2104,7 +2091,9 @@ real(pReal) dot_prod12, & dot_prod22 logical singleRun, & ! flag indicating computation for single (g,i,e) triple stateConverged, & ! flag indicating convergence of state integration - temperatureConverged ! flag indicating convergence of temperature integration + temperatureConverged, & ! flag indicating convergence of temperature integration + stateUpdateDone, & ! flag indicating successfull state update + temperatureUpdateDone ! flag indicating successfull temperature update if (present(ee) .and. present(ii) .and. present(gg)) then @@ -2126,7 +2115,7 @@ endif ! --- RESET DOTSTATE --- -!$OMP PARALLEL PRIVATE(stateConverged,temperatureConverged) +!$OMP PARALLEL !$OMP DO do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains @@ -2154,14 +2143,14 @@ endif !$OMP SINGLE crystallite_statedamper = 1.0_pReal !$OMP END SINGLE -!$OMP DO +!$OMP DO PRIVATE(stateUpdateDone,temperatureUpdateDone,stateConverged,temperatureConverged) do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains if (crystallite_todo(g,i,e)) then - call crystallite_updateState(crystallite_todo(g,i,e), stateConverged, g,i,e) ! update state - !$OMP FLUSH(crystallite_todo) - call crystallite_updateTemperature(crystallite_todo(g,i,e), temperatureConverged, g,i,e) ! update temperature - !$OMP FLUSH(crystallite_todo) - if ( .not. crystallite_localConstitution(g,i,e) .and. .not. crystallite_todo(g,i,e) ) then ! if updateState or updateTemperature signals broken non-local... + call crystallite_updateState(stateUpdateDone, stateConverged, g,i,e) ! update state + call crystallite_updateTemperature(temperatureUpdateDone, temperatureConverged, g,i,e) ! update temperature + crystallite_todo(g,i,e) = stateUpdateDone .and. temperatureUpdateDone + if ( (.not. stateUpdateDone .or. .not. temperatureUpdateDone) & + .and. .not. crystallite_localConstitution(g,i,e) ) then ! if updateState or updateTemperature signals broken non-local... !$OMP CRITICAL (checkTodo) crystallite_todo = crystallite_todo .and. crystallite_localConstitution ! ...all non-locals skipped !$OMP END CRITICAL (checkTodo) @@ -2184,16 +2173,17 @@ endif constitutive_dotState(g,i,e)%p = 0.0_pReal ! reset dotState to zero enddo; enddo; enddo !$OMP ENDDO - !$OMP END PARALLEL + ! --+>> STATE LOOP <<+-- NiterationState = 0_pInt do while (any(crystallite_todo) .and. NiterationState < nState ) ! convergence loop for crystallite NiterationState = NiterationState + 1_pInt - !$OMP PARALLEL PRIVATE(stateConverged,temperatureConverged) + !$OMP PARALLEL + ! --- STRESS INTEGRATION --- @@ -2213,13 +2203,11 @@ do while (any(crystallite_todo) .and. NiterationState < nState ) enddo; enddo; enddo !$OMP ENDDO +#ifndef _OPENMP if (debug_verbosity > 5) then - !$OMP SINGLE - !$OMP CRITICAL (write2out) - write(6,'(a,i8,a)') '<< CRYST >> ', count(crystallite_todo(:,:,:)),' grains todo after stress integration' - !$OMP END CRITICAL (write2out) - !$OMP END SINGLE + write(6,'(a,i8,a)') '<< CRYST >> ', count(crystallite_todo(:,:,:)),' grains todo after stress integration' endif +#endif ! --- DOT STATES --- @@ -2239,7 +2227,7 @@ do while (any(crystallite_todo) .and. NiterationState < nState ) !$OMP SINGLE crystallite_statedamper = 1.0_pReal !$OMP END SINGLE - !$OMP DO PRIVATE(dot_prod12,dot_prod22) + !$OMP DO PRIVATE(dot_prod12,dot_prod22,stateUpdateDone,temperatureUpdateDone,stateConverged,temperatureConverged) do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains if (crystallite_todo(g,i,e)) then @@ -2258,12 +2246,12 @@ do while (any(crystallite_todo) .and. NiterationState < nState ) ! --- updates --- - call crystallite_updateState(crystallite_todo(g,i,e), stateConverged, g,i,e) ! update state - !$OMP FLUSH(crystallite_todo) - call crystallite_updateTemperature(crystallite_todo(g,i,e), temperatureConverged, g,i,e) ! update temperature - !$OMP FLUSH(crystallite_todo) + call crystallite_updateState(stateUpdateDone, stateConverged, g,i,e) ! update state + call crystallite_updateTemperature(temperatureUpdateDone, temperatureConverged, g,i,e) ! update temperature + crystallite_todo(g,i,e) = stateUpdateDone .and. temperatureUpdateDone crystallite_converged(g,i,e) = stateConverged .and. temperatureConverged - if ( .not. crystallite_localConstitution(g,i,e) .and. .not. crystallite_todo(g,i,e)) then ! if updateState or updateTemperature signals broken non-local... + if ( (.not. stateUpdateDone .or. .not. temperatureUpdateDone) & + .and. .not. crystallite_localConstitution(g,i,e) ) then ! if updateState or updateTemperature signals broken non-local... !$OMP CRITICAL (checkTodo) crystallite_todo = crystallite_todo .and. crystallite_localConstitution ! ...all non-locals skipped !$OMP END CRITICAL (checkTodo) @@ -2295,14 +2283,14 @@ do while (any(crystallite_todo) .and. NiterationState < nState ) !$OMP ENDDO !$OMP END PARALLEL - + +#ifndef _OPENMP if (debug_verbosity > 5) then - !$OMP CRITICAL (write2out) - write(6,'(a,i8,a,i2)') '<< CRYST >> ', count(crystallite_converged(:,:,:)), & - ' grains converged after state integration no. ', NiterationState - write(6,*) - !$OMP END CRITICAL (write2out) + write(6,'(a,i8,a,i2)') '<< CRYST >> ', count(crystallite_converged(:,:,:)), & + ' grains converged after state integration no. ', NiterationState + write(6,*) endif +#endif ! --- CONVERGENCE CHECK --- @@ -2314,14 +2302,14 @@ do while (any(crystallite_todo) .and. NiterationState < nState ) endif crystallite_todo = crystallite_todo .and. .not. crystallite_converged ! skip all converged +#ifndef _OPENMP if (debug_verbosity > 5) then - !$OMP CRITICAL (write2out) - write(6,'(a,i8,a)') '<< CRYST >> ', count(crystallite_converged(:,:,:)),' grains converged after non-local check' - write(6,'(a,i8,a,i2)') '<< CRYST >> ', count(crystallite_todo(:,:,:)),' grains todo after state integration no. ',& - NiterationState - write(6,*) - !$OMP END CRITICAL (write2out) + write(6,'(a,i8,a)') '<< CRYST >> ', count(crystallite_converged(:,:,:)),' grains converged after non-local check' + write(6,'(a,i8,a,i2)') '<< CRYST >> ', count(crystallite_todo(:,:,:)),' grains todo after state integration no. ',& + NiterationState + write(6,*) endif +#endif enddo ! crystallite convergence loop @@ -2388,11 +2376,11 @@ dotState(1:mySize) = constitutive_dotState(g,i,e)%p(1:mySize) * crystallite_stat residuum = constitutive_state(g,i,e)%p(1:mySize) - constitutive_subState0(g,i,e)%p(1:mySize) & - dotState(1:mySize) * crystallite_subdt(g,i,e) if (any(residuum /= residuum)) then ! if NaN occured then return without changing the state +#ifndef _OPENMP if (debug_verbosity > 4) then - !$OMP CRITICAL (write2out) - write(6,'(a,i5,x,i2,x,i3)') '<< CRYST >> updateState encountered NaN at el ip g ',e,i,g - !$OMP END CRITICAL (write2out) - endif + write(6,'(a,i5,x,i2,x,i3)') '<< CRYST >> updateState encountered NaN at el ip g ',e,i,g + endif +#endif return endif @@ -2404,22 +2392,22 @@ done = .true. converged = all( abs(residuum) < constitutive_aTolState(g,i,e)%p(1:mySize) & .or. abs(residuum) < rTol_crystalliteState * abs(state(1:mySize)) ) +#ifndef _OPENMP if (debug_verbosity > 5 .and. ((e == debug_e .and. i == debug_i .and. g == debug_g) .or. .not. debug_selectiveDebugger)) then - !$OMP CRITICAL (write2out) - if (converged) then - write(6,'(a,i5,x,i2,x,i3)') '<< CRYST >> updateState converged at el ip g ',e,i,g - else - write(6,'(a,i5,x,i2,x,i3)') '<< CRYST >> updateState did not converge at el ip g ',e,i,g - endif - write(6,*) - write(6,'(a,f6.1)') '<< CRYST >> crystallite_statedamper ',crystallite_statedamper(g,i,e) - write(6,*) - write(6,'(a,/,(12(x),12(e12.5,x)))') '<< CRYST >> dotState',dotState(1:mySize) - write(6,*) - write(6,'(a,/,(12(x),12(e12.5,x)))') '<< CRYST >> new state',state(1:mySize) - write(6,*) - !$OMP END CRITICAL (write2out) + if (converged) then + write(6,'(a,i5,x,i2,x,i3)') '<< CRYST >> updateState converged at el ip g ',e,i,g + else + write(6,'(a,i5,x,i2,x,i3)') '<< CRYST >> updateState did not converge at el ip g ',e,i,g + endif + write(6,*) + write(6,'(a,f6.1)') '<< CRYST >> crystallite_statedamper ',crystallite_statedamper(g,i,e) + write(6,*) + write(6,'(a,/,(12(x),12(e12.5,x)))') '<< CRYST >> dotState',dotState(1:mySize) + write(6,*) + write(6,'(a,/,(12(x),12(e12.5,x)))') '<< CRYST >> new state',state(1:mySize) + write(6,*) endif +#endif !* sync global state and dotState with local copies @@ -2471,11 +2459,11 @@ residuum = crystallite_Temperature(g,i,e) - crystallite_subTemperature0(g,i,e) & - constitutive_dotTemperature(crystallite_Tstar_v(1:6,g,i,e),crystallite_Temperature(g,i,e),g,i,e) & * crystallite_subdt(g,i,e) if (residuum /= residuum) then +#ifndef _OPENMP if (debug_verbosity > 4) then - !$OMP CRITICAL (write2out) - write(6,'(a,i5,x,i2,x,i3)') '<< CRYST >> updateTemperature encountered NaN at el ip g ',e,i,g - !$OMP END CRITICAL (write2out) + write(6,'(a,i5,x,i2,x,i3)') '<< CRYST >> updateTemperature encountered NaN at el ip g ',e,i,g endif +#endif return endif @@ -2598,11 +2586,11 @@ integer(pLongInt) tick, & !* be pessimistic crystallite_integrateStress = .false. +#ifndef _OPENMP if (debug_verbosity > 5 .and. ((e == debug_e .and. i == debug_i .and. g == debug_g) .or. .not. debug_selectiveDebugger)) then - !$OMP CRITICAL (write2out) - write(6,'(a,i5,x,i2,x,i3)') '<< CRYST >> integrateStress at el ip g ',e,i,g - !$OMP END CRITICAL (write2out) + write(6,'(a,i5,x,i2,x,i3)') '<< CRYST >> integrateStress at el ip g ',e,i,g endif +#endif !* only integrate over fraction of timestep? @@ -2628,15 +2616,15 @@ Lpguess = crystallite_Lp(1:3,1:3,g,i,e) ! ... and tak invFp_current = math_inv3x3(Fp_current) if (all(invFp_current == 0.0_pReal)) then ! ... failed? +#ifndef _OPENMP if (debug_verbosity > 4) then - !$OMP CRITICAL (write2out) - write(6,'(a,i5,x,i2,x,i3)') '<< CRYST >> integrateStress failed on invFp_current inversion at el ip g ',e,i,g - if (debug_verbosity > 5 .and. ((e == debug_e .and. i == debug_i .and. g == debug_g) .or. .not. debug_selectiveDebugger)) then - write(6,*) - write(6,'(a,/,3(12(x),3(f12.7,x)/))') '<< CRYST >> invFp_new',math_transpose3x3(invFp_new(1:3,1:3)) - endif - !$OMP END CRITICAL (write2out) + write(6,'(a,i5,x,i2,x,i3)') '<< CRYST >> integrateStress failed on invFp_current inversion at el ip g ',e,i,g + if (debug_verbosity > 5 .and. ((e == debug_e .and. i == debug_i .and. g == debug_g) .or. .not. debug_selectiveDebugger)) then + write(6,*) + write(6,'(a,/,3(12(x),3(f12.7,x)/))') '<< CRYST >> invFp_new',math_transpose3x3(invFp_new(1:3,1:3)) + endif endif +#endif return endif A = math_mul33x33(transpose(invFp_current), math_mul33x33(transpose(Fg_new),math_mul33x33(Fg_new,invFp_current))) @@ -2662,12 +2650,12 @@ LpLoop: do !* too many loops required ? if (NiterationStress > nStress) then +#ifndef _OPENMP if (debug_verbosity > 4) then - !$OMP CRITICAL (write2out) - write(6,'(a,i5,x,i2,x,i3)') '<< CRYST >> integrateStress reached loop limit at el ip g ',e,i,g - write(6,*) - !$OMP END CRITICAL (write2out) + write(6,'(a,i5,x,i2,x,i3)') '<< CRYST >> integrateStress reached loop limit at el ip g ',e,i,g + write(6,*) endif +#endif return endif @@ -2700,15 +2688,15 @@ LpLoop: do !$OMP END CRITICAL (debugTimingLpTangent) endif +#ifndef _OPENMP if (debug_verbosity > 5 .and. ((e == debug_e .and. i == debug_i .and. g == debug_g) .or. .not. debug_selectiveDebugger) & .and. numerics_integrationMode == 1_pInt) then - !$OMP CRITICAL (write2out) - write(6,'(a,i3)') '<< CRYST >> iteration ', NiterationStress - write(6,*) - write(6,'(a,/,3(12(x),3(e20.7,x)/))') '<< CRYST >> Lp_constitutive', math_transpose3x3(Lp_constitutive) - write(6,'(a,/,3(12(x),3(e20.7,x)/))') '<< CRYST >> Lpguess', math_transpose3x3(Lpguess) - !$OMP END CRITICAL (write2out) + write(6,'(a,i3)') '<< CRYST >> iteration ', NiterationStress + write(6,*) + write(6,'(a,/,3(12(x),3(e20.7,x)/))') '<< CRYST >> Lp_constitutive', math_transpose3x3(Lp_constitutive) + write(6,'(a,/,3(12(x),3(e20.7,x)/))') '<< CRYST >> Lpguess', math_transpose3x3(Lpguess) endif +#endif !* update current residuum and check for convergence of loop @@ -2728,13 +2716,13 @@ LpLoop: do !* NaN occured at regular speed? -> return if (any(residuum/=residuum) .and. leapfrog == 1.0) then +#ifndef _OPENMP if (debug_verbosity > 4) then - !$OMP CRITICAL (write2out) - write(6,'(a,i5,x,i2,x,i3,a,i3,a)') '<< CRYST >> integrateStress encountered NaN at el ip g ',e,i,g,& - ' ; iteration ', NiterationStress,& - ' >> returning..!' - !$OMP END CRITICAL (write2out) + write(6,'(a,i5,x,i2,x,i3,a,i3,a)') '<< CRYST >> integrateStress encountered NaN at el ip g ',e,i,g,& + ' ; iteration ', NiterationStress,& + ' >> returning..!' endif +#endif return @@ -2746,12 +2734,12 @@ LpLoop: do any(residuum/=residuum) & ! NaN occured ) & ) then +#ifndef _OPENMP if (debug_verbosity > 5) then - !$OMP CRITICAL (write2out) - write(6,'(a,i5,x,i2,x,i3,x,a,i3)') '<< CRYST >> integrateStress encountered high-speed crash at el ip g ',e,i,g,& - '; iteration ', NiterationStress - !$OMP END CRITICAL (write2out) + write(6,'(a,i5,x,i2,x,i3,x,a,i3)') '<< CRYST >> integrateStress encountered high-speed crash at el ip g ',e,i,g,& + '; iteration ', NiterationStress endif +#endif maxleap = 0.5_pReal * leapfrog ! limit next acceleration leapfrog = 1.0_pReal ! grinding halt jacoCounter = 0_pInt ! reset counter for Jacobian update (we want to do an update next time!) @@ -2778,21 +2766,21 @@ LpLoop: do invdRdLp = 0.0_pReal call math_invert(9,dRdLp,invdRdLp,dummy,error) ! invert dR/dLp --> dLp/dR if (error) then +#ifndef _OPENMP if (debug_verbosity > 4) then - !$OMP CRITICAL (write2out) - write(6,'(a,i5,x,i2,x,i3,a,i3)') '<< CRYST >> integrateStress failed on dR/dLp inversion at el ip g ',e,i,g,& - ' ; iteration ', NiterationStress - if (debug_verbosity > 5 & - .and. ((e == debug_e .and. i == debug_i .and. g == debug_g) .or. .not. debug_selectiveDebugger)) then - write(6,*) - write(6,'(a,/,9(12(x),9(e15.3,x)/))') '<< CRYST >> dRdLp',transpose(dRdLp) - write(6,'(a,/,9(12(x),9(e15.3,x)/))') '<< CRYST >> dLpdT_constitutive',transpose(dLpdT_constitutive) - write(6,'(a,/,3(12(x),3(e20.7,x)/))') '<< CRYST >> Lp_constitutive',math_transpose3x3(Lp_constitutive) - write(6,'(a,/,3(12(x),3(e20.7,x)/))') '<< CRYST >> Lpguess',math_transpose3x3(Lpguess) - endif - !$OMP END CRITICAL (write2out) - endif - return + write(6,'(a,i5,x,i2,x,i3,a,i3)') '<< CRYST >> integrateStress failed on dR/dLp inversion at el ip g ',e,i,g,& + ' ; iteration ', NiterationStress + if (debug_verbosity > 5 & + .and. ((e == debug_e .and. i == debug_i .and. g == debug_g) .or. .not. debug_selectiveDebugger)) then + write(6,*) + write(6,'(a,/,9(12(x),9(e15.3,x)/))') '<< CRYST >> dRdLp',transpose(dRdLp) + write(6,'(a,/,9(12(x),9(e15.3,x)/))') '<< CRYST >> dLpdT_constitutive',transpose(dLpdT_constitutive) + write(6,'(a,/,3(12(x),3(e20.7,x)/))') '<< CRYST >> Lp_constitutive',math_transpose3x3(Lp_constitutive) + write(6,'(a,/,3(12(x),3(e20.7,x)/))') '<< CRYST >> Lpguess',math_transpose3x3(Lpguess) + endif + endif +#endif + return endif endif jacoCounter = jacoCounter + 1_pInt ! increase counter for jaco update @@ -2822,16 +2810,16 @@ invFp_new = math_mul33x33(invFp_current,B) invFp_new = invFp_new/math_det3x3(invFp_new)**(1.0_pReal/3.0_pReal) ! regularize by det call math_invert3x3(invFp_new,Fp_new,det,error) if (error) then +#ifndef _OPENMP if (debug_verbosity > 4) then - !$OMP CRITICAL (write2out) - write(6,'(a,i5,x,i2,x,i3,a,i3)') '<< CRYST >> integrateStress failed on invFp_new inversion at el ip g ',e,i,g, & - ' ; iteration ', NiterationStress - if (debug_verbosity > 5 .and. ((e == debug_e .and. i == debug_i .and. g == debug_g) .or. .not. debug_selectiveDebugger)) then - write(6,*) - write(6,'(a,/,3(12(x),3(f12.7,x)/))') '<< CRYST >> invFp_new',math_transpose3x3(invFp_new) - endif - !$OMP END CRITICAL (write2out) + write(6,'(a,i5,x,i2,x,i3,a,i3)') '<< CRYST >> integrateStress failed on invFp_new inversion at el ip g ',e,i,g, & + ' ; iteration ', NiterationStress + if (debug_verbosity > 5 .and. ((e == debug_e .and. i == debug_i .and. g == debug_g) .or. .not. debug_selectiveDebugger)) then + write(6,*) + write(6,'(a,/,3(12(x),3(f12.7,x)/))') '<< CRYST >> invFp_new',math_transpose3x3(invFp_new) + endif endif +#endif return endif Fe_new = math_mul33x33(Fg_new,invFp_new) ! calc resulting Fe @@ -2855,17 +2843,17 @@ crystallite_invFp(1:3,1:3,g,i,e) = invFp_new !* set return flag to true crystallite_integrateStress = .true. +#ifndef _OPENMP if (debug_verbosity > 5 .and. ((e == debug_e .and. i == debug_i .and. g == debug_g) .or. .not. debug_selectiveDebugger) & .and. numerics_integrationMode == 1_pInt) then - !$OMP CRITICAL (write2out) - write(6,'(a,/,3(12(x),3(f12.7,x)/))') '<< CRYST >> P / MPa',math_transpose3x3(crystallite_P(1:3,1:3,g,i,e))/1e6 - write(6,'(a,/,3(12(x),3(f12.7,x)/))') '<< CRYST >> Cauchy / MPa', & - math_mul33x33(crystallite_P(1:3,1:3,g,i,e), math_transpose3x3(Fg_new)) / 1e6 / math_det3x3(Fg_new) - write(6,'(a,/,3(12(x),3(f12.7,x)/))') '<< CRYST >> Fe Lp Fe^-1', & - math_transpose3x3(math_mul33x33(Fe_new, math_mul33x33(crystallite_Lp(1:3,1:3,g,i,e), math_inv3x3(Fe_new)))) ! transpose to get correct print out order - write(6,'(a,/,3(12(x),3(f12.7,x)/))') '<< CRYST >> Fp',math_transpose3x3(crystallite_Fp(1:3,1:3,g,i,e)) - !$OMP END CRITICAL (write2out) + write(6,'(a,/,3(12(x),3(f12.7,x)/))') '<< CRYST >> P / MPa',math_transpose3x3(crystallite_P(1:3,1:3,g,i,e))/1e6 + write(6,'(a,/,3(12(x),3(f12.7,x)/))') '<< CRYST >> Cauchy / MPa', & + math_mul33x33(crystallite_P(1:3,1:3,g,i,e), math_transpose3x3(Fg_new)) / 1e6 / math_det3x3(Fg_new) + write(6,'(a,/,3(12(x),3(f12.7,x)/))') '<< CRYST >> Fe Lp Fe^-1', & + math_transpose3x3(math_mul33x33(Fe_new, math_mul33x33(crystallite_Lp(1:3,1:3,g,i,e), math_inv3x3(Fe_new)))) ! transpose to get correct print out order + write(6,'(a,/,3(12(x),3(f12.7,x)/))') '<< CRYST >> Fp',math_transpose3x3(crystallite_Fp(1:3,1:3,g,i,e)) endif +#endif if (debug_verbosity > 0) then !$OMP CRITICAL (distributionStress) diff --git a/code/homogenization.f90 b/code/homogenization.f90 index 98664e526..0636bdcc2 100644 --- a/code/homogenization.f90 +++ b/code/homogenization.f90 @@ -53,165 +53,169 @@ CONTAINS !* Module initialization * !************************************** subroutine homogenization_init(Temperature) - use prec, only: pReal,pInt - use math, only: math_I3 - use debug, only: debug_verbosity - use IO, only: IO_error, IO_open_file, IO_open_jobFile - use mesh, only: mesh_maxNips,mesh_NcpElems,mesh_element,FE_Nips - use material - use constitutive, only: constitutive_maxSizePostResults - use crystallite, only: crystallite_maxSizePostResults - use homogenization_isostrain - use homogenization_RGC +use prec, only: pReal,pInt +use math, only: math_I3 +use debug, only: debug_verbosity +use IO, only: IO_error, IO_open_file, IO_open_jobFile +use mesh, only: mesh_maxNips,mesh_NcpElems,mesh_element,FE_Nips +use material +use constitutive, only: constitutive_maxSizePostResults +use crystallite, only: crystallite_maxSizePostResults +use homogenization_isostrain +use homogenization_RGC - real(pReal) Temperature - integer(pInt), parameter :: fileunit = 200 - integer(pInt) e,i,g,p,myInstance,j - integer(pInt), dimension(:,:), pointer :: thisSize - character(len=64), dimension(:,:), pointer :: thisOutput - logical knownHomogenization +implicit none - if(.not. IO_open_file(fileunit,material_configFile)) call IO_error (100) ! corrupt config file - - call homogenization_isostrain_init(fileunit) ! parse all homogenizations of this type - call homogenization_RGC_init(fileunit) - - close(fileunit) - -! write description file for homogenization output - - if(.not. IO_open_jobFile(fileunit,'outputHomogenization')) call IO_error (50) ! problems in writing file - - do p = 1,material_Nhomogenization - i = homogenization_typeInstance(p) ! which instance of this homogenization type - knownHomogenization = .true. ! assume valid - select case(homogenization_type(p)) ! split per homogenization type - case (homogenization_isostrain_label) - thisOutput => homogenization_isostrain_output - thisSize => homogenization_isostrain_sizePostResult - case (homogenization_RGC_label) - thisOutput => homogenization_RGC_output - thisSize => homogenization_RGC_sizePostResult - case default - knownHomogenization = .false. - end select - - write(fileunit,*) - write(fileunit,'(a)') '['//trim(homogenization_name(p))//']' - write(fileunit,*) - if (knownHomogenization) then - write(fileunit,'(a)') '(type)'//char(9)//trim(homogenization_type(p)) - write(fileunit,'(a,i)') '(ngrains)'//char(9),homogenization_Ngrains(p) - do e = 1,homogenization_Noutput(p) - write(fileunit,'(a,i4)') trim(thisOutput(e,i))//char(9),thisSize(e,i) - enddo - endif - enddo - - close(fileunit) +real(pReal) Temperature +integer(pInt), parameter :: fileunit = 200 +integer(pInt) e,i,g,p,myInstance,j +integer(pInt), dimension(:,:), pointer :: thisSize +character(len=64), dimension(:,:), pointer :: thisOutput +logical knownHomogenization - allocate(homogenization_state0(mesh_maxNips,mesh_NcpElems)) - allocate(homogenization_subState0(mesh_maxNips,mesh_NcpElems)) - allocate(homogenization_state(mesh_maxNips,mesh_NcpElems)) - allocate(homogenization_sizeState(mesh_maxNips,mesh_NcpElems)); homogenization_sizeState = 0_pInt - allocate(homogenization_sizePostResults(mesh_maxNips,mesh_NcpElems)); homogenization_sizePostResults = 0_pInt +! --- PARSE HOMOGENIZATIONS FROM CONFIG FILE --- - allocate(materialpoint_dPdF(3,3,3,3,mesh_maxNips,mesh_NcpElems)); materialpoint_dPdF = 0.0_pReal - allocate(materialpoint_F0(3,3,mesh_maxNips,mesh_NcpElems)); - allocate(materialpoint_F(3,3,mesh_maxNips,mesh_NcpElems)); materialpoint_F = 0.0_pReal - allocate(materialpoint_subF0(3,3,mesh_maxNips,mesh_NcpElems)); materialpoint_subF0 = 0.0_pReal - allocate(materialpoint_subF(3,3,mesh_maxNips,mesh_NcpElems)); materialpoint_subF = 0.0_pReal - allocate(materialpoint_P(3,3,mesh_maxNips,mesh_NcpElems)); materialpoint_P = 0.0_pReal - allocate(materialpoint_Temperature(mesh_maxNips,mesh_NcpElems)); materialpoint_Temperature = Temperature - allocate(materialpoint_subFrac(mesh_maxNips,mesh_NcpElems)); materialpoint_subFrac = 0.0_pReal - allocate(materialpoint_subStep(mesh_maxNips,mesh_NcpElems)); materialpoint_subStep = 0.0_pReal - allocate(materialpoint_subdt(mesh_maxNips,mesh_NcpElems)); materialpoint_subdt = 0.0_pReal - allocate(materialpoint_requested(mesh_maxNips,mesh_NcpElems)); materialpoint_requested = .false. - allocate(materialpoint_converged(mesh_maxNips,mesh_NcpElems)); materialpoint_converged = .true. - allocate(materialpoint_doneAndHappy(2,mesh_maxNips,mesh_NcpElems)); materialpoint_doneAndHappy = .true. - - forall (i = 1:mesh_maxNips,e = 1:mesh_NcpElems) - materialpoint_F0(1:3,1:3,i,e) = math_I3 - materialpoint_F(1:3,1:3,i,e) = math_I3 - end forall - - do e = 1,mesh_NcpElems ! loop over elements - myInstance = homogenization_typeInstance(mesh_element(3,e)) - do i = 1,FE_Nips(mesh_element(2,e)) ! loop over IPs - select case(homogenization_type(mesh_element(3,e))) -!* isostrain - case (homogenization_isostrain_label) - if (homogenization_isostrain_sizeState(myInstance) > 0_pInt) then - allocate(homogenization_state0(i,e)%p(homogenization_isostrain_sizeState(myInstance))) - allocate(homogenization_subState0(i,e)%p(homogenization_isostrain_sizeState(myInstance))) - allocate(homogenization_state(i,e)%p(homogenization_isostrain_sizeState(myInstance))) - homogenization_state0(i,e)%p = homogenization_isostrain_stateInit(myInstance) - homogenization_sizeState(i,e) = homogenization_isostrain_sizeState(myInstance) - endif - homogenization_sizePostResults(i,e) = homogenization_isostrain_sizePostResults(myInstance) -!* RGC homogenization - case (homogenization_RGC_label) - if (homogenization_RGC_sizeState(myInstance) > 0_pInt) then - allocate(homogenization_state0(i,e)%p(homogenization_RGC_sizeState(myInstance))) - allocate(homogenization_subState0(i,e)%p(homogenization_RGC_sizeState(myInstance))) - allocate(homogenization_state(i,e)%p(homogenization_RGC_sizeState(myInstance))) - homogenization_state0(i,e)%p = homogenization_RGC_stateInit(myInstance) - homogenization_sizeState(i,e) = homogenization_RGC_sizeState(myInstance) - endif - homogenization_sizePostResults(i,e) = homogenization_RGC_sizePostResults(myInstance) - case default - call IO_error(201,ext_msg=homogenization_type(mesh_element(3,e))) ! unknown type 201 is homogenization! - end select - enddo - enddo - - homogenization_maxSizeState = maxval(homogenization_sizeState) - homogenization_maxSizePostResults = maxval(homogenization_sizePostResults) - - materialpoint_sizeResults = 1+ 1+homogenization_maxSizePostResults + & ! grain count, homogSize, homogResult - homogenization_maxNgrains*(1+crystallite_maxSizePostResults+ & ! results count, cryst results - 1+constitutive_maxSizePostResults) ! results count, constitutive results - allocate(materialpoint_results(materialpoint_sizeResults, mesh_maxNips,mesh_NcpElems)) +if(.not. IO_open_file(fileunit,material_configFile)) call IO_error (100) ! corrupt config file +call homogenization_isostrain_init(fileunit) +call homogenization_RGC_init(fileunit) +close(fileunit) + + +! --- WRITE DESCRIPTION FILE FOR HOMOGENIZATION OUTPUT --- + +if(.not. IO_open_jobFile(fileunit,'outputHomogenization')) then ! problems in writing file + call IO_error (50) +endif +do p = 1,material_Nhomogenization + i = homogenization_typeInstance(p) ! which instance of this homogenization type + knownHomogenization = .true. ! assume valid + select case(homogenization_type(p)) ! split per homogenization type + case (homogenization_isostrain_label) + thisOutput => homogenization_isostrain_output + thisSize => homogenization_isostrain_sizePostResult + case (homogenization_RGC_label) + thisOutput => homogenization_RGC_output + thisSize => homogenization_RGC_sizePostResult + case default + knownHomogenization = .false. + end select + write(fileunit,*) + write(fileunit,'(a)') '['//trim(homogenization_name(p))//']' + write(fileunit,*) + if (knownHomogenization) then + write(fileunit,'(a)') '(type)'//char(9)//trim(homogenization_type(p)) + write(fileunit,'(a,i)') '(ngrains)'//char(9),homogenization_Ngrains(p) + do e = 1,homogenization_Noutput(p) + write(fileunit,'(a,i4)') trim(thisOutput(e,i))//char(9),thisSize(e,i) + enddo + endif +enddo +close(fileunit) + + +! --- ALLOCATE AND INITIALIZE GLOBAL VARIABLES --- + +allocate(homogenization_state0(mesh_maxNips,mesh_NcpElems)) +allocate(homogenization_subState0(mesh_maxNips,mesh_NcpElems)) +allocate(homogenization_state(mesh_maxNips,mesh_NcpElems)) +allocate(homogenization_sizeState(mesh_maxNips,mesh_NcpElems)); homogenization_sizeState = 0_pInt +allocate(homogenization_sizePostResults(mesh_maxNips,mesh_NcpElems)); homogenization_sizePostResults = 0_pInt + +allocate(materialpoint_dPdF(3,3,3,3,mesh_maxNips,mesh_NcpElems)); materialpoint_dPdF = 0.0_pReal +allocate(materialpoint_F0(3,3,mesh_maxNips,mesh_NcpElems)); +allocate(materialpoint_F(3,3,mesh_maxNips,mesh_NcpElems)); materialpoint_F = 0.0_pReal +allocate(materialpoint_subF0(3,3,mesh_maxNips,mesh_NcpElems)); materialpoint_subF0 = 0.0_pReal +allocate(materialpoint_subF(3,3,mesh_maxNips,mesh_NcpElems)); materialpoint_subF = 0.0_pReal +allocate(materialpoint_P(3,3,mesh_maxNips,mesh_NcpElems)); materialpoint_P = 0.0_pReal +allocate(materialpoint_Temperature(mesh_maxNips,mesh_NcpElems)); materialpoint_Temperature = Temperature +allocate(materialpoint_subFrac(mesh_maxNips,mesh_NcpElems)); materialpoint_subFrac = 0.0_pReal +allocate(materialpoint_subStep(mesh_maxNips,mesh_NcpElems)); materialpoint_subStep = 0.0_pReal +allocate(materialpoint_subdt(mesh_maxNips,mesh_NcpElems)); materialpoint_subdt = 0.0_pReal +allocate(materialpoint_requested(mesh_maxNips,mesh_NcpElems)); materialpoint_requested = .false. +allocate(materialpoint_converged(mesh_maxNips,mesh_NcpElems)); materialpoint_converged = .true. +allocate(materialpoint_doneAndHappy(2,mesh_maxNips,mesh_NcpElems)); materialpoint_doneAndHappy = .true. + +forall (i = 1:mesh_maxNips,e = 1:mesh_NcpElems) + materialpoint_F0(1:3,1:3,i,e) = math_I3 + materialpoint_F(1:3,1:3,i,e) = math_I3 +end forall + + +! --- ALLOCATE AND INITIALIZE GLOBAL STATE AND POSTRESULTS VARIABLES --- + +!$OMP PARALLEL DO PRIVATE(myInstance) + do e = 1,mesh_NcpElems ! loop over elements + myInstance = homogenization_typeInstance(mesh_element(3,e)) + do i = 1,FE_Nips(mesh_element(2,e)) ! loop over IPs + select case(homogenization_type(mesh_element(3,e))) + case (homogenization_isostrain_label) + if (homogenization_isostrain_sizeState(myInstance) > 0_pInt) then + allocate(homogenization_state0(i,e)%p(homogenization_isostrain_sizeState(myInstance))) + allocate(homogenization_subState0(i,e)%p(homogenization_isostrain_sizeState(myInstance))) + allocate(homogenization_state(i,e)%p(homogenization_isostrain_sizeState(myInstance))) + homogenization_state0(i,e)%p = homogenization_isostrain_stateInit(myInstance) + homogenization_sizeState(i,e) = homogenization_isostrain_sizeState(myInstance) + endif + homogenization_sizePostResults(i,e) = homogenization_isostrain_sizePostResults(myInstance) + case (homogenization_RGC_label) + if (homogenization_RGC_sizeState(myInstance) > 0_pInt) then + allocate(homogenization_state0(i,e)%p(homogenization_RGC_sizeState(myInstance))) + allocate(homogenization_subState0(i,e)%p(homogenization_RGC_sizeState(myInstance))) + allocate(homogenization_state(i,e)%p(homogenization_RGC_sizeState(myInstance))) + homogenization_state0(i,e)%p = homogenization_RGC_stateInit(myInstance) + homogenization_sizeState(i,e) = homogenization_RGC_sizeState(myInstance) + endif + homogenization_sizePostResults(i,e) = homogenization_RGC_sizePostResults(myInstance) + case default + call IO_error(201,ext_msg=homogenization_type(mesh_element(3,e))) ! unknown type 201 is homogenization! + end select + enddo + enddo +!$OMP END PARALLEL DO +homogenization_maxSizeState = maxval(homogenization_sizeState) +homogenization_maxSizePostResults = maxval(homogenization_sizePostResults) +materialpoint_sizeResults = 1 & ! grain count + + 1 + homogenization_maxSizePostResults & ! homogSize & homogResult + + homogenization_maxNgrains * (1 + crystallite_maxSizePostResults & ! crystallite size & crystallite results + + 1 + constitutive_maxSizePostResults) ! constitutive size & constitutive results +allocate(materialpoint_results(materialpoint_sizeResults,mesh_maxNips,mesh_NcpElems)) -! *** Output to MARC output file *** !$OMP CRITICAL (write2out) - write(6,*) - write(6,*) '<<<+- homogenization init -+>>>' - write(6,*) '$Id$' - write(6,*) - if (debug_verbosity > 0) then - write(6,'(a32,x,7(i5,x))') 'homogenization_state0: ', shape(homogenization_state0) - write(6,'(a32,x,7(i5,x))') 'homogenization_subState0: ', shape(homogenization_subState0) - write(6,'(a32,x,7(i5,x))') 'homogenization_state: ', shape(homogenization_state) - write(6,'(a32,x,7(i5,x))') 'homogenization_sizeState: ', shape(homogenization_sizeState) - write(6,'(a32,x,7(i5,x))') 'homogenization_sizePostResults: ', shape(homogenization_sizePostResults) - write(6,*) - write(6,'(a32,x,7(i5,x))') 'materialpoint_dPdF: ', shape(materialpoint_dPdF) - write(6,'(a32,x,7(i5,x))') 'materialpoint_F0: ', shape(materialpoint_F0) - write(6,'(a32,x,7(i5,x))') 'materialpoint_F: ', shape(materialpoint_F) - write(6,'(a32,x,7(i5,x))') 'materialpoint_subF0: ', shape(materialpoint_subF0) - write(6,'(a32,x,7(i5,x))') 'materialpoint_subF: ', shape(materialpoint_subF) - write(6,'(a32,x,7(i5,x))') 'materialpoint_P: ', shape(materialpoint_P) - write(6,'(a32,x,7(i5,x))') 'materialpoint_Temperature: ', shape(materialpoint_Temperature) - write(6,'(a32,x,7(i5,x))') 'materialpoint_subFrac: ', shape(materialpoint_subFrac) - write(6,'(a32,x,7(i5,x))') 'materialpoint_subStep: ', shape(materialpoint_subStep) - write(6,'(a32,x,7(i5,x))') 'materialpoint_subdt: ', shape(materialpoint_subdt) - write(6,'(a32,x,7(i5,x))') 'materialpoint_requested: ', shape(materialpoint_requested) - write(6,'(a32,x,7(i5,x))') 'materialpoint_converged: ', shape(materialpoint_converged) - write(6,'(a32,x,7(i5,x))') 'materialpoint_doneAndHappy: ', shape(materialpoint_doneAndHappy) - write(6,*) - write(6,'(a32,x,7(i5,x))') 'materialpoint_results: ', shape(materialpoint_results) - write(6,*) - write(6,'(a32,x,7(i5,x))') 'maxSizeState: ', homogenization_maxSizeState - write(6,'(a32,x,7(i5,x))') 'maxSizePostResults: ', homogenization_maxSizePostResults - endif - call flush(6) + write(6,*) + write(6,*) '<<<+- homogenization init -+>>>' + write(6,*) '$Id$' + write(6,*) + if (debug_verbosity > 0) then + write(6,'(a32,x,7(i5,x))') 'homogenization_state0: ', shape(homogenization_state0) + write(6,'(a32,x,7(i5,x))') 'homogenization_subState0: ', shape(homogenization_subState0) + write(6,'(a32,x,7(i5,x))') 'homogenization_state: ', shape(homogenization_state) + write(6,'(a32,x,7(i5,x))') 'homogenization_sizeState: ', shape(homogenization_sizeState) + write(6,'(a32,x,7(i5,x))') 'homogenization_sizePostResults: ', shape(homogenization_sizePostResults) + write(6,*) + write(6,'(a32,x,7(i5,x))') 'materialpoint_dPdF: ', shape(materialpoint_dPdF) + write(6,'(a32,x,7(i5,x))') 'materialpoint_F0: ', shape(materialpoint_F0) + write(6,'(a32,x,7(i5,x))') 'materialpoint_F: ', shape(materialpoint_F) + write(6,'(a32,x,7(i5,x))') 'materialpoint_subF0: ', shape(materialpoint_subF0) + write(6,'(a32,x,7(i5,x))') 'materialpoint_subF: ', shape(materialpoint_subF) + write(6,'(a32,x,7(i5,x))') 'materialpoint_P: ', shape(materialpoint_P) + write(6,'(a32,x,7(i5,x))') 'materialpoint_Temperature: ', shape(materialpoint_Temperature) + write(6,'(a32,x,7(i5,x))') 'materialpoint_subFrac: ', shape(materialpoint_subFrac) + write(6,'(a32,x,7(i5,x))') 'materialpoint_subStep: ', shape(materialpoint_subStep) + write(6,'(a32,x,7(i5,x))') 'materialpoint_subdt: ', shape(materialpoint_subdt) + write(6,'(a32,x,7(i5,x))') 'materialpoint_requested: ', shape(materialpoint_requested) + write(6,'(a32,x,7(i5,x))') 'materialpoint_converged: ', shape(materialpoint_converged) + write(6,'(a32,x,7(i5,x))') 'materialpoint_doneAndHappy: ', shape(materialpoint_doneAndHappy) + write(6,*) + write(6,'(a32,x,7(i5,x))') 'materialpoint_results: ', shape(materialpoint_results) + write(6,*) + write(6,'(a32,x,7(i5,x))') 'maxSizeState: ', homogenization_maxSizeState + write(6,'(a32,x,7(i5,x))') 'maxSizePostResults: ', homogenization_maxSizePostResults + endif + call flush(6) !$OMP END CRITICAL (write2out) - return - endsubroutine @@ -335,13 +339,13 @@ subroutine materialpoint_stressAndItsTangent(& do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) ! iterate over IPs of this element to be processed if ( materialpoint_converged(i,e) ) then +#ifndef _OPENMP if (debug_verbosity > 2 .and. ((e == debug_e .and. i == debug_i) .or. .not. debug_selectiveDebugger)) then - !$OMP CRITICAL (write2out) - write(6,'(a,x,f10.8,x,a,x,f10.8,x,a,/)') '<< HOMOG >> winding forward from', & - materialpoint_subFrac(i,e), 'to current materialpoint_subFrac', & - materialpoint_subFrac(i,e)+materialpoint_subStep(i,e),'in materialpoint_stressAndItsTangent' - !$OMP END CRITICAL (write2out) + write(6,'(a,x,f10.8,x,a,x,f10.8,x,a,/)') '<< HOMOG >> winding forward from', & + materialpoint_subFrac(i,e), 'to current materialpoint_subFrac', & + materialpoint_subFrac(i,e)+materialpoint_subStep(i,e),'in materialpoint_stressAndItsTangent' endif +#endif ! calculate new subStep and new subFrac materialpoint_subFrac(i,e) = materialpoint_subFrac(i,e) + materialpoint_subStep(i,e) @@ -386,12 +390,12 @@ subroutine materialpoint_stressAndItsTangent(& materialpoint_subStep(i,e) = subStepSizeHomog * materialpoint_subStep(i,e) ! crystallite had severe trouble, so do a significant cutback !$OMP FLUSH(materialpoint_subStep) +#ifndef _OPENMP if (debug_verbosity > 2 .and. ((e == debug_e .and. i == debug_i) .or. .not. debug_selectiveDebugger)) then - !$OMP CRITICAL (write2out) - write(6,'(a,x,f10.8,/)') '<< HOMOG >> cutback step in materialpoint_stressAndItsTangent with new materialpoint_subStep:',& - materialpoint_subStep(i,e) - !$OMP END CRITICAL (write2out) + write(6,'(a,x,f10.8,/)') '<< HOMOG >> cutback step in materialpoint_stressAndItsTangent with new materialpoint_subStep:',& + materialpoint_subStep(i,e) endif +#endif ! restore... crystallite_Temperature(1:myNgrains,i,e) = crystallite_partionedTemperature0(1:myNgrains,i,e) ! ...temperatures @@ -495,8 +499,7 @@ subroutine materialpoint_stressAndItsTangent(& enddo ! cutback loop - if (.not. terminallyIll ) then - + if (.not. terminallyIll ) then call crystallite_orientations() ! calculate crystal orientations !$OMP PARALLEL DO do e = FEsolving_execElem(1),FEsolving_execElem(2) ! iterate over elements to be processed @@ -505,21 +508,12 @@ subroutine materialpoint_stressAndItsTangent(& call homogenization_averageTemperature(i,e) enddo; enddo !$OMP END PARALLEL DO - - if (debug_verbosity > 2) then - !$OMP CRITICAL (write2out) - write (6,*) - write (6,'(a)') '<< HOMOG >> Material Point end' - write (6,*) - !$OMP END CRITICAL (write2out) - endif else !$OMP CRITICAL (write2out) - write (6,*) - write (6,'(a)') '<< HOMOG >> Material Point terminally ill' - write (6,*) + write (6,*) + write (6,'(a)') '<< HOMOG >> Material Point terminally ill' + write (6,*) !$OMP END CRITICAL (write2out) - endif return @@ -540,25 +534,31 @@ subroutine materialpoint_postResults(dt) implicit none real(pReal), intent(in) :: dt - integer(pInt) g,i,e,c,d,myNgrains,myCrystallite + integer(pInt) g,i,e,pos,size,myNgrains,myCrystallite - !$OMP PARALLEL DO PRIVATE(myNgrains,myCrystallite,c,d) + !$OMP PARALLEL DO PRIVATE(myNgrains,myCrystallite,pos,size) do e = FEsolving_execElem(1),FEsolving_execElem(2) ! iterate over elements to be processed myNgrains = homogenization_Ngrains(mesh_element(3,e)) myCrystallite = microstructure_crystallite(mesh_element(4,e)) do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) ! iterate over IPs of this element to be processed - c = 0_pInt - materialpoint_results(c+1,i,e) = myNgrains; c = c+1_pInt ! tell number of grains at materialpoint - d = homogenization_sizePostResults(i,e) - materialpoint_results(c+1,i,e) = d; c = c+1_pInt ! tell size of homogenization results - if (d > 0_pInt) then ! any homogenization results to mention? - materialpoint_results(c+1:c+d,i,e) = & ! tell homogenization results - homogenization_postResults(i,e); c = c+d + pos = 0_pInt + + materialpoint_results(pos+1,i,e) = myNgrains ! tell number of grains at materialpoint + pos = pos + 1_pInt + + size = homogenization_sizePostResults(i,e) + materialpoint_results(pos+1,i,e) = size ! tell size of homogenization results + pos = pos + 1_pInt + + if (size > 0_pInt) then ! any homogenization results to mention? + materialpoint_results(pos+1:pos+size,i,e) = homogenization_postResults(i,e) ! tell homogenization results + pos = pos + size endif + do g = 1,myNgrains ! loop over all grains - d = 1+crystallite_sizePostResults(myCrystallite) + 1+constitutive_sizePostResults(g,i,e) - materialpoint_results(c+1:c+d,i,e) = & ! tell crystallite results - crystallite_postResults(dt,g,i,e); c = c+d + size = (1 + crystallite_sizePostResults(myCrystallite)) + (1 + constitutive_sizePostResults(g,i,e)) + materialpoint_results(pos+1:pos+size,i,e) = crystallite_postResults(dt,g,i,e) ! tell crystallite results + pos = pos + size enddo enddo enddo diff --git a/code/mpie_cpfem_abaqus_exp.f b/code/mpie_cpfem_abaqus_exp.f index 7b0a3c249..b42025838 100644 --- a/code/mpie_cpfem_abaqus_exp.f +++ b/code/mpie_cpfem_abaqus_exp.f @@ -14,7 +14,7 @@ ! !******************************************************************** -include "prec.f90" ! uses nothing else +#include "prec.f90" MODULE mpie_interface @@ -73,25 +73,25 @@ end function END MODULE -include "IO.f90" ! uses prec -include "numerics.f90" ! uses prec, IO -include "debug.f90" ! uses prec, numerics -include "math.f90" ! uses prec, numerics, debug -include "FEsolving.f90" ! uses prec, IO, debug -include "mesh.f90" ! uses prec, math, IO, FEsolving, debug -include "material.f90" ! uses prec, math, IO, mesh, debug -include "lattice.f90" ! uses prec, math, IO, material, debug -include "constitutive_j2.f90" ! uses prec, math, IO, lattice, material, debug -include "constitutive_phenopowerlaw.f90" ! uses prec, math, IO, lattice, material, debug -include "constitutive_titanmod.f90" ! uses prec, math, IO, lattice, material, debug -include "constitutive_dislotwin.f90" ! uses prec, math, IO, lattice, material, debug -include "constitutive_nonlocal.f90" ! uses prec, math, IO, lattice, material, debug -include "constitutive.f90" ! uses prec, IO, math, lattice, mesh, debug -include "crystallite.f90" ! uses prec, math, IO, numerics, Fesolving, material, mesh, constitutive, debug -include "homogenization_isostrain.f90" ! uses prec, math, IO, debug -include "homogenization_RGC.f90" ! uses prec, math, IO, numerics, mesh, debug -include "homogenization.f90" ! uses prec, math, IO, numerics, debug -include "CPFEM.f90" ! uses prec, math, IO, numerics, debug, FEsolving, mesh, lattice, constitutive, crystallite, debug +#include "IO.f90" +#include "numerics.f90" +#include "debug.f90" +#include "math.f90" +#include "FEsolving.f90" +#include "mesh.f90" +#include "material.f90" +#include "lattice.f90" +#include "constitutive_j2.f90" +#include "constitutive_phenopowerlaw.f90" +#include "constitutive_titanmod.f90" +#include "constitutive_dislotwin.f90" +#include "constitutive_nonlocal.f90" +#include "constitutive.f90" +#include "crystallite.f90" +#include "homogenization_isostrain.f90" +#include "homogenization_RGC.f90" +#include "homogenization.f90" +#include "CPFEM.f90" subroutine vumat (jblock, ndir, nshr, nstatev, nfieldv, nprops, lanneal, & stepTime, totalTime, dt, cmname, coordMp, charLength, & diff --git a/code/mpie_cpfem_abaqus_std.f b/code/mpie_cpfem_abaqus_std.f index 3e2eba764..12a689a41 100644 --- a/code/mpie_cpfem_abaqus_std.f +++ b/code/mpie_cpfem_abaqus_std.f @@ -14,7 +14,7 @@ ! !******************************************************************** -include "prec.f90" ! uses nothing else +#include "prec.f90" MODULE mpie_interface @@ -75,25 +75,25 @@ end function END MODULE -include "IO.f90" ! uses prec -include "numerics.f90" ! uses prec, IO -include "debug.f90" ! uses prec, numerics -include "math.f90" ! uses prec, numerics, debug -include "FEsolving.f90" ! uses prec, IO, debug -include "mesh.f90" ! uses prec, math, IO, FEsolving, debug -include "material.f90" ! uses prec, math, IO, mesh, debug -include "lattice.f90" ! uses prec, math, IO, material, debug -include "constitutive_j2.f90" ! uses prec, math, IO, lattice, material, debug -include "constitutive_phenopowerlaw.f90" ! uses prec, math, IO, lattice, material, debug -include "constitutive_titanmod.f90" ! uses prec, math, IO, lattice, material, debug -include "constitutive_dislotwin.f90" ! uses prec, math, IO, lattice, material, debug -include "constitutive_nonlocal.f90" ! uses prec, math, IO, lattice, material, debug -include "constitutive.f90" ! uses prec, IO, math, lattice, mesh, debug -include "crystallite.f90" ! uses prec, math, IO, numerics, Fesolving, material, mesh, constitutive, debug -include "homogenization_isostrain.f90" ! uses prec, math, IO, debug -include "homogenization_RGC.f90" ! uses prec, math, IO, numerics, mesh, debug -include "homogenization.f90" ! uses prec, math, IO, numerics, debug -include "CPFEM.f90" ! uses prec, math, IO, numerics, debug, FEsolving, mesh, lattice, constitutive, crystallite, debug +#include "IO.f90" +#include "numerics.f90" +#include "debug.f90" +#include "math.f90" +#include "FEsolving.f90" +#include "mesh.f90" +#include "material.f90" +#include "lattice.f90" +#include "constitutive_j2.f90" +#include "constitutive_phenopowerlaw.f90" +#include "constitutive_titanmod.f90" +#include "constitutive_dislotwin.f90" +#include "constitutive_nonlocal.f90" +#include "constitutive.f90" +#include "crystallite.f90" +#include "homogenization_isostrain.f90" +#include "homogenization_RGC.f90" +#include "homogenization.f90" +#include "CPFEM.f90" subroutine UMAT(STRESS,STATEV,DDSDDE,SSE,SPD,SCD,& RPL,DDSDDT,DRPLDE,DRPLDT,STRAN,DSTRAN,& diff --git a/code/mpie_cpfem_marc.f90 b/code/mpie_cpfem_marc.f90 index 6bb0f9d10..7ac121b5e 100644 --- a/code/mpie_cpfem_marc.f90 +++ b/code/mpie_cpfem_marc.f90 @@ -36,7 +36,7 @@ ! - creeps: timinc !******************************************************************** ! -include "prec.f90" ! uses nothing else +#include "prec.f90" MODULE mpie_interface @@ -100,25 +100,25 @@ end function END MODULE -include "IO.f90" ! uses prec -include "numerics.f90" ! uses prec, IO -include "debug.f90" ! uses prec, numerics -include "math.f90" ! uses prec, numerics, debug -include "FEsolving.f90" ! uses prec, IO, debug -include "mesh.f90" ! uses prec, math, IO, FEsolving, debug -include "material.f90" ! uses prec, math, IO, mesh, debug -include "lattice.f90" ! uses prec, math, IO, material, debug -include "constitutive_j2.f90" ! uses prec, math, IO, lattice, material, debug -include "constitutive_phenopowerlaw.f90" ! uses prec, math, IO, lattice, material, debug -include "constitutive_titanmod.f90" ! uses prec, math, IO, lattice, material, debug -include "constitutive_dislotwin.f90" ! uses prec, math, IO, lattice, material, debug -include "constitutive_nonlocal.f90" ! uses prec, math, IO, lattice, material, debug -include "constitutive.f90" ! uses prec, IO, math, lattice, mesh, debug -include "crystallite.f90" ! uses prec, math, IO, numerics, Fesolving, material, mesh, constitutive, debug -include "homogenization_isostrain.f90" ! uses prec, math, IO, debug -include "homogenization_RGC.f90" ! uses prec, math, IO, numerics, mesh, debug -include "homogenization.f90" ! uses prec, math, IO, numerics, debug -include "CPFEM.f90" ! uses prec, math, IO, numerics, debug, FEsolving, mesh, lattice, constitutive, crystallite, debug +#include "IO.f90" +#include "numerics.f90" +#include "debug.f90" +#include "math.f90" +#include "FEsolving.f90" +#include "mesh.f90" +#include "material.f90" +#include "lattice.f90" +#include "constitutive_j2.f90" +#include "constitutive_phenopowerlaw.f90" +#include "constitutive_titanmod.f90" +#include "constitutive_dislotwin.f90" +#include "constitutive_nonlocal.f90" +#include "constitutive.f90" +#include "crystallite.f90" +#include "homogenization_isostrain.f90" +#include "homogenization_RGC.f90" +#include "homogenization.f90" +#include "CPFEM.f90" !********************************************************************