From 82326ed8126cb6a7bc6e91a1fe77d6ef1e446044 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Thu, 11 Jun 2020 08:21:11 +0200 Subject: [PATCH 1/8] drop support for ping-pong scheme --- src/DAMASK_marc.f90 | 33 +-------------------------------- src/crystallite.f90 | 2 +- src/numerics.f90 | 5 ----- 3 files changed, 2 insertions(+), 38 deletions(-) diff --git a/src/DAMASK_marc.f90 b/src/DAMASK_marc.f90 index fa0711b02..dae77fb5d 100644 --- a/src/DAMASK_marc.f90 +++ b/src/DAMASK_marc.f90 @@ -305,36 +305,6 @@ subroutine hypela2(d,g,e,de,s,t,dt,ngens,m,nn,kcus,matus,ndi,nshear,disp, & endif ! convergence treatment end - if (usePingPong) then - calcMode(nn,cp_en) = .not. calcMode(nn,cp_en) ! ping pong (calc <--> collect) - if (calcMode(nn,cp_en)) then ! now --- CALC --- - computationMode = CPFEM_CALCRESULTS - if (lastLovl /= lovl) then ! first after ping pong - call debug_reset() ! resets debugging - outdatedFFN1 = .false. - cycleCounter = cycleCounter + 1 - !mesh_cellnode = mesh_build_cellnodes() ! update cell node coordinates - !call mesh_build_ipCoordinates() ! update ip coordinates - endif - if (outdatedByNewInc) then - computationMode = ior(computationMode,CPFEM_AGERESULTS) ! calc and age results - outdatedByNewInc = .false. ! reset flag - endif - else ! now --- COLLECT --- - computationMode = CPFEM_COLLECT ! plain collect - if (lastLovl /= lovl .and. & .not. terminallyIll) & - call debug_info() ! first after ping pong reports (meaningful) debugging - if (lastIncConverged) then - computationMode = ior(computationMode,CPFEM_BACKUPJACOBIAN) ! collect and backup Jacobian after convergence - lastIncConverged = .false. ! reset flag - endif - !do node = 1,theMesh%elem%nNodes - !CPnodeID = mesh_element(4+node,cp_en) - !mesh_node(1:ndeg,CPnodeID) = mesh_node0(1:ndeg,CPnodeID) + numerics_unitlength * dispt(1:ndeg,node) - !enddo - endif - - else ! --- PLAIN MODE --- computationMode = CPFEM_CALCRESULTS ! always calc if (lastLovl /= lovl) then if (.not. terminallyIll) & @@ -353,7 +323,6 @@ subroutine hypela2(d,g,e,de,s,t,dt,ngens,m,nn,kcus,matus,ndi,nshear,disp, & computationMode = ior(computationMode,CPFEM_BACKUPJACOBIAN) ! backup Jacobian after convergence lastIncConverged = .false. ! reset flag endif - endif theTime = cptim ! record current starting time theDelta = timinc ! record current time increment @@ -362,7 +331,7 @@ subroutine hypela2(d,g,e,de,s,t,dt,ngens,m,nn,kcus,matus,ndi,nshear,disp, & endif lastLovl = lovl ! record lovl - call CPFEM_general(computationMode,usePingPong,ffn,ffn1,t(1),timinc,m(1),nn,stress,ddsdde) + call CPFEM_general(computationMode,.false.,ffn,ffn1,t(1),timinc,m(1),nn,stress,ddsdde) d = ddsdde(1:ngens,1:ngens) s = stress(1:ndi+nshear) diff --git a/src/crystallite.f90 b/src/crystallite.f90 index 929fec862..c9fef3a51 100644 --- a/src/crystallite.f90 +++ b/src/crystallite.f90 @@ -245,7 +245,7 @@ subroutine crystallite_init enddo !$OMP END PARALLEL DO - if(any(plasticState%nonlocal) .and. .not. usePingPong) call IO_error(601) ! exit if nonlocal but no ping-pong ToDo: Why not check earlier? or in nonlocal? + !if(any(plasticState%nonlocal) .and. .not. usePingPong) call IO_error(601) crystallite_partionedFp0 = crystallite_Fp0 crystallite_partionedFi0 = crystallite_Fi0 diff --git a/src/numerics.f90 b/src/numerics.f90 index b0163aee3..77c3397de 100644 --- a/src/numerics.f90 +++ b/src/numerics.f90 @@ -28,8 +28,6 @@ module numerics numerics_unitlength = 1.0_pReal, & !< determines the physical length of one computational length unit charLength = 1.0_pReal, & !< characteristic length scale for gradient problems residualStiffness = 1.0e-6_pReal !< non-zero residual damage - logical, protected, public :: & - usePingPong = .true. !-------------------------------------------------------------------------------------------------- ! field parameters: @@ -133,8 +131,6 @@ subroutine numerics_init defgradTolerance = IO_floatValue(line,chunkPos,2) case ('ijacostiffness') iJacoStiffness = IO_intValue(line,chunkPos,2) - case ('usepingpong') - usepingpong = IO_intValue(line,chunkPos,2) > 0 case ('unitlength') numerics_unitlength = IO_floatValue(line,chunkPos,2) @@ -221,7 +217,6 @@ subroutine numerics_init ! writing parameters to output write(6,'(a24,1x,es8.1)') ' defgradTolerance: ',defgradTolerance write(6,'(a24,1x,i8)') ' iJacoStiffness: ',iJacoStiffness - write(6,'(a24,1x,L8)') ' use ping pong scheme: ',usepingpong write(6,'(a24,1x,es8.1,/)')' unitlength: ',numerics_unitlength !-------------------------------------------------------------------------------------------------- From b353129ba88000e3c21542097e4380af2bfc947a Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Thu, 11 Jun 2020 08:36:21 +0200 Subject: [PATCH 2/8] cleaning --- src/CPFEM.f90 | 74 +++------------------------------------------ src/DAMASK_marc.f90 | 4 +-- 2 files changed, 6 insertions(+), 72 deletions(-) diff --git a/src/CPFEM.f90 b/src/CPFEM.f90 index 7791d2b12..7ab6a2ece 100644 --- a/src/CPFEM.f90 +++ b/src/CPFEM.f90 @@ -124,7 +124,7 @@ end subroutine CPFEM_init !-------------------------------------------------------------------------------------------------- !> @brief perform initialization at first call, update variables and call the actual material model !-------------------------------------------------------------------------------------------------- -subroutine CPFEM_general(mode, parallelExecution, ffn, ffn1, temperature_inp, dt, elFE, ip, cauchyStress, jacobian) +subroutine CPFEM_general(mode, ffn, ffn1, temperature_inp, dt, elFE, ip, cauchyStress, jacobian) integer(pInt), intent(in) :: elFE, & !< FE element number ip !< integration point number @@ -133,17 +133,14 @@ subroutine CPFEM_general(mode, parallelExecution, ffn, ffn1, temperature_inp, dt ffn1 !< deformation gradient for t=t1 integer(pInt), intent(in) :: mode !< computation mode 1: regular computation plus aging of results real(pReal), intent(in) :: temperature_inp !< temperature - logical, intent(in) :: parallelExecution !< flag indicating parallel computation of requested IPs real(pReal), dimension(6), intent(out) :: cauchyStress !< stress as 6 vector real(pReal), dimension(6,6), intent(out) :: jacobian !< jacobian as 66 tensor (Consistent tangent dcs/dE) real(pReal) J_inverse, & ! inverse of Jacobian rnd - real(pReal), dimension (3,3) :: Kirchhoff, & ! Piola-Kirchhoff stress - cauchyStress33 ! stress vector + real(pReal), dimension (3,3) :: Kirchhoff ! Piola-Kirchhoff stress real(pReal), dimension (3,3,3,3) :: H_sym, & - H, & - jacobian3333 ! jacobian in Matrix notation + H integer(pInt) elCP, & ! crystal plasticity element number i, j, k, l, m, n, ph, homog, mySource @@ -178,9 +175,6 @@ subroutine CPFEM_general(mode, parallelExecution, ffn, ffn1, temperature_inp, dt if (iand(mode, CPFEM_AGERESULTS) /= 0_pInt) call CPFEM_forward - !*** collection of FEM input with returning of randomize odd stress and jacobian - !* If no parallel execution is required, there is no need to collect FEM input - if (.not. parallelExecution) then chosenThermal1: select case (thermal_type(material_homogenizationAt(elCP))) case (THERMAL_conduction_ID) chosenThermal1 temperature(material_homogenizationAt(elCP))%p(thermalMapping(material_homogenizationAt(elCP))%p(ip,elCP)) = & @@ -189,38 +183,11 @@ subroutine CPFEM_general(mode, parallelExecution, ffn, ffn1, temperature_inp, dt materialpoint_F0(1:3,1:3,ip,elCP) = ffn materialpoint_F(1:3,1:3,ip,elCP) = ffn1 - elseif (iand(mode, CPFEM_COLLECT) /= 0_pInt) then - call random_number(rnd) - if (rnd < 0.5_pReal) rnd = rnd - 1.0_pReal - CPFEM_cs(1:6,ip,elCP) = rnd * ODD_STRESS - CPFEM_dcsde(1:6,1:6,ip,elCP) = ODD_JACOBIAN * math_identity2nd(6) - chosenThermal2: select case (thermal_type(material_homogenizationAt(elCP))) - case (THERMAL_conduction_ID) chosenThermal2 - temperature(material_homogenizationAt(elCP))%p(thermalMapping(material_homogenizationAt(elCP))%p(ip,elCP)) = & - temperature_inp - end select chosenThermal2 - materialpoint_F0(1:3,1:3,ip,elCP) = ffn - materialpoint_F(1:3,1:3,ip,elCP) = ffn1 - CPFEM_calc_done = .false. - endif - - !*** calculation of stress and jacobian if (iand(mode, CPFEM_CALCRESULTS) /= 0_pInt) then !*** deformation gradient outdated or any actual deformation gradient differs more than relevantStrain from the stored one - validCalculation: if (terminallyIll & - .or. outdatedFFN1 & - .or. any(abs(ffn1 - materialpoint_F(1:3,1:3,ip,elCP)) > defgradTolerance)) then - if (any(abs(ffn1 - materialpoint_F(1:3,1:3,ip,elCP)) > defgradTolerance)) then - if (iand(debug_level(debug_CPFEM), debug_levelBasic) /= 0_pInt) then - write(6,'(a,1x,i8,1x,i2)') '<< CPFEM >> OUTDATED at elFE ip',elFE,ip - write(6,'(a,/,3(12x,3(f10.6,1x),/))') '<< CPFEM >> FFN1 old:',& - transpose(materialpoint_F(1:3,1:3,ip,elCP)) - write(6,'(a,/,3(12x,3(f10.6,1x),/))') '<< CPFEM >> FFN1 now:',transpose(ffn1) - endif - outdatedFFN1 = .true. - endif + validCalculation: if (terminallyIll .or. outdatedFFN1) then call random_number(rnd) if (rnd < 0.5_pReal) rnd = rnd - 1.0_pReal CPFEM_cs(1:6,ip,elCP) = ODD_STRESS * rnd @@ -229,23 +196,12 @@ subroutine CPFEM_general(mode, parallelExecution, ffn, ffn1, temperature_inp, dt !*** deformation gradient is not outdated else validCalculation updateJaco = mod(cycleCounter,iJacoStiffness) == 0 - !* no parallel computation, so we use just one single elFE and ip for computation - if (.not. parallelExecution) then FEsolving_execElem = elCP FEsolving_execIP = ip if (iand(debug_level(debug_CPFEM), debug_levelExtensive) /= 0_pInt) & write(6,'(a,i8,1x,i2)') '<< CPFEM >> calculation for elFE ip ',elFE,ip call materialpoint_stressAndItsTangent(updateJaco, dt) - !* parallel computation and calulation not yet done - elseif (.not. CPFEM_calc_done) then - if (iand(debug_level(debug_CPFEM), debug_levelExtensive) /= 0_pInt) & - write(6,'(a,i8,a,i8)') '<< CPFEM >> calculation for elements ',FEsolving_execElem(1),& - ' to ',FEsolving_execElem(2) - call materialpoint_stressAndItsTangent(updateJaco, dt) - CPFEM_calc_done = .true. - endif - !* map stress and stiffness (or return odd values if terminally ill) terminalIllness: if (terminallyIll) then @@ -300,28 +256,6 @@ subroutine CPFEM_general(mode, parallelExecution, ffn, ffn1, temperature_inp, dt cauchyStress = CPFEM_cs (1:6, ip,elCP) jacobian = CPFEM_dcsdE(1:6,1:6,ip,elCP) - - !*** remember extreme values of stress ... - cauchyStress33 = math_6toSym33(CPFEM_cs(1:6,ip,elCP),weighted=.false.) - if (maxval(cauchyStress33) > debug_stressMax) then - debug_stressMaxLocation = [elCP, ip] - debug_stressMax = maxval(cauchyStress33) - endif - if (minval(cauchyStress33) < debug_stressMin) then - debug_stressMinLocation = [elCP, ip] - debug_stressMin = minval(cauchyStress33) - endif - !*** ... and Jacobian - jacobian3333 = math_66toSym3333(CPFEM_dcsdE(1:6,1:6,ip,elCP),weighted=.false.) - if (maxval(jacobian3333) > debug_jacobianMax) then - debug_jacobianMaxLocation = [elCP, ip] - debug_jacobianMax = maxval(jacobian3333) - endif - if (minval(jacobian3333) < debug_jacobianMin) then - debug_jacobianMinLocation = [elCP, ip] - debug_jacobianMin = minval(jacobian3333) - endif - end subroutine CPFEM_general diff --git a/src/DAMASK_marc.f90 b/src/DAMASK_marc.f90 index dae77fb5d..fbff127e9 100644 --- a/src/DAMASK_marc.f90 +++ b/src/DAMASK_marc.f90 @@ -261,7 +261,7 @@ subroutine hypela2(d,g,e,de,s,t,dt,ngens,m,nn,kcus,matus,ndi,nshear,disp, & endif !$ defaultNumThreadsInt = omp_get_num_threads() ! remember number of threads set by Marc - !$ call omp_set_num_threads(DAMASK_NumThreadsInt) ! set number of threads for parallel execution set by DAMASK_NUM_THREADS + !$ call omp_set_num_threads(1) ! set number of threads for parallel execution set by DAMASK_NUM_THREADS if (.not. CPFEM_init_done) call CPFEM_initAll(m(1),nn) @@ -331,7 +331,7 @@ subroutine hypela2(d,g,e,de,s,t,dt,ngens,m,nn,kcus,matus,ndi,nshear,disp, & endif lastLovl = lovl ! record lovl - call CPFEM_general(computationMode,.false.,ffn,ffn1,t(1),timinc,m(1),nn,stress,ddsdde) + call CPFEM_general(computationMode,ffn,ffn1,t(1),timinc,m(1),nn,stress,ddsdde) d = ddsdde(1:ngens,1:ngens) s = stress(1:ndi+nshear) From 579ced6a52aab319c3b4630459648438be12500b Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Thu, 11 Jun 2020 08:44:24 +0200 Subject: [PATCH 3/8] removed global public variables --- src/DAMASK_marc.f90 | 3 --- src/crystallite.f90 | 3 --- src/debug.f90 | 55 +-------------------------------------------- 3 files changed, 1 insertion(+), 60 deletions(-) diff --git a/src/DAMASK_marc.f90 b/src/DAMASK_marc.f90 index fbff127e9..a653bcc00 100644 --- a/src/DAMASK_marc.f90 +++ b/src/DAMASK_marc.f90 @@ -307,9 +307,6 @@ subroutine hypela2(d,g,e,de,s,t,dt,ngens,m,nn,kcus,matus,ndi,nshear,disp, & computationMode = CPFEM_CALCRESULTS ! always calc if (lastLovl /= lovl) then - if (.not. terminallyIll) & - call debug_info() ! first reports (meaningful) debugging - call debug_reset() ! and resets debugging outdatedFFN1 = .false. cycleCounter = cycleCounter + 1 !mesh_cellnode = mesh_build_cellnodes() ! update cell node coordinates diff --git a/src/crystallite.f90 b/src/crystallite.f90 index c9fef3a51..14b624b0a 100644 --- a/src/crystallite.f90 +++ b/src/crystallite.f90 @@ -276,9 +276,6 @@ subroutine crystallite_init write(6,'(a42,1x,i10)') 'max # of constituents/integration point: ', cMax flush(6) endif - - call debug_info - call debug_reset #endif end subroutine crystallite_init diff --git a/src/debug.f90 b/src/debug.f90 index a0947c410..7564c3037 100644 --- a/src/debug.f90 +++ b/src/debug.f90 @@ -49,26 +49,11 @@ module debug debug_i = 1, & debug_g = 1 - integer, dimension(2), public :: & - debug_stressMaxLocation = 0, & - debug_stressMinLocation = 0, & - debug_jacobianMaxLocation = 0, & - debug_jacobianMinLocation = 0 - - - real(pReal), public :: & - debug_stressMax = -huge(1.0_pReal), & - debug_stressMin = huge(1.0_pReal), & - debug_jacobianMax = -huge(1.0_pReal), & - debug_jacobianMin = huge(1.0_pReal) - #ifdef PETSc character(len=1024), parameter, public :: & PETSCDEBUG = ' -snes_view -snes_monitor ' #endif - public :: debug_init, & - debug_reset, & - debug_info + public :: debug_init contains @@ -230,42 +215,4 @@ subroutine debug_init end subroutine debug_init - -!-------------------------------------------------------------------------------------------------- -!> @brief resets all debug values -!-------------------------------------------------------------------------------------------------- -subroutine debug_reset - - debug_stressMaxLocation = 0 - debug_stressMinLocation = 0 - debug_jacobianMaxLocation = 0 - debug_jacobianMinLocation = 0 - debug_stressMax = -huge(1.0_pReal) - debug_stressMin = huge(1.0_pReal) - debug_jacobianMax = -huge(1.0_pReal) - debug_jacobianMin = huge(1.0_pReal) - -end subroutine debug_reset - - -!-------------------------------------------------------------------------------------------------- -!> @brief writes debug statements to standard out -!-------------------------------------------------------------------------------------------------- -subroutine debug_info - - !$OMP CRITICAL (write2out) - debugOutputCPFEM: if (iand(debug_level(debug_CPFEM),debug_LEVELBASIC) /= 0 & - .and. any(debug_stressMinLocation /= 0) & - .and. any(debug_stressMaxLocation /= 0) ) then - write(6,'(2/,a,/)') ' Extreme values of returned stress and Jacobian' - write(6,'(a39)') ' value el ip' - write(6,'(a14,1x,e12.3,1x,i8,1x,i4)') ' stress min :', debug_stressMin, debug_stressMinLocation - write(6,'(a14,1x,e12.3,1x,i8,1x,i4)') ' max :', debug_stressMax, debug_stressMaxLocation - write(6,'(a14,1x,e12.3,1x,i8,1x,i4)') ' Jacobian min :', debug_jacobianMin, debug_jacobianMinLocation - write(6,'(a14,1x,e12.3,1x,i8,1x,i4,/)') ' max :', debug_jacobianMax, debug_jacobianMaxLocation - endif debugOutputCPFEM - !$OMP END CRITICAL (write2out) - -end subroutine debug_info - end module debug From e5c9380bac2be818b6886f01a006def560b7a7bc Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Thu, 11 Jun 2020 08:51:17 +0200 Subject: [PATCH 4/8] cleaning --- src/CPFEM.f90 | 12 ++++-------- src/DAMASK_marc.f90 | 10 +--------- src/marc/discretization_marc.f90 | 10 ++-------- 3 files changed, 7 insertions(+), 25 deletions(-) diff --git a/src/CPFEM.f90 b/src/CPFEM.f90 index 7ab6a2ece..5acfa9aa1 100644 --- a/src/CPFEM.f90 +++ b/src/CPFEM.f90 @@ -43,7 +43,6 @@ module CPFEM theTime = 0.0_pReal, & !< needs description theDelta = 0.0_pReal logical, public :: & - outdatedFFN1 = .false., & !< needs description lastIncConverged = .false., & !< needs description outdatedByNewInc = .false. !< needs description @@ -68,12 +67,9 @@ contains !-------------------------------------------------------------------------------------------------- -!> @brief call (thread safe) all module initializations +!> @brief call all module initializations !-------------------------------------------------------------------------------------------------- -subroutine CPFEM_initAll(el,ip) - - integer(pInt), intent(in) :: el, & !< FE el number - ip !< FE integration point number +subroutine CPFEM_initAll CPFEM_init_done = .true. call DAMASK_interface_init @@ -88,7 +84,7 @@ subroutine CPFEM_initAll(el,ip) call YAML_init call HDF5_utilities_init call results_init(.false.) - call discretization_marc_init(ip, el) + call discretization_marc_init call lattice_init call material_init(.false.) call constitutive_init @@ -187,7 +183,7 @@ subroutine CPFEM_general(mode, ffn, ffn1, temperature_inp, dt, elFE, ip, cauchyS if (iand(mode, CPFEM_CALCRESULTS) /= 0_pInt) then !*** deformation gradient outdated or any actual deformation gradient differs more than relevantStrain from the stored one - validCalculation: if (terminallyIll .or. outdatedFFN1) then + validCalculation: if (terminallyIll) then call random_number(rnd) if (rnd < 0.5_pReal) rnd = rnd - 1.0_pReal CPFEM_cs(1:6,ip,elCP) = ODD_STRESS * rnd diff --git a/src/DAMASK_marc.f90 b/src/DAMASK_marc.f90 index a653bcc00..6de01273d 100644 --- a/src/DAMASK_marc.f90 +++ b/src/DAMASK_marc.f90 @@ -43,9 +43,6 @@ module DAMASK_interface logical, protected, public :: symmetricSolver character(len=*), parameter, public :: INPUTFILEEXTENSION = '.dat' - logical, dimension(:,:), public, allocatable :: & - calcMode !< calculate or collect (ping pong scheme) - public :: & DAMASK_interface_init, & getSolverJobName @@ -263,7 +260,7 @@ subroutine hypela2(d,g,e,de,s,t,dt,ngens,m,nn,kcus,matus,ndi,nshear,disp, & !$ defaultNumThreadsInt = omp_get_num_threads() ! remember number of threads set by Marc !$ call omp_set_num_threads(1) ! set number of threads for parallel execution set by DAMASK_NUM_THREADS - if (.not. CPFEM_init_done) call CPFEM_initAll(m(1),nn) + if (.not. CPFEM_init_done) call CPFEM_initAll computationMode = 0 ! save initialization value, since it does not result in any calculation if (lovl == 4 ) then ! jacobian requested by marc @@ -277,20 +274,17 @@ subroutine hypela2(d,g,e,de,s,t,dt,ngens,m,nn,kcus,matus,ndi,nshear,disp, & if (inc == 0) then ! >> start of analysis << lastIncConverged = .false. ! no Jacobian backup outdatedByNewInc = .false. ! no aging of state - calcMode = .false. ! pretend last step was collection lastLovl = lovl ! pretend that this is NOT the first after a lovl change write(6,'(a,i6,1x,i2)') '<< HYPELA2 >> start of analysis..! ',m(1),nn flush(6) else if (inc - theInc > 1) then ! >> restart of broken analysis << lastIncConverged = .false. ! no Jacobian backup outdatedByNewInc = .false. ! no aging of state - calcMode = .true. ! pretend last step was calculation write(6,'(a,i6,1x,i2)') '<< HYPELA2 >> restart of analysis..! ',m(1),nn flush(6) else ! >> just the next inc << lastIncConverged = .true. ! request Jacobian backup outdatedByNewInc = .true. ! request aging of state - calcMode = .true. ! assure last step was calculation write(6,'(a,i6,1x,i2)') '<< HYPELA2 >> new increment..! ',m(1),nn flush(6) endif @@ -299,7 +293,6 @@ subroutine hypela2(d,g,e,de,s,t,dt,ngens,m,nn,kcus,matus,ndi,nshear,disp, & outdatedByNewInc = .false. ! no aging of state terminallyIll = .false. cycleCounter = -1 ! first calc step increments this to cycle = 0 - calcMode = .true. ! pretend last step was calculation write(6,'(a,i6,1x,i2)') '<< HYPELA2 >> cutback detected..! ',m(1),nn flush(6) endif ! convergence treatment end @@ -307,7 +300,6 @@ subroutine hypela2(d,g,e,de,s,t,dt,ngens,m,nn,kcus,matus,ndi,nshear,disp, & computationMode = CPFEM_CALCRESULTS ! always calc if (lastLovl /= lovl) then - outdatedFFN1 = .false. cycleCounter = cycleCounter + 1 !mesh_cellnode = mesh_build_cellnodes() ! update cell node coordinates !call mesh_build_ipCoordinates() ! update ip coordinates diff --git a/src/marc/discretization_marc.f90 b/src/marc/discretization_marc.f90 index 0221ad693..d1026d568 100644 --- a/src/marc/discretization_marc.f90 +++ b/src/marc/discretization_marc.f90 @@ -45,10 +45,8 @@ contains !> @brief initializes the mesh by calling all necessary private routines the mesh module !! Order and routines strongly depend on type of solver !-------------------------------------------------------------------------------------------------- -subroutine discretization_marc_init(ip,el) +subroutine discretization_marc_init - integer, intent(in) :: el, ip - real(pReal), dimension(:,:), allocatable :: & node0_elem, & !< node x,y,z coordinates (initially!) node0_cell @@ -70,7 +68,7 @@ subroutine discretization_marc_init(ip,el) real(pReal), dimension(:,:,:,:),allocatable :: & unscaledNormals - write(6,'(/,a)') ' <<<+- mesh init -+>>>'; flush(6) + write(6,'(/,a)') ' <<<+- discretization_marc init -+>>>'; flush(6) mesh_unitlength = numerics_unitlength ! set physical extent of a length unit in mesh @@ -83,10 +81,6 @@ subroutine discretization_marc_init(ip,el) FEsolving_execElem = [1,nElems] FEsolving_execIP = [1,elem%nIPs] - allocate(calcMode(elem%nIPs,nElems),source=.false.) ! pretend to have collected what first call is asking (F = I) - calcMode(ip,mesh_FEM2DAMASK_elem(el)) = .true. ! first ip,el needs to be already pingponged to "calc" - - allocate(cellNodeDefinition(elem%nNodes-1)) allocate(connectivity_cell(elem%NcellNodesPerCell,elem%nIPs,nElems)) call buildCells(connectivity_cell,cellNodeDefinition,& From e952ab71278602a9215b6afe7bf2f3b405aa0ed4 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Mon, 15 Jun 2020 23:12:49 +0200 Subject: [PATCH 5/8] bugfix do not access unitinialized memory --- src/homogenization.f90 | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/homogenization.f90 b/src/homogenization.f90 index 62150899e..4fdded94e 100644 --- a/src/homogenization.f90 +++ b/src/homogenization.f90 @@ -249,7 +249,8 @@ subroutine materialpoint_stressAndItsTangent(updateJaco,dt) NiterationHomog = 0 cutBackLooping: do while (.not. terminallyIll .and. & - any(subStep(:,FEsolving_execELem(1):FEsolving_execElem(2)) > num%subStepMinHomog)) + any(subStep(FEsolving_execIP(1):FEsolving_execIP(2),& + FEsolving_execElem(1):FEsolving_execElem(2)) > num%subStepMinHomog)) !$OMP PARALLEL DO PRIVATE(myNgrains) elementLooping1: do e = FEsolving_execElem(1),FEsolving_execElem(2) From 0a9902818c8a3731a81647a152493c57850a7998 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Tue, 16 Jun 2020 06:34:12 +0200 Subject: [PATCH 6/8] polishing --- src/DAMASK_marc.f90 | 272 ++++++++++++++++++++++---------------------- 1 file changed, 134 insertions(+), 138 deletions(-) diff --git a/src/DAMASK_marc.f90 b/src/DAMASK_marc.f90 index 6de01273d..5349608c9 100644 --- a/src/DAMASK_marc.f90 +++ b/src/DAMASK_marc.f90 @@ -172,62 +172,62 @@ subroutine hypela2(d,g,e,de,s,t,dt,ngens,m,nn,kcus,matus,ndi,nshear,disp, & dispt,coord,ffn,frotn,strechn,eigvn,ffn1,frotn1, & strechn1,eigvn1,ncrd,itel,ndeg,ndm,nnode, & jtype,lclass,ifr,ifu) - use prec - use DAMASK_interface - use numerics - use FEsolving - use debug - use discretization_marc - use homogenization - use CPFEM + use prec + use DAMASK_interface + use numerics + use FEsolving + use debug + use discretization_marc + use homogenization + use CPFEM - implicit none -!$ include "omp_lib.h" ! the openMP function library - integer, intent(in) :: & ! according to MSC.Marc 2012 Manual D - ngens, & !< size of stress-strain law - nn, & !< integration point number - ndi, & !< number of direct components - nshear, & !< number of shear components - ncrd, & !< number of coordinates - itel, & !< dimension of F and R, either 2 or 3 - ndeg, & !< number of degrees of freedom - ndm, & !< not specified in MSC.Marc 2012 Manual D - nnode, & !< number of nodes per element - jtype, & !< element type - ifr, & !< set to 1 if R has been calculated - ifu !< set to 1 if stretch has been calculated - integer, dimension(2), intent(in) :: & ! according to MSC.Marc 2012 Manual D - m, & !< (1) user element number, (2) internal element number - matus, & !< (1) user material identification number, (2) internal material identification number - kcus, & !< (1) layer number, (2) internal layer number - lclass !< (1) element class, (2) 0: displacement, 1: low order Herrmann, 2: high order Herrmann - real(pReal), dimension(*), intent(in) :: & ! has dimension(1) according to MSC.Marc 2012 Manual D, but according to example hypela2.f dimension(*) - e, & !< total elastic strain - de, & !< increment of strain - dt !< increment of state variables - real(pReal), dimension(itel), intent(in) :: & ! according to MSC.Marc 2012 Manual D - strechn, & !< square of principal stretch ratios, lambda(i) at t=n - strechn1 !< square of principal stretch ratios, lambda(i) at t=n+1 - real(pReal), dimension(3,3), intent(in) :: & ! has dimension(itel,*) according to MSC.Marc 2012 Manual D, but we alway assume dimension(3,3) - ffn, & !< deformation gradient at t=n - ffn1 !< deformation gradient at t=n+1 - real(pReal), dimension(itel,*), intent(in) :: & ! according to MSC.Marc 2012 Manual D - frotn, & !< rotation tensor at t=n - eigvn, & !< i principal direction components for j eigenvalues at t=n - frotn1, & !< rotation tensor at t=n+1 - eigvn1 !< i principal direction components for j eigenvalues at t=n+1 - real(pReal), dimension(ndeg,*), intent(in) :: & ! according to MSC.Marc 2012 Manual D - disp, & !< incremental displacements - dispt !< displacements at t=n (at assembly, lovl=4) and displacements at t=n+1 (at stress recovery, lovl=6) - real(pReal), dimension(ncrd,*), intent(in) :: & ! according to MSC.Marc 2012 Manual D - coord !< coordinates - real(pReal), dimension(*), intent(inout) :: & ! according to MSC.Marc 2012 Manual D - t !< state variables (comes in at t=n, must be updated to have state variables at t=n+1) - real(pReal), dimension(ndi+nshear), intent(out) :: & ! has dimension(*) according to MSC.Marc 2012 Manual D, but we need to loop over it - s, & !< stress - should be updated by user - g !< change in stress due to temperature effects - real(pReal), dimension(ngens,ngens), intent(out) :: & ! according to MSC.Marc 2012 Manual D, but according to example hypela2.f dimension(ngens,*) - d !< stress-strain law to be formed + implicit none + include "omp_lib.h" ! the openMP function library + integer, intent(in) :: & ! according to MSC.Marc 2012 Manual D + ngens, & !< size of stress-strain law + nn, & !< integration point number + ndi, & !< number of direct components + nshear, & !< number of shear components + ncrd, & !< number of coordinates + itel, & !< dimension of F and R, either 2 or 3 + ndeg, & !< number of degrees of freedom + ndm, & !< not specified in MSC.Marc 2012 Manual D + nnode, & !< number of nodes per element + jtype, & !< element type + ifr, & !< set to 1 if R has been calculated + ifu !< set to 1 if stretch has been calculated + integer, dimension(2), intent(in) :: & ! according to MSC.Marc 2012 Manual D + m, & !< (1) user element number, (2) internal element number + matus, & !< (1) user material identification number, (2) internal material identification number + kcus, & !< (1) layer number, (2) internal layer number + lclass !< (1) element class, (2) 0: displacement, 1: low order Herrmann, 2: high order Herrmann + real(pReal), dimension(*), intent(in) :: & ! has dimension(1) according to MSC.Marc 2012 Manual D, but according to example hypela2.f dimension(*) + e, & !< total elastic strain + de, & !< increment of strain + dt !< increment of state variables + real(pReal), dimension(itel), intent(in) :: & ! according to MSC.Marc 2012 Manual D + strechn, & !< square of principal stretch ratios, lambda(i) at t=n + strechn1 !< square of principal stretch ratios, lambda(i) at t=n+1 + real(pReal), dimension(3,3), intent(in) :: & ! has dimension(itel,*) according to MSC.Marc 2012 Manual D, but we alway assume dimension(3,3) + ffn, & !< deformation gradient at t=n + ffn1 !< deformation gradient at t=n+1 + real(pReal), dimension(itel,*), intent(in) :: & ! according to MSC.Marc 2012 Manual D + frotn, & !< rotation tensor at t=n + eigvn, & !< i principal direction components for j eigenvalues at t=n + frotn1, & !< rotation tensor at t=n+1 + eigvn1 !< i principal direction components for j eigenvalues at t=n+1 + real(pReal), dimension(ndeg,*), intent(in) :: & ! according to MSC.Marc 2012 Manual D + disp, & !< incremental displacements + dispt !< displacements at t=n (at assembly, lovl=4) and displacements at t=n+1 (at stress recovery, lovl=6) + real(pReal), dimension(ncrd,*), intent(in) :: & ! according to MSC.Marc 2012 Manual D + coord !< coordinates + real(pReal), dimension(*), intent(inout) :: & ! according to MSC.Marc 2012 Manual D + t !< state variables (comes in at t=n, must be updated to have state variables at t=n+1) + real(pReal), dimension(ndi+nshear), intent(out) :: & ! has dimension(*) according to MSC.Marc 2012 Manual D, but we need to loop over it + s, & !< stress - should be updated by user + g !< change in stress due to temperature effects + real(pReal), dimension(ngens,ngens), intent(out) :: & ! according to MSC.Marc 2012 Manual D, but according to example hypela2.f dimension(ngens,*) + d !< stress-strain law to be formed !-------------------------------------------------------------------------------------------------- ! Marc common blocks are in fixed format so they have to be reformated to free format (f90) @@ -236,98 +236,94 @@ subroutine hypela2(d,g,e,de,s,t,dt,ngens,m,nn,kcus,matus,ndi,nshear,disp, & #include QUOTE(PASTE(./marc/include/concom,Marc4DAMASK)) ! concom is needed for inc, lovl #include QUOTE(PASTE(./marc/include/creeps,Marc4DAMASK)) ! creeps is needed for timinc (time increment) - logical :: cutBack - real(pReal), dimension(6) :: stress - real(pReal), dimension(6,6) :: ddsdde - integer :: computationMode, i, cp_en, node, CPnodeID - !$ integer(4) :: defaultNumThreadsInt !< default value set by Marc + logical :: cutBack + real(pReal), dimension(6) :: stress + real(pReal), dimension(6,6) :: ddsdde + integer :: computationMode, i, cp_en, node, CPnodeID + integer(4) :: defaultNumThreadsInt !< default value set by Marc + + if (iand(debug_level(debug_MARC),debug_LEVELBASIC) /= 0) then + write(6,'(a,/,i8,i8,i2)') ' MSC.MARC information on shape of element(2), IP:', m, nn + write(6,'(a,2(i1))') ' Jacobian: ', ngens,ngens + write(6,'(a,i1)') ' Direct stress: ', ndi + write(6,'(a,i1)') ' Shear stress: ', nshear + write(6,'(a,i2)') ' DoF: ', ndeg + write(6,'(a,i2)') ' Coordinates: ', ncrd + write(6,'(a,i12)') ' Nodes: ', nnode + write(6,'(a,i1)') ' Deformation gradient: ', itel + write(6,'(/,a,/,3(3(f12.7,1x)/))',advance='no') ' Deformation gradient at t=n:', & + transpose(ffn) + write(6,'(/,a,/,3(3(f12.7,1x)/))',advance='no') ' Deformation gradient at t=n+1:', & + transpose(ffn1) + endif - if (iand(debug_level(debug_MARC),debug_LEVELBASIC) /= 0) then - write(6,'(a,/,i8,i8,i2)') ' MSC.MARC information on shape of element(2), IP:', m, nn - write(6,'(a,2(i1))') ' Jacobian: ', ngens,ngens - write(6,'(a,i1)') ' Direct stress: ', ndi - write(6,'(a,i1)') ' Shear stress: ', nshear - write(6,'(a,i2)') ' DoF: ', ndeg - write(6,'(a,i2)') ' Coordinates: ', ncrd - write(6,'(a,i12)') ' Nodes: ', nnode - write(6,'(a,i1)') ' Deformation gradient: ', itel - write(6,'(/,a,/,3(3(f12.7,1x)/))',advance='no') ' Deformation gradient at t=n:', & - transpose(ffn) - write(6,'(/,a,/,3(3(f12.7,1x)/))',advance='no') ' Deformation gradient at t=n+1:', & - transpose(ffn1) - endif + defaultNumThreadsInt = omp_get_num_threads() ! remember number of threads set by Marc + call omp_set_num_threads(1) ! no openMP - !$ defaultNumThreadsInt = omp_get_num_threads() ! remember number of threads set by Marc - !$ call omp_set_num_threads(1) ! set number of threads for parallel execution set by DAMASK_NUM_THREADS + if (.not. CPFEM_init_done) call CPFEM_initAll - if (.not. CPFEM_init_done) call CPFEM_initAll + computationMode = 0 ! save initialization value, since it does not result in any calculation + if (lovl == 4 ) then ! jacobian requested by marc + if (timinc < theDelta .and. theInc == inc .and. lastLovl /= lovl) & ! first after cutback + computationMode = CPFEM_RESTOREJACOBIAN + elseif (lovl == 6) then ! stress requested by marc + computationMode = CPFEM_CALCRESULTS ! always calc + cp_en = mesh_FEM2DAMASK_elem(m(1)) + if (cptim > theTime .or. inc /= theInc) then ! reached "convergence" + terminallyIll = .false. + cycleCounter = -1 ! first calc step increments this to cycle = 0 + if (inc == 0) then ! >> start of analysis << + lastIncConverged = .false. + outdatedByNewInc = .false. + lastLovl = lovl ! pretend that this is NOT the first after a lovl change + write(6,'(a,i6,1x,i2)') '<< HYPELA2 >> start of analysis..! ',m(1),nn + else if (inc - theInc > 1) then ! >> restart of broken analysis << + lastIncConverged = .false. + outdatedByNewInc = .false. + write(6,'(a,i6,1x,i2)') '<< HYPELA2 >> restart of analysis..! ',m(1),nn + else ! >> just the next inc << + lastIncConverged = .true. + outdatedByNewInc = .true. + write(6,'(a,i6,1x,i2)') '<< HYPELA2 >> new increment..! ',m(1),nn + endif + else if ( timinc < theDelta ) then ! >> cutBack << + lastIncConverged = .false. + outdatedByNewInc = .false. + terminallyIll = .false. + cycleCounter = -1 ! first calc step increments this to cycle = 0 + write(6,'(a,i6,1x,i2)') '<< HYPELA2 >> cutback detected..! ',m(1),nn + endif ! convergence treatment end + flush(6) - computationMode = 0 ! save initialization value, since it does not result in any calculation - if (lovl == 4 ) then ! jacobian requested by marc - if (timinc < theDelta .and. theInc == inc .and. lastLovl /= lovl) & ! first after cutback - computationMode = CPFEM_RESTOREJACOBIAN - elseif (lovl == 6) then ! stress requested by marc - cp_en = mesh_FEM2DAMASK_elem(m(1)) - if (cptim > theTime .or. inc /= theInc) then ! reached "convergence" - terminallyIll = .false. - cycleCounter = -1 ! first calc step increments this to cycle = 0 - if (inc == 0) then ! >> start of analysis << - lastIncConverged = .false. ! no Jacobian backup - outdatedByNewInc = .false. ! no aging of state - lastLovl = lovl ! pretend that this is NOT the first after a lovl change - write(6,'(a,i6,1x,i2)') '<< HYPELA2 >> start of analysis..! ',m(1),nn - flush(6) - else if (inc - theInc > 1) then ! >> restart of broken analysis << - lastIncConverged = .false. ! no Jacobian backup - outdatedByNewInc = .false. ! no aging of state - write(6,'(a,i6,1x,i2)') '<< HYPELA2 >> restart of analysis..! ',m(1),nn - flush(6) - else ! >> just the next inc << - lastIncConverged = .true. ! request Jacobian backup - outdatedByNewInc = .true. ! request aging of state - write(6,'(a,i6,1x,i2)') '<< HYPELA2 >> new increment..! ',m(1),nn - flush(6) - endif - else if ( timinc < theDelta ) then ! >> cutBack << - lastIncConverged = .false. ! no Jacobian backup - outdatedByNewInc = .false. ! no aging of state - terminallyIll = .false. - cycleCounter = -1 ! first calc step increments this to cycle = 0 - write(6,'(a,i6,1x,i2)') '<< HYPELA2 >> cutback detected..! ',m(1),nn - flush(6) - endif ! convergence treatment end + if (lastLovl /= lovl) then + cycleCounter = cycleCounter + 1 + !mesh_cellnode = mesh_build_cellnodes() ! update cell node coordinates + !call mesh_build_ipCoordinates() ! update ip coordinates + endif + if (outdatedByNewInc) then + computationMode = ior(computationMode,CPFEM_AGERESULTS) + outdatedByNewInc = .false. + endif + if (lastIncConverged) then + computationMode = ior(computationMode,CPFEM_BACKUPJACOBIAN) + lastIncConverged = .false. + endif + theTime = cptim + theDelta = timinc + theInc = inc - computationMode = CPFEM_CALCRESULTS ! always calc - if (lastLovl /= lovl) then - cycleCounter = cycleCounter + 1 - !mesh_cellnode = mesh_build_cellnodes() ! update cell node coordinates - !call mesh_build_ipCoordinates() ! update ip coordinates - endif - if (outdatedByNewInc) then - computationMode = ior(computationMode,CPFEM_AGERESULTS) - outdatedByNewInc = .false. ! reset flag - endif - if (lastIncConverged) then - computationMode = ior(computationMode,CPFEM_BACKUPJACOBIAN) ! backup Jacobian after convergence - lastIncConverged = .false. ! reset flag - endif + endif + lastLovl = lovl - theTime = cptim ! record current starting time - theDelta = timinc ! record current time increment - theInc = inc ! record current increment number + call CPFEM_general(computationMode,ffn,ffn1,t(1),timinc,m(1),nn,stress,ddsdde) - endif - lastLovl = lovl ! record lovl + d = ddsdde(1:ngens,1:ngens) + s = stress(1:ndi+nshear) + g = 0.0_pReal + if(symmetricSolver) d = 0.5_pReal*(d+transpose(d)) - call CPFEM_general(computationMode,ffn,ffn1,t(1),timinc,m(1),nn,stress,ddsdde) - - d = ddsdde(1:ngens,1:ngens) - s = stress(1:ndi+nshear) - g = 0.0_pReal - if(symmetricSolver) d = 0.5_pReal*(d+transpose(d)) - - !$ call omp_set_num_threads(defaultNumThreadsInt) ! reset number of threads to stored default value + call omp_set_num_threads(defaultNumThreadsInt) ! reset number of threads to stored default value end subroutine hypela2 From 54aa5a67ff57ca7db298d06441226dcb3a08810f Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Tue, 16 Jun 2020 06:41:53 +0200 Subject: [PATCH 7/8] polishing --- src/CPFEM.f90 | 38 +++++++++++++------------------------- 1 file changed, 13 insertions(+), 25 deletions(-) diff --git a/src/CPFEM.f90 b/src/CPFEM.f90 index 5acfa9aa1..8c53eb9a4 100644 --- a/src/CPFEM.f90 +++ b/src/CPFEM.f90 @@ -48,15 +48,12 @@ module CPFEM logical, public, protected :: & CPFEM_init_done = .false. !< remember whether init has been done already - logical, private :: & - CPFEM_calc_done = .false. !< remember whether first ip has already calced the results integer(pInt), parameter, public :: & - CPFEM_COLLECT = 2_pInt**0_pInt, & - CPFEM_CALCRESULTS = 2_pInt**1_pInt, & - CPFEM_AGERESULTS = 2_pInt**2_pInt, & - CPFEM_BACKUPJACOBIAN = 2_pInt**3_pInt, & - CPFEM_RESTOREJACOBIAN = 2_pInt**4_pInt + CPFEM_CALCRESULTS = 2_pInt**0_pInt, & + CPFEM_AGERESULTS = 2_pInt**1_pInt, & + CPFEM_BACKUPJACOBIAN = 2_pInt**2_pInt, & + CPFEM_RESTOREJACOBIAN = 2_pInt**3_pInt public :: & CPFEM_general, & @@ -134,7 +131,7 @@ subroutine CPFEM_general(mode, ffn, ffn1, temperature_inp, dt, elFE, ip, cauchyS real(pReal) J_inverse, & ! inverse of Jacobian rnd - real(pReal), dimension (3,3) :: Kirchhoff ! Piola-Kirchhoff stress + real(pReal), dimension (3,3) :: Kirchhoff ! Piola-Kirchhoff stress real(pReal), dimension (3,3,3,3) :: H_sym, & H @@ -142,8 +139,8 @@ subroutine CPFEM_general(mode, ffn, ffn1, temperature_inp, dt, elFE, ip, cauchyS i, j, k, l, m, n, ph, homog, mySource logical updateJaco ! flag indicating if Jacobian has to be updated - real(pReal), parameter :: ODD_STRESS = 1e15_pReal, & !< return value for stress in case of ping pong dummy cycle - ODD_JACOBIAN = 1e50_pReal !< return value for jacobian in case of ping pong dummy cycle + real(pReal), parameter :: ODD_STRESS = 1e15_pReal, & !< return value for stress if terminallyIll + ODD_JACOBIAN = 1e50_pReal !< return value for jacobian if terminallyIll elCP = mesh_FEM2DAMASK_elem(elFE) @@ -167,10 +164,8 @@ subroutine CPFEM_general(mode, ffn, ffn1, temperature_inp, dt, elFE, ip, cauchyS if (iand(mode, CPFEM_RESTOREJACOBIAN) /= 0_pInt) & CPFEM_dcsde = CPFEM_dcsde_knownGood - !*** age results if (iand(mode, CPFEM_AGERESULTS) /= 0_pInt) call CPFEM_forward - chosenThermal1: select case (thermal_type(material_homogenizationAt(elCP))) case (THERMAL_conduction_ID) chosenThermal1 temperature(material_homogenizationAt(elCP))%p(thermalMapping(material_homogenizationAt(elCP))%p(ip,elCP)) = & @@ -179,26 +174,22 @@ subroutine CPFEM_general(mode, ffn, ffn1, temperature_inp, dt, elFE, ip, cauchyS materialpoint_F0(1:3,1:3,ip,elCP) = ffn materialpoint_F(1:3,1:3,ip,elCP) = ffn1 - !*** calculation of stress and jacobian if (iand(mode, CPFEM_CALCRESULTS) /= 0_pInt) then - !*** deformation gradient outdated or any actual deformation gradient differs more than relevantStrain from the stored one validCalculation: if (terminallyIll) then call random_number(rnd) if (rnd < 0.5_pReal) rnd = rnd - 1.0_pReal CPFEM_cs(1:6,ip,elCP) = ODD_STRESS * rnd CPFEM_dcsde(1:6,1:6,ip,elCP) = ODD_JACOBIAN * math_identity2nd(6) - !*** deformation gradient is not outdated else validCalculation updateJaco = mod(cycleCounter,iJacoStiffness) == 0 - FEsolving_execElem = elCP - FEsolving_execIP = ip - if (iand(debug_level(debug_CPFEM), debug_levelExtensive) /= 0_pInt) & - write(6,'(a,i8,1x,i2)') '<< CPFEM >> calculation for elFE ip ',elFE,ip - call materialpoint_stressAndItsTangent(updateJaco, dt) + FEsolving_execElem = elCP + FEsolving_execIP = ip + if (iand(debug_level(debug_CPFEM), debug_levelExtensive) /= 0_pInt) & + write(6,'(a,i8,1x,i2)') '<< CPFEM >> calculation for elFE ip ',elFE,ip + call materialpoint_stressAndItsTangent(updateJaco, dt) - !* map stress and stiffness (or return odd values if terminally ill) terminalIllness: if (terminallyIll) then call random_number(rnd) @@ -208,7 +199,7 @@ subroutine CPFEM_general(mode, ffn, ffn1, temperature_inp, dt, elFE, ip, cauchyS else terminalIllness - ! translate from P to CS + ! translate from P to sigma Kirchhoff = matmul(materialpoint_P(1:3,1:3,ip,elCP), transpose(materialpoint_F(1:3,1:3,ip,elCP))) J_inverse = 1.0_pReal / math_det33(materialpoint_F(1:3,1:3,ip,elCP)) CPFEM_cs(1:6,ip,elCP) = math_sym33to6(J_inverse * Kirchhoff,weighted=.false.) @@ -232,7 +223,6 @@ subroutine CPFEM_general(mode, ffn, ffn1, temperature_inp, dt, elFE, ip, cauchyS endif terminalIllness endif validCalculation - !* report stress and stiffness if ((iand(debug_level(debug_CPFEM), debug_levelExtensive) /= 0_pInt) & .and. ((debug_e == elCP .and. debug_i == ip) & .or. .not. iand(debug_level(debug_CPFEM), debug_levelSelective) /= 0_pInt)) then @@ -245,10 +235,8 @@ subroutine CPFEM_general(mode, ffn, ffn1, temperature_inp, dt, elFE, ip, cauchyS endif - !*** warn if stiffness close to zero if (all(abs(CPFEM_dcsdE(1:6,1:6,ip,elCP)) < 1e-10_pReal)) call IO_warning(601,elCP,ip) - !*** copy to output if using commercial FEM solver cauchyStress = CPFEM_cs (1:6, ip,elCP) jacobian = CPFEM_dcsdE(1:6,1:6,ip,elCP) From 06f6e15123a674efce1b5b1c9d2f3adfa1c145f4 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Tue, 16 Jun 2020 07:02:29 +0200 Subject: [PATCH 8/8] avoid public variables --- src/CPFEM.f90 | 18 ++---------------- src/DAMASK_marc.f90 | 19 +++++++++++++++++-- 2 files changed, 19 insertions(+), 18 deletions(-) diff --git a/src/CPFEM.f90 b/src/CPFEM.f90 index 8c53eb9a4..d285a1e01 100644 --- a/src/CPFEM.f90 +++ b/src/CPFEM.f90 @@ -35,19 +35,9 @@ module CPFEM CPFEM_dcsdE !< Cauchy stress tangent real(pReal), dimension (:,:,:,:), allocatable, private :: & CPFEM_dcsdE_knownGood !< known good tangent - integer(pInt), public :: & - cycleCounter = 0_pInt, & !< needs description - theInc = -1_pInt, & !< needs description - lastLovl = 0_pInt !< lovl in previous call to marc hypela2 - real(pReal), public :: & - theTime = 0.0_pReal, & !< needs description - theDelta = 0.0_pReal - logical, public :: & - lastIncConverged = .false., & !< needs description - outdatedByNewInc = .false. !< needs description - logical, public, protected :: & - CPFEM_init_done = .false. !< remember whether init has been done already + integer(pInt), public :: & + cycleCounter = 0_pInt !< needs description integer(pInt), parameter, public :: & CPFEM_CALCRESULTS = 2_pInt**0_pInt, & @@ -68,7 +58,6 @@ contains !-------------------------------------------------------------------------------------------------- subroutine CPFEM_initAll - CPFEM_init_done = .true. call DAMASK_interface_init call prec_init call IO_init @@ -149,9 +138,6 @@ subroutine CPFEM_general(mode, ffn, ffn1, temperature_inp, dt, elFE, ip, cauchyS write(6,'(/,a)') '#############################################' write(6,'(a1,a22,1x,i8,a13)') '#','element', elCP, '#' write(6,'(a1,a22,1x,i8,a13)') '#','ip', ip, '#' - write(6,'(a1,a22,1x,f15.7,a6)') '#','theTime', theTime, '#' - write(6,'(a1,a22,1x,f15.7,a6)') '#','theDelta', theDelta, '#' - write(6,'(a1,a22,1x,i8,a13)') '#','theInc', theInc, '#' write(6,'(a1,a22,1x,i8,a13)') '#','cycleCounter', cycleCounter, '#' write(6,'(a1,a22,1x,i8,a13)') '#','computationMode',mode, '#' if (terminallyIll) & diff --git a/src/DAMASK_marc.f90 b/src/DAMASK_marc.f90 index 5349608c9..efa054dbf 100644 --- a/src/DAMASK_marc.f90 +++ b/src/DAMASK_marc.f90 @@ -42,6 +42,7 @@ module DAMASK_interface logical, protected, public :: symmetricSolver character(len=*), parameter, public :: INPUTFILEEXTENSION = '.dat' + public :: & DAMASK_interface_init, & @@ -241,6 +242,17 @@ subroutine hypela2(d,g,e,de,s,t,dt,ngens,m,nn,kcus,matus,ndi,nshear,disp, & real(pReal), dimension(6,6) :: ddsdde integer :: computationMode, i, cp_en, node, CPnodeID integer(4) :: defaultNumThreadsInt !< default value set by Marc + + integer(pInt), save :: & + theInc = -1_pInt, & !< needs description + lastLovl = 0_pInt !< lovl in previous call to marc hypela2 + real(pReal), save :: & + theTime = 0.0_pReal, & !< needs description + theDelta = 0.0_pReal + logical, save :: & + lastIncConverged = .false., & !< needs description + outdatedByNewInc = .false., & !< needs description + CPFEM_init_done = .false. !< remember whether init has been done already if (iand(debug_level(debug_MARC),debug_LEVELBASIC) /= 0) then write(6,'(a,/,i8,i8,i2)') ' MSC.MARC information on shape of element(2), IP:', m, nn @@ -260,14 +272,17 @@ subroutine hypela2(d,g,e,de,s,t,dt,ngens,m,nn,kcus,matus,ndi,nshear,disp, & defaultNumThreadsInt = omp_get_num_threads() ! remember number of threads set by Marc call omp_set_num_threads(1) ! no openMP - if (.not. CPFEM_init_done) call CPFEM_initAll + if (.not. CPFEM_init_done) then + CPFEM_init_done = .true. + call CPFEM_initAll + endif computationMode = 0 ! save initialization value, since it does not result in any calculation if (lovl == 4 ) then ! jacobian requested by marc if (timinc < theDelta .and. theInc == inc .and. lastLovl /= lovl) & ! first after cutback computationMode = CPFEM_RESTOREJACOBIAN elseif (lovl == 6) then ! stress requested by marc - computationMode = CPFEM_CALCRESULTS ! always calc + computationMode = CPFEM_CALCRESULTS cp_en = mesh_FEM2DAMASK_elem(m(1)) if (cptim > theTime .or. inc /= theInc) then ! reached "convergence" terminallyIll = .false.