From 40b1478dacf8dd29ebb5717537e090a0cae32cc2 Mon Sep 17 00:00:00 2001 From: Denny Tjahjanto Date: Wed, 24 Mar 2010 13:20:12 +0000 Subject: [PATCH] the latest RGC model + corrections for "element homogeneous" feature --- code/CPFEM.f90 | 21 ++- code/FEsolving.f90 | 2 +- code/constitutive.f90 | 45 ++++- code/constitutive_nonlocal.f90 | 17 +- code/crystallite.f90 | 86 +++++----- code/debug.f90 | 2 - code/homogenization.f90 | 3 + code/homogenization_RGC.f90 | 304 ++++++++++++++++++++++++--------- code/material.f90 | 2 + code/math.f90 | 26 ++- code/mesh.f90 | 2 +- code/mpie_cpfem_marc.f90 | 1 - code/numerics.f90 | 15 +- 13 files changed, 364 insertions(+), 162 deletions(-) diff --git a/code/CPFEM.f90 b/code/CPFEM.f90 index 55941ff93..2f6959d05 100644 --- a/code/CPFEM.f90 +++ b/code/CPFEM.f90 @@ -86,7 +86,6 @@ subroutine CPFEM_general(mode, ffn, ffn1, Temperature, dt, element, IP, cauchySt cycleCounter, & theInc, & theTime, & - theDelta, & FEsolving_execElem, & FEsolving_execIP use math, only: math_init, & @@ -102,11 +101,14 @@ subroutine CPFEM_general(mode, ffn, ffn1, Temperature, dt, element, IP, cauchySt mesh_maxNips use lattice, only: lattice_init use material, only: material_init, & - homogenization_maxNgrains + homogenization_maxNgrains, & + microstructure_elemhomo use constitutive, only: constitutive_init,& constitutive_state0,constitutive_state use crystallite, only: crystallite_init, & crystallite_F0, & + crystallite_rexParm, & + crystallite_critVal, & crystallite_partionedF, & crystallite_Fp0, & crystallite_Fp, & @@ -122,6 +124,7 @@ subroutine CPFEM_general(mode, ffn, ffn1, Temperature, dt, element, IP, cauchySt materialpoint_F0, & materialpoint_P, & materialpoint_dPdF, & + materialpoint_results, & materialpoint_Temperature, & materialpoint_stressAndItsTangent, & materialpoint_postResults @@ -286,6 +289,16 @@ subroutine CPFEM_general(mode, ffn, ffn1, Temperature, dt, element, IP, cauchySt elseif (.not. CPFEM_calc_done) then 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(:,:,i,e) = materialpoint_P(:,:,1,e) + materialpoint_F(:,:,i,e) = materialpoint_F(:,:,1,e) + materialpoint_dPdF(:,:,:,:,i,e) = materialpoint_dPdF(:,:,:,:,1,e) + materialpoint_results(:,i,e) = materialpoint_results(:,1,e) + end forall + endif + enddo CPFEM_calc_done = .true. endif @@ -321,8 +334,8 @@ subroutine CPFEM_general(mode, ffn, ffn1, Temperature, dt, element, IP, cauchySt else if (mode == 5) then CPFEM_dcsde = CPFEM_dcsde_knownGood ! --+>> RESTORE CONSISTENT JACOBIAN FROM FORMER CONVERGED INC end if - call random_number(rnd) - if (rnd < 0.5_pReal) rnd = 1.0_pReal - rnd + call random_number(rnd) + if (rnd < 0.5_pReal) rnd = 1.0_pReal - rnd materialpoint_Temperature(IP,cp_en) = Temperature materialpoint_F0(:,:,IP,cp_en) = ffn materialpoint_F(:,:,IP,cp_en) = ffn1 diff --git a/code/FEsolving.f90 b/code/FEsolving.f90 index c08efadc4..6dd30548d 100644 --- a/code/FEsolving.f90 +++ b/code/FEsolving.f90 @@ -1,4 +1,4 @@ -!* $Id$ +!* 'Id' !############################################################## MODULE FEsolving !############################################################## diff --git a/code/constitutive.f90 b/code/constitutive.f90 index aeb277bad..f6f10f924 100644 --- a/code/constitutive.f90 +++ b/code/constitutive.f90 @@ -32,6 +32,7 @@ CONTAINS !**************************************** !* - constitutive_init !* - constitutive_homogenizedC +!* - constitutive_averageBurgers !* - constitutive_microstructure !* - constitutive_LpAndItsTangent !* - constitutive_collectDotState @@ -268,8 +269,48 @@ function constitutive_homogenizedC(ipc,ip,el) return endfunction +function constitutive_averageBurgers(ipc,ip,el) +!********************************************************************* +!* This function returns the average length of Burgers vector * +!* INPUT: * +!* - state : state variables * +!* - ipc : component-ID of current integration point * +!* - ip : current integration point * +!* - el : current element * +!********************************************************************* + use prec, only: pReal,pInt + use material, only: phase_constitution,material_phase + use constitutive_j2 + use constitutive_phenopowerlaw + use constitutive_dislotwin + use constitutive_nonlocal + implicit none -subroutine constitutive_microstructure(Temperature,Tstar_v,Fe,Fp,ipc,ip,el) + !* Definition of variables + integer(pInt) ipc,ip,el + real(pReal) :: constitutive_averageBurgers + + select case (phase_constitution(material_phase(ipc,ip,el))) + + case (constitutive_j2_label) + constitutive_averageBurgers = 2.5e-10_pReal !constitutive_j2_averageBurgers(constitutive_state,ipc,ip,el) + + case (constitutive_phenopowerlaw_label) + constitutive_averageBurgers = 2.5e-10_pReal !constitutive_phenopowerlaw_averageBurgers(constitutive_state,ipc,ip,el) + + case (constitutive_dislotwin_label) + constitutive_averageBurgers = 2.5e-10_pReal !constitutive_dislotwin_averageBurgers(constitutive_state,ipc,ip,el) + + case (constitutive_nonlocal_label) + constitutive_averageBurgers = 2.5e-10_pReal !constitutive_nonlocal_averageBurgers(constitutive_state,ipc,ip,el) + + end select + +return +endfunction + + +subroutine constitutive_microstructure(Temperature,Fe,Fp,ipc,ip,el) !********************************************************************* !* This function calculates from state needed variables * !* INPUT: * @@ -538,4 +579,4 @@ return endfunction -END MODULE \ No newline at end of file +END MODULE diff --git a/code/constitutive_nonlocal.f90 b/code/constitutive_nonlocal.f90 index def0e5817..2592e12f2 100644 --- a/code/constitutive_nonlocal.f90 +++ b/code/constitutive_nonlocal.f90 @@ -73,10 +73,10 @@ real(pReal), dimension(:,:), allocatable :: constitutive_nonlocal_ constitutive_nonlocal_lambda0PerSlipSystem, & ! mean free path prefactor for each slip system and instance constitutive_nonlocal_burgersPerSlipFamily, & ! absolute length of burgers vector [m] for each family and instance constitutive_nonlocal_burgersPerSlipSystem, & ! absolute length of burgers vector [m] for each slip system and instance - constitutive_nonlocal_dLowerEdgePerSlipFamily, & ! minimum stable edge dipole height for each family and instance - constitutive_nonlocal_dLowerEdgePerSlipSystem, & ! minimum stable edge dipole height for each slip system and instance - constitutive_nonlocal_dLowerScrewPerSlipFamily, & ! minimum stable screw dipole height for each family and instance - constitutive_nonlocal_dLowerScrewPerSlipSystem, & ! minimum stable screw dipole height for each slip system and instance + constitutive_nonlocal_dLowerEdgePerSlipFamily, & ! minimum stable edge dipole height for each family and instance + constitutive_nonlocal_dLowerEdgePerSlipSystem, & ! minimum stable edge dipole height for each slip system and instance + constitutive_nonlocal_dLowerScrewPerSlipFamily, & ! minimum stable screw dipole height for each family and instance + constitutive_nonlocal_dLowerScrewPerSlipSystem, & ! minimum stable screw dipole height for each slip system and instance constitutive_nonlocal_interactionSlipSlip ! coefficients for slip-slip interaction for each interaction type and instance real(pReal), dimension(:,:,:,:,:), allocatable :: constitutive_nonlocal_v, & ! dislocation velocity constitutive_nonlocal_rhoDotFlux ! dislocation convection term @@ -1230,17 +1230,9 @@ integer(pInt) myInstance, & ! current neighboring_ip, & ! integration point of my neighbor c, & ! character of dislocation n, & ! index of my current neighbor - opposite_n, & ! index of my opposite neighbor - opposite_ip, & ! ip of my opposite neighbor - opposite_el, & ! element index of my opposite neighbor t, & ! type of dislocation s, & ! index of my current slip system sLattice ! index of my current slip system according to lattice order -real(pReal), dimension(constitutive_nonlocal_totalNslip(phase_constitutionInstance(material_phase(g,ip,el))),8) :: & - rhoSgl, & ! current single dislocation densities (positive/negative screw and edge without dipoles) - previousRhoSgl, & ! previous single dislocation densities (positive/negative screw and edge without dipoles) - totalRhoDotSgl, & ! total rate of change of single dislocation densities - thisRhoDotSgl ! rate of change of single dislocation densities for this mechanism real(pReal), dimension(constitutive_nonlocal_totalNslip(phase_constitutionInstance(material_phase(g,ip,el))),4) :: & fluxdensity, & ! flux density at central material point neighboring_fluxdensity, &! flux density at neighbroing material point @@ -1428,7 +1420,6 @@ fluxdensity = rhoSgl(:,1:4) * constitutive_nonlocal_v(:,:,g,ip,el) do n = 1,FE_NipNeighbors(mesh_element(2,el)) ! loop through my neighbors opposite_n = n - 1_pInt + 2_pInt*mod(n,2_pInt) - neighboring_el = mesh_ipNeighborhood(1,n,ip,el) neighboring_ip = mesh_ipNeighborhood(2,n,ip,el) diff --git a/code/crystallite.f90 b/code/crystallite.f90 index 3e4ff6d9e..d3c721cf3 100644 --- a/code/crystallite.f90 +++ b/code/crystallite.f90 @@ -291,7 +291,8 @@ subroutine crystallite_init(Temperature) write(6,'(a35,x,7(i5,x))') 'crystallite_subTstar0_v: ', shape(crystallite_subTstar0_v) write(6,'(a35,x,7(i5,x))') 'crystallite_dPdF: ', shape(crystallite_dPdF) write(6,'(a35,x,7(i5,x))') 'crystallite_fallbackdPdF: ', shape(crystallite_fallbackdPdF) - write(6,'(a35,x,7(i5,x))') 'crystallite_orientation: ', shape(crystallite_orientation) + write(6,'(a35,x,7(i5,x))') 'crystallite_R: ', shape(crystallite_R) + write(6,'(a35,x,7(i5,x))') 'crystallite_eulerangles: ', shape(crystallite_eulerangles) write(6,'(a35,x,7(i5,x))') 'crystallite_misorientation: ', shape(crystallite_misorientation) write(6,'(a35,x,7(i5,x))') 'crystallite_dt: ', shape(crystallite_dt) write(6,'(a35,x,7(i5,x))') 'crystallite_subdt: ', shape(crystallite_subdt) @@ -299,12 +300,11 @@ subroutine crystallite_init(Temperature) write(6,'(a35,x,7(i5,x))') 'crystallite_subStep: ', shape(crystallite_subStep) write(6,'(a35,x,7(i5,x))') 'crystallite_localConstitution: ', shape(crystallite_localConstitution) write(6,'(a35,x,7(i5,x))') 'crystallite_requested: ', shape(crystallite_requested) - write(6,'(a35,x,7(i5,x))') 'crystallite_todo: ', shape(crystallite_todo) + write(6,'(a35,x,7(i5,x))') 'crystallite_onTrack: ', shape(crystallite_onTrack) write(6,'(a35,x,7(i5,x))') 'crystallite_converged: ', shape(crystallite_converged) write(6,'(a35,x,7(i5,x))') 'crystallite_stateConverged: ', shape(crystallite_stateConverged) write(6,'(a35,x,7(i5,x))') 'crystallite_temperatureConverged: ', shape(crystallite_temperatureConverged) - write(6,'(a35,x,7(i5,x))') 'crystallite_sizePostResults: ', shape(crystallite_sizePostResults) - write(6,'(a35,x,7(i5,x))') 'crystallite_sizePostResult: ', shape(crystallite_sizePostResult) + write(6,'(a35,x,7(i5,x))') 'crystallite_todo: ', shape(crystallite_todo) write(6,*) write(6,*) 'Number of nonlocal grains: ',count(.not. crystallite_localConstitution) call flush(6) @@ -333,14 +333,8 @@ subroutine crystallite_stressAndItsTangent(updateJaco) pert_Fg, & pert_method, & nState, & - nCryst, & - iJacoStiffness + nCryst use debug, only: debugger, & - selectiveDebugger, & - verboseDebugger, & - debug_e, & - debug_i, & - debug_g, & debug_CrystalliteLoopDistribution, & debug_CrystalliteStateLoopDistribution, & debug_StiffnessStateLoopDistribution @@ -351,7 +345,8 @@ subroutine crystallite_stressAndItsTangent(updateJaco) math_Mandel6to33, & math_Mandel33to6, & math_I3, & - math_Plain3333to99 + math_Plain3333to99, & + math_EulerToR use FEsolving, only: FEsolving_execElem, & FEsolving_execIP, & theInc, & @@ -385,7 +380,8 @@ subroutine crystallite_stressAndItsTangent(updateJaco) !*** local variables ***! real(pReal) myTemperature, & ! local copy of the temperature - myPert ! perturbation with correct sign + myPert, & ! perturbation with correct sign + rexCritStrain ! perturbation with correct sign real(pReal), dimension(3,3) :: invFp, & ! inverse of the plastic deformation gradient Fe_guess, & ! guess for elastic deformation gradient Tstar ! 2nd Piola-Kirchhoff stress tensor @@ -401,8 +397,7 @@ subroutine crystallite_stressAndItsTangent(updateJaco) myNgrains, & mySizeState, & mySizeDotState - integer(pInt), dimension(2,9,homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems) :: & - kl + integer(pInt), dimension(2,9) :: kl logical onTrack, & ! flag indicating whether we are still on track temperatureConverged, & ! flag indicating if temperature converged stateConverged, & ! flag indicating if state converged @@ -412,8 +407,7 @@ subroutine crystallite_stressAndItsTangent(updateJaco) real(pReal), dimension(constitutive_maxSizeDotState) :: delta_dotState1, & ! difference between current and previous dotstate delta_dotState2 ! difference between previousDotState and previousDotState2 real(pReal) dot_prod12, & - dot_prod22, & - formerSubStep + dot_prod22 real(pReal), dimension(3,3,homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems) :: & storedF, & storedFp, & @@ -432,9 +426,6 @@ subroutine crystallite_stressAndItsTangent(updateJaco) logical, dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems) :: & storedConvergenceFlag logical, dimension(3,3) :: mask - logical forceLocalStiffnessCalculation ! flag indicating that stiffness calculation is always done locally - forceLocalStiffnessCalculation = .false. - ! ------ initialize to starting condition ------ @@ -451,17 +442,17 @@ subroutine crystallite_stressAndItsTangent(updateJaco) crystallite_subStep = 0.0_pReal !$OMP PARALLEL DO - do e = FEsolving_execElem(1),FEsolving_execElem(2) ! iterate over elements to be processed + do e = FEsolving_execElem(1),FEsolving_execElem(2) ! iterate over elements to be processed myNgrains = homogenization_Ngrains(mesh_element(3,e)) - do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) ! iterate over IPs of this element to be processed + do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) ! iterate over IPs of this element to be processed do g = 1,myNgrains - if (crystallite_requested(g,i,e)) then ! initialize restoration point of ... - crystallite_subTemperature0(g,i,e) = crystallite_partionedTemperature0(g,i,e) ! ...temperature - constitutive_subState0(g,i,e)%p = constitutive_partionedState0(g,i,e)%p ! ...microstructure - crystallite_subFp0(:,:,g,i,e) = crystallite_partionedFp0(:,:,g,i,e) ! ...plastic def grad - crystallite_subLp0(:,:,g,i,e) = crystallite_partionedLp0(:,:,g,i,e) ! ...plastic velocity grad - crystallite_subF0(:,:,g,i,e) = crystallite_partionedF0(:,:,g,i,e) ! ...def grad - crystallite_subTstar0_v(:,g,i,e) = crystallite_partionedTstar0_v(:,g,i,e) !...2nd PK stress + if (crystallite_requested(g,i,e)) then ! initialize restoration point of ... + crystallite_subTemperature0(g,i,e) = crystallite_partionedTemperature0(g,i,e) ! ...temperature + constitutive_subState0(g,i,e)%p = constitutive_partionedState0(g,i,e)%p ! ...microstructure + crystallite_subFp0(:,:,g,i,e) = crystallite_partionedFp0(:,:,g,i,e) ! ...plastic def grad + crystallite_subLp0(:,:,g,i,e) = crystallite_partionedLp0(:,:,g,i,e) ! ...plastic velocity grad + crystallite_subF0(:,:,g,i,e) = crystallite_partionedF0(:,:,g,i,e) ! ...def grad + crystallite_subTstar0_v(:,g,i,e) = crystallite_partionedTstar0_v(:,g,i,e) ! ...2nd PK stress crystallite_subFrac(g,i,e) = 0.0_pReal crystallite_subStep(g,i,e) = 1.0_pReal/subStepSizeCryst @@ -472,7 +463,6 @@ subroutine crystallite_stressAndItsTangent(updateJaco) enddo enddo !$OMPEND PARALLEL DO - ! --+>> crystallite loop <<+-- @@ -553,7 +543,7 @@ subroutine crystallite_stressAndItsTangent(updateJaco) !$OMP PARALLEL DO do e = FEsolving_execElem(1),FEsolving_execElem(2) ! iterate over elements to be processed myNgrains = homogenization_Ngrains(mesh_element(3,e)) - do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) ! iterate over IPs of this element to be processed + do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) ! iterate over IPs of this element to be processed do g = 1,myNgrains if (crystallite_todo(g,i,e)) then ! all undone crystallites call constitutive_microstructure(crystallite_Temperature(g,i,e), crystallite_Tstar_v(:,g,i,e), crystallite_Fe, & @@ -583,7 +573,7 @@ subroutine crystallite_stressAndItsTangent(updateJaco) !$OMP PARALLEL DO do e = FEsolving_execElem(1),FEsolving_execElem(2) ! iterate over elements to be processed myNgrains = homogenization_Ngrains(mesh_element(3,e)) - do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) ! iterate over IPs of this element to be processed + do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) ! iterate over IPs of this element to be processed do g = 1,myNgrains selectiveDebugger = (e == debug_e .and. i == debug_i .and. g == debug_g) if (crystallite_todo(g,i,e)) then ! all undone crystallites @@ -612,9 +602,9 @@ subroutine crystallite_stressAndItsTangent(updateJaco) ! to account for substepping within _integrateStress ! results in crystallite_Fp,.._Lp !$OMP PARALLEL DO - do e = FEsolving_execElem(1),FEsolving_execElem(2) ! iterate over elements to be processed + do e = FEsolving_execElem(1),FEsolving_execElem(2) ! iterate over elements to be processed myNgrains = homogenization_Ngrains(mesh_element(3,e)) - do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) ! iterate over IPs of this element to be processed + do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) ! iterate over IPs of this element to be processed do g = 1,myNgrains selectiveDebugger = (e == debug_e .and. i == debug_i .and. g == debug_g) if (crystallite_todo(g,i,e)) then ! all undone crystallites @@ -645,10 +635,10 @@ subroutine crystallite_stressAndItsTangent(updateJaco) myNgrains = homogenization_Ngrains(mesh_element(3,e)) do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) ! iterate over IPs of this element to be processed do g = 1,myNgrains - if (crystallite_todo(g,i,e)) then ! all undone crystallites + if (crystallite_todo(g,i,e)) then ! all undone crystallites constitutive_previousDotState2(g,i,e)%p = constitutive_previousDotState(g,i,e)%p constitutive_previousDotState(g,i,e)%p = constitutive_dotState(g,i,e)%p - constitutive_dotState(g,i,e)%p = 0.0_pReal ! zero out dotState + constitutive_dotState(g,i,e)%p = 0.0_pReal ! zero out dotState endif enddo; enddo; enddo !$OMPEND PARALLEL DO @@ -658,7 +648,7 @@ subroutine crystallite_stressAndItsTangent(updateJaco) !$OMP PARALLEL DO do e = FEsolving_execElem(1),FEsolving_execElem(2) ! iterate over elements to be processed myNgrains = homogenization_Ngrains(mesh_element(3,e)) - do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) ! iterate over IPs of this element to be processed + do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) ! iterate over IPs of this element to be processed do g = 1,myNgrains selectiveDebugger = (e == debug_e .and. i == debug_i .and. g == debug_g) if (crystallite_todo(g,i,e)) then ! all undone crystallites @@ -679,9 +669,9 @@ subroutine crystallite_stressAndItsTangent(updateJaco) !$OMPEND PARALLEL DO !$OMP PARALLEL DO - do e = FEsolving_execElem(1),FEsolving_execElem(2) ! iterate over elements to be processed + do e = FEsolving_execElem(1),FEsolving_execElem(2) ! iterate over elements to be processed myNgrains = homogenization_Ngrains(mesh_element(3,e)) - do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) ! iterate over IPs of this element to be processed + do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) ! iterate over IPs of this element to be processed do g = 1,myNgrains selectiveDebugger = (e == debug_e .and. i == debug_i .and. g == debug_g) if (crystallite_todo(g,i,e)) then ! all undone crystallites @@ -728,18 +718,18 @@ subroutine crystallite_stressAndItsTangent(updateJaco) !$OMPEND CRITICAL (write2out) endif - enddo ! crystallite convergence loop + enddo ! crystallite convergence loop NiterationCrystallite = NiterationCrystallite + 1 - enddo ! cutback loop + enddo ! cutback loop ! ------ check for non-converged crystallites ------ !$OMP PARALLEL DO - do e = FEsolving_execElem(1),FEsolving_execElem(2) ! iterate over elements to be processed + do e = FEsolving_execElem(1),FEsolving_execElem(2) ! iterate over elements to be processed myNgrains = homogenization_Ngrains(mesh_element(3,e)) - do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) ! iterate over IPs of this element to be processed + do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) ! iterate over IPs of this element to be processed do g = 1,myNgrains - if (.not. crystallite_converged(g,i,e)) then ! respond fully elastically (might be not required due to becoming terminally ill anyway) + if (.not. crystallite_converged(g,i,e)) then ! respond fully elastically (might be not required due to becoming terminally ill anyway) ! call IO_warning(600,e,i,g) invFp = math_inv3x3(crystallite_partionedFp0(:,:,g,i,e)) Fe_guess = math_mul33x33(crystallite_partionedF(:,:,g,i,e),invFp) @@ -861,11 +851,11 @@ subroutine crystallite_stressAndItsTangent(updateJaco) enddo ! perturbation direction select case(pert_method) case (1) - crystallite_dPdF(:,:,:,:,g,i,e) = dPdF_perturbation(:,:,:,:,1) + crystallite_dPdF(:,:,:,:,g,i,e) = dPdF_perturbation(:,:,:,:,1) case (2) - crystallite_dPdF(:,:,:,:,g,i,e) = dPdF_perturbation(:,:,:,:,2) + crystallite_dPdF(:,:,:,:,g,i,e) = dPdF_perturbation(:,:,:,:,2) case (3) - crystallite_dPdF(:,:,:,:,g,i,e) = 0.5_pReal*(dPdF_perturbation(:,:,:,:,1)+dPdF_perturbation(:,:,:,:,2)) + crystallite_dPdF(:,:,:,:,g,i,e) = 0.5_pReal*(dPdF_perturbation(:,:,:,:,1)+dPdF_perturbation(:,:,:,:,2)) end select else ! grain did not converge crystallite_dPdF(:,:,:,:,g,i,e) = crystallite_fallbackdPdF(:,:,:,:,g,i,e) ! use (elastic) fallback @@ -876,7 +866,7 @@ subroutine crystallite_stressAndItsTangent(updateJaco) enddo ! element loop !$OMPEND PARALLEL DO - elseif (any(.not. crystallite_localConstitution)) then ! if any nonlocal grain present, we have to do a full loop over all grains after each perturbance + elseif (any(.not. crystallite_localConstitution)) then ! if any nonlocal grain present, we have to do a full loop over all grains after each perturbance do e = FEsolving_execElem(1),FEsolving_execElem(2) myNgrains = homogenization_Ngrains(mesh_element(3,e)) diff --git a/code/debug.f90 b/code/debug.f90 index 39228116d..c7a8f1d51 100644 --- a/code/debug.f90 +++ b/code/debug.f90 @@ -41,8 +41,6 @@ subroutine debug_init() nHomog implicit none - if (.not. debugger) verboseDebugger = .false. - write(6,*) write(6,*) '<<<+- debug init -+>>>' write(6,*) '$Id$' diff --git a/code/homogenization.f90 b/code/homogenization.f90 index a68793c17..74b2ec5dd 100644 --- a/code/homogenization.f90 +++ b/code/homogenization.f90 @@ -263,6 +263,7 @@ subroutine materialpoint_stressAndItsTangent(& debug_i, & debug_MaterialpointLoopDistribution, & debug_MaterialpointStateLoopDistribution + use math, only: math_pDecomposition implicit none @@ -270,6 +271,7 @@ subroutine materialpoint_stressAndItsTangent(& logical, intent(in) :: updateJaco integer(pInt) NiterationHomog,NiterationMPstate integer(pInt) g,i,e,myNgrains + logical error ! ------ initialize to starting condition ------ @@ -283,6 +285,7 @@ subroutine materialpoint_stressAndItsTangent(& write (6,'(a,/,3(3(f12.7,x)/))') 'Lp0 of 1 1 1',crystallite_Lp0(1:3,:,1,1,1) endif + !$OMP PARALLEL DO do e = FEsolving_execElem(1),FEsolving_execElem(2) ! iterate over elements to be processed myNgrains = homogenization_Ngrains(mesh_element(3,e)) diff --git a/code/homogenization_RGC.f90 b/code/homogenization_RGC.f90 index 47c249239..42ee91a48 100644 --- a/code/homogenization_RGC.f90 +++ b/code/homogenization_RGC.f90 @@ -18,14 +18,15 @@ MODULE homogenization_RGC character (len=*), parameter :: homogenization_RGC_label = 'rgc' - integer(pInt), dimension(:), allocatable :: homogenization_RGC_sizeState, & - homogenization_RGC_sizePostResults + integer(pInt), dimension(:), allocatable :: homogenization_RGC_sizeState, & + homogenization_RGC_sizePostResults integer(pInt), dimension(:,:), allocatable,target :: homogenization_RGC_sizePostResult - integer(pInt), dimension(:,:), allocatable :: homogenization_RGC_Ngrains - real(pReal), dimension(:,:), allocatable :: homogenization_RGC_xiAlpha, & - homogenization_RGC_ciAlpha - real(pReal), dimension(:), allocatable :: homogenization_RGC_maxVol0, & - homogenization_RGC_vPower0 + integer(pInt), dimension(:,:), allocatable :: homogenization_RGC_Ngrains + real(pReal), dimension(:,:), allocatable :: homogenization_RGC_dAlpha, & + homogenization_RGC_angles + real(pReal), dimension(:,:,:,:), allocatable :: homogenization_RGC_orientation + real(pReal), dimension(:), allocatable :: homogenization_RGC_xiAlpha, & + homogenization_RGC_ciAlpha character(len=64), dimension(:,:), allocatable,target :: homogenization_RGC_output ! name of each post result output CONTAINS @@ -71,11 +72,15 @@ subroutine homogenization_RGC_init(& allocate(homogenization_RGC_sizePostResult(maxval(homogenization_Noutput), & maxNinstance)); homogenization_RGC_sizePostResult = 0_pInt allocate(homogenization_RGC_Ngrains(3,maxNinstance)); homogenization_RGC_Ngrains = 0_pInt - allocate(homogenization_RGC_ciAlpha(3,maxNinstance)); homogenization_RGC_ciAlpha = 0.0_pReal - allocate(homogenization_RGC_xiAlpha(3,maxNinstance)); homogenization_RGC_xiAlpha = 0.0_pReal - allocate(homogenization_RGC_maxVol0(maxNinstance)); homogenization_RGC_maxVol0 = 0.0_pReal - allocate(homogenization_RGC_vPower0(maxNinstance)); homogenization_RGC_vPower0 = 0.0_pReal + allocate(homogenization_RGC_ciAlpha(maxNinstance)); homogenization_RGC_ciAlpha = 0.0_pReal + allocate(homogenization_RGC_xiAlpha(maxNinstance)); homogenization_RGC_xiAlpha = 0.0_pReal + allocate(homogenization_RGC_dAlpha(3,maxNinstance)); homogenization_RGC_dAlpha = 0.0_pReal + allocate(homogenization_RGC_angles(3,maxNinstance)); homogenization_RGC_angles = 400.0_pReal allocate(homogenization_RGC_output(maxval(homogenization_Noutput),maxNinstance)); homogenization_RGC_output = '' + allocate(homogenization_RGC_orientation(3,3,mesh_maxNips,mesh_NcpElems)) + forall (i = 1:mesh_maxNips,e = 1:mesh_NcpElems) + homogenization_RGC_orientation(:,:,i,e) = math_I3 + end forall rewind(file) line = '' @@ -105,23 +110,55 @@ subroutine homogenization_RGC_init(& homogenization_RGC_Ngrains(1,i) = IO_intValue(line,positions,2) homogenization_RGC_Ngrains(2,i) = IO_intValue(line,positions,3) homogenization_RGC_Ngrains(3,i) = IO_intValue(line,positions,4) - case ('grainsizeparameter') - homogenization_RGC_xiAlpha(1,i) = IO_floatValue(line,positions,2) - homogenization_RGC_xiAlpha(2,i) = IO_floatValue(line,positions,3) - homogenization_RGC_xiAlpha(3,i) = IO_floatValue(line,positions,4) + case ('scalingparameter') + homogenization_RGC_xiAlpha(i) = IO_floatValue(line,positions,2) case ('overproportionality') - homogenization_RGC_ciAlpha(1,i) = IO_floatValue(line,positions,2) - homogenization_RGC_ciAlpha(2,i) = IO_floatValue(line,positions,3) - homogenization_RGC_ciAlpha(3,i) = IO_floatValue(line,positions,4) - case ('maxvoldiscrepancy') - homogenization_RGC_maxVol0(i) = IO_floatValue(line,positions,2) - case ('discrepancypower') - homogenization_RGC_vPower0(i) = IO_floatValue(line,positions,2) + homogenization_RGC_ciAlpha(i) = IO_floatValue(line,positions,2) + case ('grainsize') + homogenization_RGC_dAlpha(1,i) = IO_floatValue(line,positions,2) + homogenization_RGC_dAlpha(2,i) = IO_floatValue(line,positions,3) + homogenization_RGC_dAlpha(3,i) = IO_floatValue(line,positions,4) + case ('clusterorientation') + homogenization_RGC_angles(1,i) = IO_floatValue(line,positions,2) + homogenization_RGC_angles(2,i) = IO_floatValue(line,positions,3) + homogenization_RGC_angles(3,i) = IO_floatValue(line,positions,4) end select endif enddo -100 do i = 1,maxNinstance ! sanity checks +!*** assigning cluster orientations + do e = 1,mesh_NcpElems + if (homogenization_type(mesh_element(3,e)) == homogenization_RGC_label) then + myInstance = homogenization_typeInstance(mesh_element(3,e)) + if (all (homogenization_RGC_angles(:,myInstance) >= 399.9_pReal)) then + homogenization_RGC_orientation(:,:,1,e) = math_EulerToR(math_sampleRandomOri()) + do i = 1,FE_Nips(mesh_element(2,e)) + if (microstructure_elemhomo(mesh_element(4,e))) then + homogenization_RGC_orientation(:,:,i,e) = homogenization_RGC_orientation(:,:,1,e) + else + homogenization_RGC_orientation(:,:,i,e) = math_EulerToR(math_sampleRandomOri()) + endif + enddo + else + do i = 1,FE_Nips(mesh_element(2,e)) + homogenization_RGC_orientation(:,:,i,e) = math_EulerToR(homogenization_RGC_angles(:,myInstance)*inRad) + enddo + endif + endif + enddo + +100 do i = 1,maxNinstance ! sanity checks + write(6,*) + write(6,*) '<<<+- homogenization_RGC init -+>>>' + write(6,*) '$Id$' + write(6,*) + write(6,'(a15,x,i4)') 'instance: ', i + write(6,*) + write(6,'(a25,3(x,i8))') 'cluster size: ',(homogenization_RGC_Ngrains(j,i),j=1,3) + write(6,'(a25,x,e10.3)') 'scaling parameter: ', homogenization_RGC_xiAlpha(i) + write(6,'(a25,x,e10.3)') 'over-proportionality: ', homogenization_RGC_ciAlpha(i) + write(6,'(a25,3(x,e10.3))') 'grain size: ',(homogenization_RGC_dAlpha(j,i),j=1,3) + write(6,'(a25,3(x,e10.3))') 'cluster orientation: ',(homogenization_RGC_angles(j,i),j=1,3) enddo do i = 1,maxNinstance @@ -130,13 +167,17 @@ subroutine homogenization_RGC_init(& case('constitutivework') mySize = 1 case('magnitudemismatch') - mySize = 1 + mySize = 3 case('penaltyenergy') mySize = 1 case('volumediscrepancy') mySize = 1 case default mySize = 0 + case('averagerelaxrate') + mySize = 1 + case('maximumrelaxrate') + mySize = 1 end select if (mySize > 0_pInt) then ! any meaningful output found @@ -146,13 +187,12 @@ subroutine homogenization_RGC_init(& endif enddo - homogenization_RGC_sizeState(i) & = 3*(homogenization_RGC_Ngrains(1,i)-1)*homogenization_RGC_Ngrains(2,i)*homogenization_RGC_Ngrains(3,i) & + 3*homogenization_RGC_Ngrains(1,i)*(homogenization_RGC_Ngrains(2,i)-1)*homogenization_RGC_Ngrains(3,i) & + 3*homogenization_RGC_Ngrains(1,i)*homogenization_RGC_Ngrains(2,i)*(homogenization_RGC_Ngrains(3,i)-1) & - + 4_pInt ! (1) Average constitutive work, (2) Overall mismatch, (3) Average penalty energy, - ! (4) Volume discrepancy + + 8_pInt ! (1) Average constitutive work, (2-4) Overall mismatch, (5) Average penalty energy, + ! (6) Volume discrepancy, (7) Avg relaxation rate component, (8) Max relaxation rate component enddo return @@ -235,7 +275,7 @@ subroutine homogenization_RGC_partitionDeformation(& do iFace = 1,nFace call homogenization_RGC_getInterface(intFace,iFace,iGrain3) call homogenization_RGC_relaxationVector(aVect,intFace,state,homID) - call homogenization_RGC_interfaceNormal(nVect,intFace) + call homogenization_RGC_interfaceNormal(nVect,intFace,ip,el) forall (i=1:3,j=1:3) & F(i,j,iGrain) = F(i,j,iGrain) + aVect(i)*nVect(j) ! effective relaxations enddo @@ -279,7 +319,7 @@ function homogenization_RGC_updateState(& use material, only: homogenization_maxNgrains,homogenization_typeInstance, & homogenization_Ngrains use numerics, only: absTol_RGC,relTol_RGC,absMax_RGC,relMax_RGC,pPert_RGC, & - maxdRelax_RGC,ratePower_RGC,viscModus_RGC + maxdRelax_RGC,viscPower_RGC,viscModus_RGC,refRelaxRate_RGC use FEsolving, only: theInc,cycleCounter,theTime implicit none @@ -299,7 +339,7 @@ function homogenization_RGC_updateState(& integer(pInt), dimension (2) :: residLoc integer(pInt) homID,i1,i2,i3,iNum,i,j,nIntFaceTot,iGrN,iGrP,iMun,iFace,k,l,ival,ipert,iGrain,nGrain real(pReal), dimension (3,3,homogenization_maxNgrains) :: R,pF,pR,D,pD - real(pReal), dimension (homogenization_maxNgrains) :: NN,pNN + real(pReal), dimension (3,homogenization_maxNgrains) :: NN,pNN real(pReal), dimension (3) :: normP,normN,mornP,mornN real(pReal) residMax,stresMax,constitutiveWork,penaltyEnergy,volDiscrep,penDiscrep logical error,RGCdebug,RGCdebugJacobi,RGCcheck @@ -336,7 +376,7 @@ function homogenization_RGC_updateState(& endif !* Stress-like penalty related to mismatch or incompatibility at interfaces - call homogenization_RGC_stressPenalty(R,NN,F,ip,el,homID) + call homogenization_RGC_stressPenalty(R,NN,avgF,F,ip,el,homID) !* Stress-like penalty related to overall volume discrepancy call homogenization_RGC_volumePenalty(D,volDiscrep,F,avgF,ip,el,homID) @@ -344,7 +384,7 @@ function homogenization_RGC_updateState(& !* Debugging the mismatch, stress and penalties of grains if (RGCdebug) then do iGrain = 1,nGrain - write(6,'(x,a30,x,i3,x,a4,x,e14.8)')'Mismatch magnitude of grain(',iGrain,') :',NN(iGrain) + write(6,'(x,a30,x,i3,x,a4,3(x,e14.8))')'Mismatch magnitude of grain(',iGrain,') :',NN(1,iGrain),NN(2,iGrain),NN(3,iGrain) write(6,*)' ' write(6,'(x,a30,x,i3)')'Stress and penalties of grain: ',iGrain do i = 1,3 @@ -364,15 +404,15 @@ function homogenization_RGC_updateState(& iGr3N = faceID(2:4) call homogenization_RGC_grain3to1(iGrN,iGr3N,homID) call homogenization_RGC_getInterface(intFaceN,2*faceID(1),iGr3N) - call homogenization_RGC_interfaceNormal(normN,intFaceN) ! get the interface normal + call homogenization_RGC_interfaceNormal(normN,intFaceN,ip,el) ! get the interface normal !* Identify the right/up/front grain (+|P) iGr3P = iGr3N iGr3P(faceID(1)) = iGr3N(faceID(1))+1 call homogenization_RGC_grain3to1(iGrP,iGr3P,homID) call homogenization_RGC_getInterface(intFaceP,2*faceID(1)-1,iGr3P) - call homogenization_RGC_interfaceNormal(normP,intFaceP) ! get the interface normal + call homogenization_RGC_interfaceNormal(normP,intFaceP,ip,el) ! get the interface normal do i = 1,3 ! compute the traction balance at the interface - tract(iNum,i) = sign(viscModus_RGC*(abs(drelax(i+3*(iNum-1))/dt))**ratePower_RGC, & + tract(iNum,i) = sign(viscModus_RGC*(abs(drelax(i+3*(iNum-1)))/(refRelaxRate_RGC*dt))**viscPower_RGC, & drelax(i+3*(iNum-1))) do j = 1,3 tract(iNum,i) = tract(iNum,i) + (P(i,j,iGrP) + R(i,j,iGrP) + D(i,j,iGrP))*normP(j) & @@ -416,7 +456,7 @@ function homogenization_RGC_updateState(& !* Then compute/update the state for postResult, i.e., ... !* ... all energy densities computed by time-integration constitutiveWork = state%p(3*nIntFaceTot+1) - penaltyEnergy = state%p(3*nIntFaceTot+3) + penaltyEnergy = state%p(3*nIntFaceTot+5) do iGrain = 1,homogenization_Ngrains(mesh_element(3,el)) do i = 1,3 do j = 1,3 @@ -428,16 +468,23 @@ function homogenization_RGC_updateState(& !* ... the bulk mechanical/constitutive work state%p(3*nIntFaceTot+1) = constitutiveWork !* ... the overall mismatch - state%p(3*nIntFaceTot+2) = sum(NN)/dble(nGrain) - state%p(3*nIntFaceTot+3) = penaltyEnergy + state%p(3*nIntFaceTot+2) = sum(NN(1,:))/dble(nGrain) + state%p(3*nIntFaceTot+3) = sum(NN(2,:))/dble(nGrain) + state%p(3*nIntFaceTot+4) = sum(NN(3,:))/dble(nGrain) + state%p(3*nIntFaceTot+5) = penaltyEnergy !* ... the volume discrepancy - state%p(3*nIntFaceTot+4) = volDiscrep + state%p(3*nIntFaceTot+6) = volDiscrep + state%p(3*nIntFaceTot+7) = sum(abs(drelax))/dt/dble(3*nIntFaceTot) + state%p(3*nIntFaceTot+8) = maxval(abs(drelax))/dt if (el == 1 .and. ip == 1) then write(6,'(x,a30,x,e14.8)')'Constitutive work: ',constitutiveWork - write(6,'(x,a30,x,e14.8)')'Magnitude mismatch: ',sum(NN) + write(6,'(x,a30,3(x,e14.8))')'Magnitude mismatch: ',sum(NN(1,:))/dble(nGrain),sum(NN(2,:))/dble(nGrain),sum(NN(3,:))/dble(nGrain) write(6,'(x,a30,x,e14.8)')'Penalty energy: ',penaltyEnergy write(6,'(x,a30,x,e14.8)')'Volume discrepancy: ',volDiscrep write(6,*)'' + write(6,'(x,a30,x,e14.8)')'Maximum relaxation rate: ',maxval(abs(drelax))/dt + write(6,'(x,a30,x,e14.8)')'Average relaxation rate: ',sum(abs(drelax))/dt/dble(3*nIntFaceTot) + write(6,*)'' call flush(6) endif deallocate(tract,resid,relax,drelax) @@ -446,8 +493,8 @@ function homogenization_RGC_updateState(& elseif (residMax > relMax_RGC*stresMax .or. residMax > absMax_RGC) then !* Try to restart when residual blows up exceeding maximum bound homogenization_RGC_updateState = (/.true.,.false./) ! ... with direct cut-back - write(6,'(x,a,x,i3,x,a,x,i3,x,a)')'RGC_updateState: ip',ip,'| el',el,'enforces cutback' - write(6,'(x,a,x,e14.8,x,a,x,e14.8)')'due to high residual =',residMax,'for stress =',stresMax +! write(6,'(x,a,x,i3,x,a,x,i3,x,a)')'RGC_updateState: ip',ip,'| el',el,'enforces cutback' +! write(6,'(x,a,x,e14.8,x,a,x,e14.8)')'due to high residual =',residMax,'for stress =',stresMax ! state%p(1:3*nIntFaceTot) = 0.0_pReal ! ... with full Taylor assumption ! write(6,'(x,a,x,i3,x,a,x,i3,x,a)')'RGC_updateState: ip',ip,'| el',el,'relaxation vectors reset to zero' if (RGCcheck) then @@ -475,10 +522,10 @@ function homogenization_RGC_updateState(& iGr3N = faceID(2:4) call homogenization_RGC_grain3to1(iGrN,iGr3N,homID) call homogenization_RGC_getInterface(intFaceN,2*faceID(1),iGr3N) - call homogenization_RGC_interfaceNormal(normN,intFaceN) ! get the interface normal + call homogenization_RGC_interfaceNormal(normN,intFaceN,ip,el) ! get the interface normal do iFace = 1,nFace call homogenization_RGC_getInterface(intFaceN,iFace,iGr3N) - call homogenization_RGC_interfaceNormal(mornN,intFaceN) ! get influencing interfaces normal + call homogenization_RGC_interfaceNormal(mornN,intFaceN,ip,el) ! get influencing interfaces normal call homogenization_RGC_interface4to1(iMun,intFaceN,homID) if (iMun .gt. 0) then ! get the tangent forall(i=1:3,j=1:3,k=1:3,l=1:3) & @@ -490,10 +537,10 @@ function homogenization_RGC_updateState(& iGr3P(faceID(1)) = iGr3N(faceID(1))+1 call homogenization_RGC_grain3to1(iGrP,iGr3P,homID) call homogenization_RGC_getInterface(intFaceP,2*faceID(1)-1,iGr3P) - call homogenization_RGC_interfaceNormal(normP,intFaceP) ! get the interface normal + call homogenization_RGC_interfaceNormal(normP,intFaceP,ip,el) ! get the interface normal do iFace = 1,nFace call homogenization_RGC_getInterface(intFaceP,iFace,iGr3P) - call homogenization_RGC_interfaceNormal(mornP,intFaceP) ! get influencing interfaces normal + call homogenization_RGC_interfaceNormal(mornP,intFaceP,ip,el) ! get influencing interfaces normal call homogenization_RGC_interface4to1(iMun,intFaceP,homID) if (iMun .gt. 0) then ! get the tangent forall(i=1:3,j=1:3,k=1:3,l=1:3) & @@ -519,8 +566,8 @@ function homogenization_RGC_updateState(& p_relax = relax p_relax(ipert) = relax(ipert) + pPert_RGC ! perturb the relaxation vector state%p(1:3*nIntFaceTot) = p_relax - call homogenization_RGC_grainDeformation(pF,F0,avgF,state,el) - call homogenization_RGC_stressPenalty(pR,pNN,pF,ip,el,homID) + call homogenization_RGC_grainDeformation(pF,F0,avgF,state,ip,el) + call homogenization_RGC_stressPenalty(pR,pNN,avgF,pF,ip,el,homID) call homogenization_RGC_volumePenalty(pD,volDiscrep,pF,avgF,ip,el,homID) p_resid = 0.0_pReal do iNum = 1,nIntFaceTot @@ -529,13 +576,13 @@ function homogenization_RGC_updateState(& iGr3N = faceID(2:4) call homogenization_RGC_grain3to1(iGrN,iGr3N,homID) call homogenization_RGC_getInterface(intFaceN,2*faceID(1),iGr3N) - call homogenization_RGC_interfaceNormal(normN,intFaceN) ! get the corresponding normal + call homogenization_RGC_interfaceNormal(normN,intFaceN,ip,el) ! get the corresponding normal !* Identify the right/up/front grain (+|P) iGr3P = iGr3N iGr3P(faceID(1)) = iGr3N(faceID(1))+1 call homogenization_RGC_grain3to1(iGrP,iGr3P,homID) call homogenization_RGC_getInterface(intFaceP,2*faceID(1)-1,iGr3P) - call homogenization_RGC_interfaceNormal(normP,intFaceP) ! get the corresponding normal + call homogenization_RGC_interfaceNormal(normP,intFaceP,ip,el) ! get the corresponding normal !* Compute the perturbed traction at interface do i = 1,3 do j = 1,3 @@ -561,7 +608,8 @@ function homogenization_RGC_updateState(& !* Construct the Jacobian matrix of the numerical viscosity tangent allocate(rmatrix(3*nIntFaceTot,3*nIntFaceTot)); rmatrix = 0.0_pReal forall (i=1:3*nIntFaceTot) & - rmatrix(i,i) = viscModus_RGC*ratePower_RGC/dt*(abs(drelax(i)/dt))**(ratePower_RGC - 1.0_pReal) + rmatrix(i,i) = viscModus_RGC*viscPower_RGC/(refRelaxRate_RGC*dt)* & + (abs(drelax(i))/(refRelaxRate_RGC*dt))**(viscPower_RGC - 1.0_pReal) !* Debugging the global Jacobian matrix of numerical viscosity tangent if (RGCdebugJacobi) then write(6,'(x,a30)')'Jacobian matrix of penalty' @@ -609,9 +657,9 @@ function homogenization_RGC_updateState(& ! state%p(1:3*nIntFaceTot) = 0.0_pReal ! write(6,'(x,a,x,i3,x,a,x,i3,x,a)')'RGC_updateState: ip',ip,'| el',el,'relaxation vectors reset to zero' homogenization_RGC_updateState = (/.true.,.false./) - write(6,'(x,a,x,i3,x,a,x,i3,x,a)')'RGC_updateState: ip',ip,'| el',el,'enforces cutback' - write(6,'(x,a,x,e14.8)')'due to large relaxation change =',maxval(abs(drelax)) - call flush(6) +! write(6,'(x,a,x,i3,x,a,x,i3,x,a)')'RGC_updateState: ip',ip,'| el',el,'enforces cutback' +! write(6,'(x,a,x,e14.8)')'due to large relaxation change =',maxval(abs(drelax)) +! call flush(6) endif !* Debugging the return state if (RGCdebugJacobi) then @@ -744,13 +792,21 @@ pure function homogenization_RGC_postResults(& homogenization_RGC_postResults(c+1) = state%p(3*nIntFaceTot+1) c = c + 1 case('magnitudemismatch') - homogenization_RGC_postResults(c+1) = state%p(3*nIntFaceTot+3) - c = c + 1 - case('penaltyenergy') homogenization_RGC_postResults(c+1) = state%p(3*nIntFaceTot+2) + homogenization_RGC_postResults(c+2) = state%p(3*nIntFaceTot+3) + homogenization_RGC_postResults(c+3) = state%p(3*nIntFaceTot+4) + c = c + 3 + case('penaltyenergy') + homogenization_RGC_postResults(c+1) = state%p(3*nIntFaceTot+5) c = c + 1 case('volumediscrepancy') - homogenization_RGC_postResults(c+1) = state%p(3*nIntFaceTot+4) + homogenization_RGC_postResults(c+1) = state%p(3*nIntFaceTot+6) + c = c + 1 + case('averagerelaxrate') + homogenization_RGC_postResults(c+1) = state%p(3*nIntFaceTot+7) + c = c + 1 + case('maximumrelaxrate') + homogenization_RGC_postResults(c+1) = state%p(3*nIntFaceTot+8) c = c + 1 end select enddo @@ -766,6 +822,7 @@ subroutine homogenization_RGC_stressPenalty(& rPen, & ! stress-like penalty nMis, & ! total amount of mismatch ! + avgF, & ! initial effective stretch tensor fDef, & ! deformation gradients ip, & ! integration point el, & ! element @@ -774,22 +831,23 @@ subroutine homogenization_RGC_stressPenalty(& use prec, only: pReal,pInt,p_vec use mesh, only: mesh_element use constitutive, only: constitutive_homogenizedC - use math, only: math_civita + use math, only: math_civita,math_invert3x3 use material, only: homogenization_maxNgrains,homogenization_Ngrains use numerics, only: xSmoo_RGC implicit none !* Definition of variables real(pReal), dimension (3,3,homogenization_maxNgrains), intent(out) :: rPen - real(pReal), dimension (homogenization_maxNgrains), intent(out) :: nMis + real(pReal), dimension (3,homogenization_maxNgrains), intent(out) :: nMis real(pReal), dimension (3,3,homogenization_maxNgrains), intent(in) :: fDef + real(pReal), dimension (3,3), intent(in) :: avgF integer(pInt), intent(in) :: ip,el integer(pInt), dimension (4) :: intFace integer(pInt), dimension (3) :: iGrain3,iGNghb3,nGDim - real(pReal), dimension (3,3) :: gDef,nDef - real(pReal), dimension (3) :: nVect - integer(pInt) homID,iGrain,iGNghb,iFace,i,j,k,l - real(pReal) muGrain,muGNghb,nDefNorm + real(pReal), dimension (3,3) :: gDef,nDef,avgC + real(pReal), dimension (3) :: nVect,surfCorr + integer(pInt) homID,iGrain,iGNghb,iFace,i,j,k,l,m + real(pReal) muGrain,muGNghb,nDefNorm,xiAlpha,ciAlpha,bgGrain,bgGNghb,detF ! integer(pInt), parameter :: nFace = 6 real(pReal), parameter :: nDefToler = 1.0e-10 @@ -798,13 +856,23 @@ subroutine homogenization_RGC_stressPenalty(& rPen = 0.0_pReal nMis = 0.0_pReal + +!* Get the penalty correction factor representing surface evolution + call homogenization_RGC_surfaceCorrection(surfCorr,avgF,ip,el) + +!* Debugging the surface correction factor +! if (ip == 1 .and. el == 1) then +! write(6,'(x,a20,2(x,i3))')'Correction factor: ',ip,el +! write(6,'(x,3(e10.4,x))')(surfCorr(i), i = 1,3) +! endif + do iGrain = 1,homogenization_Ngrains(mesh_element(3,el)) - call homogenization_RGC_equivalentModuli(muGrain,constitutive_homogenizedC(iGrain,ip,el)) + call homogenization_RGC_equivalentModuli(muGrain,bgGrain,iGrain,ip,el) call homogenization_RGC_grain1to3(iGrain3,iGrain,homID) !* Compute the mismatch tensor at all six interfaces do iFace = 1,nFace call homogenization_RGC_getInterface(intFace,iFace,iGrain3) - call homogenization_RGC_interfaceNormal(nVect,intFace) ! get the interface normal + call homogenization_RGC_interfaceNormal(nVect,intFace,ip,el) ! get the interface normal iGNghb3 = iGrain3 ! identify the grain neighbor iGNghb3(abs(intFace(1))) = iGNghb3(abs(intFace(1))) + int(dble(intFace(1))/dble(abs(intFace(1)))) !* The grain periodicity along e1 @@ -817,7 +885,7 @@ subroutine homogenization_RGC_stressPenalty(& if (iGNghb3(3) < 1) iGNghb3(3) = nGDim(3) if (iGNghb3(3) > nGDim(3)) iGNghb3(3) = 1 call homogenization_RGC_grain3to1(iGNghb,iGNghb3,homID) ! get the grain neighbor - call homogenization_RGC_equivalentModuli(muGNghb,constitutive_homogenizedC(iGNghb,ip,el)) + call homogenization_RGC_equivalentModuli(muGNghb,bgGNghb,iGNghb,ip,el) gDef = 0.5_pReal*(fDef(:,:,iGNghb) - fDef(:,:,iGrain)) ! difference in F with the neighbor nDefNorm = 0.0_pReal nDef = 0.0_pReal @@ -845,8 +913,9 @@ subroutine homogenization_RGC_stressPenalty(& do j = 1,3 do k = 1,3 do l = 1,3 - rPen(i,j,iGrain) = rPen(i,j,iGrain) + 0.5_pReal*(muGrain + muGNghb)/homogenization_RGC_xiAlpha(abs(intFace(1)),homID) & - *cosh(homogenization_RGC_ciAlpha(abs(intFace(1)),homID)*nDefNorm) & + rPen(i,j,iGrain) = rPen(i,j,iGrain) + 0.5_pReal*(muGrain*bgGrain + muGNghb*bgGNghb)*homogenization_RGC_xiAlpha(homID) & + *surfCorr(abs(intFace(1)))/homogenization_RGC_dAlpha(abs(intFace(1)),homID) & + *cosh(homogenization_RGC_ciAlpha(homID)*nDefNorm) & *0.5_pReal*nVect(l)*nDef(i,k)/nDefNorm*math_civita(k,l,j) & *tanh(nDefNorm/xSmoo_RGC) enddo @@ -854,7 +923,7 @@ subroutine homogenization_RGC_stressPenalty(& enddo enddo !* Total amount of mismatch experienced by the grain (at all six interfaces) - nMis(iGrain) = nMis(iGrain) + nDefNorm + nMis(abs(intFace(1)),iGrain) = nMis(abs(intFace(1)),iGrain) + nDefNorm enddo !* Debugging the stress-like penalty ! if (ip == 1 .and. el == 1) then @@ -884,7 +953,6 @@ subroutine homogenization_RGC_volumePenalty(& ) use prec, only: pReal,pInt,p_vec use mesh, only: mesh_element - use constitutive, only: constitutive_homogenizedC use math, only: math_det3x3,math_inv3x3 use material, only: homogenization_maxNgrains,homogenization_Ngrains use numerics, only: maxVolDiscr_RGC,volDiscrMod_RGC,volDiscrPow_RGC @@ -903,18 +971,18 @@ subroutine homogenization_RGC_volumePenalty(& nGrain = homogenization_Ngrains(mesh_element(3,el)) !* Compute volumes of grain and the effective volume and the total volume discrepancy - vDiscrep = math_det3x3(fAvg) + vDiscrep = math_det3x3(fAvg) do iGrain = 1,nGrain gVol(iGrain) = math_det3x3(fDef(:,:,iGrain)) vDiscrep = vDiscrep - gVol(iGrain)/dble(nGrain) enddo !* Calculate the stress and penalty due to volume discrepancy - vPen = 0.0_pReal + vPen = 0.0_pReal do iGrain = 1,nGrain - vPen(:,:,iGrain) = -1.0_pReal*sign(volDiscrMod_RGC*volDiscrPow_RGC/maxVolDiscr_RGC* & - (abs(vDiscrep)/maxVolDiscr_RGC)**(volDiscrPow_RGC - 1.0)*gVol(iGrain)/dble(nGrain),vDiscrep)* & - transpose(math_inv3x3(fDef(:,:,iGrain))) + vPen(:,:,iGrain) = -1.0_pReal/dble(nGrain)*volDiscrMod_RGC*volDiscrPow_RGC/maxVolDiscr_RGC* & + sign((abs(vDiscrep)/maxVolDiscr_RGC)**(volDiscrPow_RGC - 1.0),vDiscrep)* & + gVol(iGrain)*transpose(math_inv3x3(fDef(:,:,iGrain))) !* Debugging the stress-like penalty of volume discrepancy ! if (ip == 1 .and. el == 1) then @@ -929,25 +997,79 @@ subroutine homogenization_RGC_volumePenalty(& endsubroutine +!******************************************************************** +! subroutine to compute the correction factor due to surface area evolution +!******************************************************************** +subroutine homogenization_RGC_surfaceCorrection(& + vSurf, & ! surface correction factor +! + avgF, & ! average deformation gradient + ip, & ! my IP + el & ! my element + ) + + use prec, only: pReal,pInt,p_vec + use math, only: math_invert3x3,math_mul33x33 + implicit none + +!* Definition of variables + real(pReal), dimension(3,3), intent(in) :: avgF + real(pReal), dimension(3), intent(out) :: vSurf + integer(pInt), intent(in) :: ip,el + real(pReal), dimension(3,3) :: invC,avgC + real(pReal), dimension(3) :: nVect + real(pReal) detF + integer(pInt), dimension(4) :: intFace + integer(pInt) i,j,iBase + logical error + +!* Compute the correction factor accouted for surface evolution (area change) + avgC = 0.0_pReal + avgC = math_mul33x33(transpose(avgF),avgF) + invC = 0.0_pReal + call math_invert3x3(avgC,invC,detF,error) + vSurf = 0.0_pReal + do iBase = 1,3 + intFace = (/iBase,1_pInt,1_pInt,1_pInt/) + call homogenization_RGC_interfaceNormal(nVect,intFace,ip,el) + do i = 1,3 + do j = 1,3 + vSurf(iBase) = vSurf(iBase) + invC(i,j)*nVect(i)*nVect(j) + enddo + enddo + vSurf(iBase) = sqrt(vSurf(iBase))*detF + enddo + + return + +endsubroutine + !******************************************************************** ! subroutine to compute the equivalent shear and bulk moduli from the elasticity tensor !******************************************************************** subroutine homogenization_RGC_equivalentModuli(& shearMod, & ! equivalent (isotropic) shear modulus + vBurgers, & ! length of burgers vector ! - elasTens & ! elasticity tensor in Mandel notation + grainID, & ! grain ID + ip, & ! IP number + el & ! element number ) use prec, only: pReal,pInt,p_vec + use constitutive, only: constitutive_homogenizedC,constitutive_averageBurgers use mesh, only: mesh_element,mesh_NcpElems,mesh_maxNips use material, only: homogenization_typeInstance implicit none !* Definition of variables - real(pReal), dimension (6,6), intent(in) :: elasTens - real(pReal), intent(out) :: shearMod + integer(pInt), intent(in) :: grainID,ip,el + real(pReal), dimension (6,6) :: elasTens + real(pReal), intent(out) :: shearMod,vBurgers real(pReal) cEquiv_11,cEquiv_12,cEquiv_44 + elasTens = constitutive_homogenizedC(grainID,ip,el) + !* Compute the equivalent shear modulus using Turterltaub and Suiker, JMPS (2005) cEquiv_11 = (elasTens(1,1) + elasTens(2,2) + elasTens(3,3))/3.0_pReal cEquiv_12 = (elasTens(1,2) + elasTens(2,3) + elasTens(3,1) + & @@ -955,6 +1077,8 @@ subroutine homogenization_RGC_equivalentModuli(& cEquiv_44 = (elasTens(4,4) + elasTens(5,5) + elasTens(6,6))/3.0_pReal shearMod = 0.2_pReal*(cEquiv_11 - cEquiv_12) + 0.6_pReal*cEquiv_44 + vBurgers = constitutive_averageBurgers(grainID,ip,el) + return endsubroutine @@ -998,23 +1122,36 @@ endsubroutine subroutine homogenization_RGC_interfaceNormal(& nVect, & ! interface normal ! - intFace & ! interface ID in 4D array (normal and position) + intFace, & ! interface ID in 4D array (normal and position) + ip, & ! my IP + el & ! my element ) use prec, only: pReal,pInt,p_vec + use math, only: math_mul33x3 use mesh, only: mesh_element,mesh_NcpElems,mesh_maxNips implicit none !* Definition of variables real(pReal), dimension (3), intent(out) :: nVect integer(pInt), dimension (4), intent(in) :: intFace - integer(pInt) nPos + integer(pInt), intent(in) :: ip,el + integer(pInt) nPos,i !* Get the normal of the interface, identified from the value of intFace(1) nVect = 0.0_pReal nPos = abs(intFace(1)) nVect(nPos) = intFace(1)/abs(intFace(1)) + nVect = math_mul33x3(homogenization_RGC_orientation(:,:,ip,el),nVect) + +! if (ip == 1 .and. el == 1) then +! write(6,'(x,a32,3(x,i3))')'Interface normal: ',intFace(1) +! write(6,'(x,3(e14.8,x))')(nVect(i), i = 1,3) +! write(6,*)' ' +! call flush(6) +! endif + return endsubroutine @@ -1211,6 +1348,7 @@ subroutine homogenization_RGC_grainDeformation(& F0, & ! initial partioned def grad per grain avgF, & ! my average def grad state, & ! my state + ip, & ! my IP el & ! my element ) use prec, only: pReal,pInt,p_vec @@ -1223,7 +1361,7 @@ subroutine homogenization_RGC_grainDeformation(& real(pReal), dimension (3,3,homogenization_maxNgrains), intent(in) :: F0 real(pReal), dimension (3,3), intent(in) :: avgF type(p_vec), intent(in) :: state - integer(pInt), intent(in) :: el + integer(pInt), intent(in) :: el,ip ! real(pReal), dimension (3) :: aVect,nVect integer(pInt), dimension (4) :: intFace @@ -1240,7 +1378,7 @@ subroutine homogenization_RGC_grainDeformation(& do iFace = 1,nFace call homogenization_RGC_getInterface(intFace,iFace,iGrain3) call homogenization_RGC_relaxationVector(aVect,intFace,state,homID) - call homogenization_RGC_interfaceNormal(nVect,intFace) + call homogenization_RGC_interfaceNormal(nVect,intFace,ip,el) forall (i=1:3,j=1:3) & F(i,j,iGrain) = F(i,j,iGrain) + aVect(i)*nVect(j) ! effective relaxations enddo diff --git a/code/material.f90 b/code/material.f90 index f7a1a23af..d841001f0 100644 --- a/code/material.f90 +++ b/code/material.f90 @@ -535,6 +535,7 @@ subroutine material_populateGrains() use math, only: math_sampleRandomOri, math_sampleGaussOri, math_sampleFiberOri, math_symmetricEulers, inDeg use mesh, only: mesh_element, mesh_maxNips, mesh_NcpElems, mesh_ipVolume, FE_Nips use IO, only: IO_error, IO_hybridIA + use FEsolving, only: FEsolving_execIP implicit none integer(pInt), dimension (:,:), allocatable :: Ngrains @@ -711,6 +712,7 @@ subroutine material_populateGrains() material_phase(g,i,e) = phaseOfGrain(grain+g) material_EulerAngles(:,g,i,e) = orientationOfGrain(:,grain+g) end forall + FEsolving_execIP(2,e) = 1_pInt ! restrict calculation to first IP only, since all other results are to be copied from this grain = grain + dGrains ! wind forward by NgrainsPerIP else forall (i = 1:FE_Nips(mesh_element(2,e)), g = 1:dGrains) ! loop over IPs and grains diff --git a/code/math.f90 b/code/math.f90 index 0276cf4f0..9065d5560 100644 --- a/code/math.f90 +++ b/code/math.f90 @@ -176,8 +176,6 @@ subroutine math_misorientation(dQ, Q1, Q2, symmetryType) return endif - allocate(mySymOperations(4,NsymOperations(symmetryType))) - mySymOperations = symOperations(:,sum(NsymOperations(1:symmetryType-1))+1:sum(NsymOperations(1:symmetryType))) ! choose symmetry operations according to crystal symmetry dQ(1) = -1.0_pReal ! start with maximum misorientation angle do s = 1,NsymOperations(symmetryType) ! loop ver symmetry operations @@ -999,6 +997,30 @@ pure function math_transpose3x3(A) ENDFUNCTION +!******************************************************************** +! equivalent scalar quantity of a full strain tensor +!******************************************************************** + PURE FUNCTION math_equivStrain33(m) + + use prec, only: pReal,pInt + implicit none + + real(pReal), dimension(3,3), intent(in) :: m + real(pReal) math_equivStrain33,e11,e22,e33,s12,s23,s31 + + e11 = (2.0_pReal*m(1,1)-m(2,2)-m(3,3))/3.0_pReal + e22 = (2.0_pReal*m(2,2)-m(3,3)-m(1,1))/3.0_pReal + e33 = (2.0_pReal*m(3,3)-m(1,1)-m(2,2))/3.0_pReal + s12 = 2.0_pReal*m(1,2) + s23 = 2.0_pReal*m(2,3) + s31 = 2.0_pReal*m(3,1) + + math_equivStrain33 = 2.0_pReal*(1.50_pReal*(e11**2.0_pReal+e22**2.0_pReal+e33**2.0_pReal) + & + 0.75_pReal*(s12**2.0_pReal+s23**2.0_pReal+s31**2.0_pReal))**(0.5_pReal)/3.0_pReal + return + + ENDFUNCTION + !******************************************************************** ! determinant of a 3x3 matrix !******************************************************************** diff --git a/code/mesh.f90 b/code/mesh.f90 index fd9dd9e8a..677dcc316 100644 --- a/code/mesh.f90 +++ b/code/mesh.f90 @@ -240,7 +240,7 @@ call mesh_build_ipAreas() call mesh_tell_statistics() - parallelExecution = (mesh_Nelems == mesh_NcpElems) ! plus potential STOP if non-local constitutive present + parallelExecution = (parallelExecution .and. (mesh_Nelems == mesh_NcpElems)) ! plus potential killer from non-local constitutive else call IO_error(101) ! cannot open input file endif diff --git a/code/mpie_cpfem_marc.f90 b/code/mpie_cpfem_marc.f90 index c478ec6f2..5a5b783af 100644 --- a/code/mpie_cpfem_marc.f90 +++ b/code/mpie_cpfem_marc.f90 @@ -212,7 +212,6 @@ subroutine hypela2(& outdatedByNewInc = .true. terminallyIll = .false. cycleCounter = 0 - !$OMP CRITICAL (write2out) write (6,'(i6,x,i2,x,a)') n(1),nn,'<< hypela2 >> lastIncConverged + outdated'; call flush(6) !$OMP END CRITICAL (write2out) diff --git a/code/numerics.f90 b/code/numerics.f90 index 0a6f08453..437a27a83 100644 --- a/code/numerics.f90 +++ b/code/numerics.f90 @@ -35,8 +35,9 @@ real(pReal) relevantStrain, & ! strain relMax_RGC, & ! relative maximum of RGC residuum pPert_RGC, & ! perturbation for computing RGC penalty tangent xSmoo_RGC, & ! RGC penalty smoothing parameter (hyperbolic tangent) - ratePower_RGC, & ! power (sensitivity rate) of numerical viscosity in RGC scheme + viscPower_RGC, & ! power (sensitivity rate) of numerical viscosity in RGC scheme viscModus_RGC, & ! stress modulus of RGC numerical viscosity + refRelaxRate_RGC, & ! reference relaxation rate in RGC viscosity maxdRelax_RGC, & ! threshold of maximum relaxation vector increment (if exceed this then cutback) maxVolDiscr_RGC, & ! threshold of maximum volume discrepancy allowed volDiscrMod_RGC, & ! stiffness of RGC volume discrepancy (zero = without volume discrepancy constraint) @@ -111,8 +112,9 @@ subroutine numerics_init() relMax_RGC = 1.0e+2 pPert_RGC = 1.0e-7 xSmoo_RGC = 1.0e-5 - ratePower_RGC = 1.0e+0 ! Newton viscosity (linear model) + viscPower_RGC = 1.0e+0 ! Newton viscosity (linear model) viscModus_RGC = 0.0e+0 ! No viscosity is applied + refRelaxRate_RGC = 1.0e-3 maxdRelax_RGC = 1.0e+0 maxVolDiscr_RGC = 1.0e-5 ! tolerance for volume discrepancy allowed volDiscrMod_RGC = 1.0e+12 @@ -189,10 +191,12 @@ subroutine numerics_init() pPert_RGC = IO_floatValue(line,positions,2) case ('relevantmismatch_rgc') xSmoo_RGC = IO_floatValue(line,positions,2) - case ('viscosityrate_rgc') - ratePower_RGC = IO_floatValue(line,positions,2) + case ('viscositypower_rgc') + viscPower_RGC = IO_floatValue(line,positions,2) case ('viscositymodulus_rgc') viscModus_RGC = IO_floatValue(line,positions,2) + case ('refrelaxationrate_rgc') + refRelaxRate_RGC = IO_floatValue(line,positions,2) case ('maxrelaxation_rgc') maxdRelax_RGC = IO_floatValue(line,positions,2) case ('maxvoldiscrepancy_rgc') @@ -292,8 +296,9 @@ subroutine numerics_init() if (relMax_RGC <= 0.0_pReal) call IO_error(275) if (pPert_RGC <= 0.0_pReal) call IO_error(276) !! oops !! if (xSmoo_RGC <= 0.0_pReal) call IO_error(277) - if (ratePower_RGC < 0.0_pReal) call IO_error(278) + if (viscPower_RGC < 0.0_pReal) call IO_error(278) if (viscModus_RGC < 0.0_pReal) call IO_error(278) + if (refRelaxRate_RGC <= 0.0_pReal) call IO_error(278) if (maxdRelax_RGC <= 0.0_pReal) call IO_error(288) if (maxVolDiscr_RGC <= 0.0_pReal) call IO_error(289) if (volDiscrMod_RGC < 0.0_pReal) call IO_error(289)