From ada92a9b7481eb944e31db45656ddc6ac3ac7466 Mon Sep 17 00:00:00 2001 From: Christoph Kords Date: Mon, 15 Jun 2009 13:11:21 +0000 Subject: [PATCH] added a new module called numerics.f90 which reads in all numerical "parameters" from the file numerics.config (also being added). From now on this file has to be located in the working directory of the FEM-model and has to contain all necessary parameters. --- trunk/CPFEM.f90 | 440 ++++++++++++++---------- trunk/IO.f90 | 176 ++++++---- trunk/constitutive_j2.f90 | 14 +- trunk/constitutive_phenomenological.f90 | 1 - trunk/crystallite.f90 | 22 +- trunk/debug.f90 | 55 ++- trunk/homogenization.f90 | 30 +- trunk/homogenization_isostrain.f90 | 12 +- trunk/lattice.f90 | 56 +-- trunk/material.f90 | 18 +- trunk/mesh.f90 | 181 +++++----- trunk/mpie_cpfem_marc.f90 | 39 ++- trunk/numerics.config | 18 + trunk/numerics.f90 | 169 +++++++++ trunk/prec.f90 | 24 +- 15 files changed, 782 insertions(+), 473 deletions(-) create mode 100644 trunk/numerics.config create mode 100644 trunk/numerics.f90 diff --git a/trunk/CPFEM.f90 b/trunk/CPFEM.f90 index 4ad3dd916..1420f964c 100644 --- a/trunk/CPFEM.f90 +++ b/trunk/CPFEM.f90 @@ -1,221 +1,297 @@ !############################################################## - MODULE CPFEM +MODULE CPFEM !############################################################## ! *** CPFEM engine *** ! - use prec, only: pReal,pInt - implicit none -! -! **************************************************************** -! *** General variables for the material behaviour calculation *** -! **************************************************************** - real(pReal), dimension (:,:,:), allocatable :: CPFEM_cs ! Cauchy stress - real(pReal), dimension (:,:,:,:), allocatable :: CPFEM_dcsdE ! Cauchy stress tangent - real(pReal), dimension (:,:,:,:), allocatable :: CPFEM_dcsdE_knownGood ! known good tangent +use prec, only: pReal, & + pInt +implicit none + +real(pReal), parameter :: CPFEM_odd_stress = 1e15_pReal, & + CPFEM_odd_jacobian = 1e50_pReal - logical :: CPFEM_init_done = .false. ! remember whether init has been done already - logical :: CPFEM_calc_done = .false. ! remember whether first IP has already calced the results +real(pReal), dimension (:,:,:), allocatable :: CPFEM_cs ! Cauchy stress +real(pReal), dimension (:,:,:,:), allocatable :: CPFEM_dcsdE ! Cauchy stress tangent +real(pReal), dimension (:,:,:,:), allocatable :: CPFEM_dcsdE_knownGood ! known good tangent + +logical :: CPFEM_init_done = .false., & ! remember whether init has been done already + CPFEM_calc_done = .false. ! remember whether first IP has already calced the results + + +CONTAINS - real(pReal), parameter :: CPFEM_odd_stress = 1e15_pReal, CPFEM_odd_jacobian = 1e50_pReal -! - CONTAINS -! !********************************************************* !*** allocate the arrays defined in module CPFEM *** !*** and initialize them *** !********************************************************* - SUBROUTINE CPFEM_init() +subroutine CPFEM_init() - use prec, only: pInt,pReal - use FEsolving, only: parallelExecution,symmetricSolver,FEsolving_execElem,FEsolving_execIP - use mesh, only: mesh_element,mesh_NcpElems,mesh_maxNips,FE_Nips - use material, only: homogenization_maxNgrains - use constitutive, only: constitutive_maxSizePostResults - use crystallite, only: crystallite_Nresults - use homogenization, only: homogenization_maxSizePostResults + use prec, only: pInt + use FEsolving, only: parallelExecution, & + symmetricSolver + use mesh, only: mesh_NcpElems, & + mesh_maxNips - implicit none - integer(pInt) e,i,g + implicit none - allocate(CPFEM_cs(6,mesh_maxNips,mesh_NcpElems)) ; CPFEM_cs = 0.0_pReal - allocate(CPFEM_dcsdE(6,6,mesh_maxNips,mesh_NcpElems)) ; CPFEM_dcsde = 0.0_pReal - allocate(CPFEM_dcsdE_knownGood(6,6,mesh_maxNips,mesh_NcpElems)) ; CPFEM_dcsde_knownGood = 0.0_pReal + ! initialize stress and jacobian to zero + allocate(CPFEM_cs(6,mesh_maxNips,mesh_NcpElems)) ; CPFEM_cs = 0.0_pReal + allocate(CPFEM_dcsdE(6,6,mesh_maxNips,mesh_NcpElems)) ; CPFEM_dcsde = 0.0_pReal + allocate(CPFEM_dcsdE_knownGood(6,6,mesh_maxNips,mesh_NcpElems)) ; CPFEM_dcsde_knownGood = 0.0_pReal + + !$OMP CRITICAL (write2out) + write(6,*) + write(6,*) '<<<+- cpfem init -+>>>' + write(6,*) + write(6,'(a32,x,6(i5,x))') 'CPFEM_cs: ', shape(CPFEM_cs) + write(6,'(a32,x,6(i5,x))') 'CPFEM_dcsde: ', shape(CPFEM_dcsde) + write(6,'(a32,x,6(i5,x))') 'CPFEM_dcsde_knownGood: ', shape(CPFEM_dcsde_knownGood) + write(6,*) + write(6,*) 'parallelExecution: ', parallelExecution + write(6,*) 'symmetricSolver: ', symmetricSolver + call flush(6) + !$OMP END CRITICAL (write2out) + return + +endsubroutine -! *** Output to MARC output file *** -!$OMP CRITICAL (write2out) - write(6,*) - write(6,*) '<<<+- cpfem init -+>>>' - write(6,*) - write(6,'(a32,x,6(i5,x))') 'CPFEM_cs: ', shape(CPFEM_cs) - write(6,'(a32,x,6(i5,x))') 'CPFEM_dcsde: ', shape(CPFEM_dcsde) - write(6,'(a32,x,6(i5,x))') 'CPFEM_dcsde_knownGood: ', shape(CPFEM_dcsde_knownGood) - write(6,*) - write(6,*) 'parallelExecution: ', parallelExecution - write(6,*) 'symmetricSolver: ', symmetricSolver - call flush(6) -!$OMP END CRITICAL (write2out) - return -! - END SUBROUTINE -! -! !*********************************************************************** !*** perform initialization at first call, update variables and *** !*** call the actual material model *** -! -! CPFEM_mode computation mode (regular, collection, recycle) -! ffn deformation gradient for t=t0 -! ffn1 deformation gradient for t=t1 -! Temperature temperature -! CPFEM_dt time increment -! CPFEM_en element number -! CPFEM_in intergration point number -! CPFEM_stress stress vector in Mandel notation -! CPFEM_updateJaco flag to initiate computation of Jacobian -! CPFEM_jaco jacobian in Mandel notation -! CPFEM_ngens size of stress strain law !*********************************************************************** -subroutine CPFEM_general(CPFEM_mode, ffn, ffn1, Temperature, CPFEM_dt,& - CPFEM_en, CPFEM_in, CPFEM_stress, CPFEM_updateJaco, CPFEM_jaco, CPFEM_ngens) -! note: CPFEM_stress = Cauchy stress cs(6) and CPFEM_jaco = Consistent tangent dcs/de -! - use prec, only: pReal,pInt - use FEsolving - use debug - use math - use mesh, only: mesh_init,& - mesh_FEasCP,mesh_element,mesh_NcpElems,mesh_maxNips,FE_Nips - use lattice, only: lattice_init - use material, only: material_init, homogenization_maxNgrains - use constitutive, only: constitutive_init,& - constitutive_state0,constitutive_state - use crystallite - use homogenization - implicit none +subroutine CPFEM_general(mode, ffn, ffn1, Temperature, dt, element, IP, cauchyStress, jacobian, ngens) + ! note: cauchyStress = Cauchy stress cs(6) and jacobian = Consistent tangent dcs/de - integer(pInt) CPFEM_en, CPFEM_in, cp_en, CPFEM_ngens, i,j,k,l,m,n - real(pReal), dimension (3,3) :: ffn,ffn1,Kirchhoff - real(pReal), dimension (3,3,3,3) :: H, H_sym - real(pReal), dimension(CPFEM_ngens) :: CPFEM_stress - real(pReal), dimension(CPFEM_ngens,CPFEM_ngens) :: CPFEM_jaco - real(pReal) Temperature,CPFEM_dt,J_inverse - integer(pInt) CPFEM_mode ! 1: regular computation with aged results& - ! 2: regular computation& - ! 3: collection of FEM data& - ! 4: recycling of former results (MARC speciality)& - ! 5: record tangent from former converged inc& - ! 6: restore tangent from former converged inc - integer(pInt) e - logical CPFEM_updateJaco + !*** variables and functions from other modules ***! + use prec, only: pReal, & + pInt + use numerics, only: numerics_init, & + relevantStrain, & + iJacoStiffness + use debug, only: debug_init + use FEsolving, only: FE_init, & + parallelExecution, & + outdatedFFN1, & + cycleCounter, & + theInc, & + theCycle, & + theLovl, & + theTime, & + FEsolving_execElem, & + FEsolving_execIP + use math, only: math_init, & + math_identity2nd, & + math_mul33x33, & + math_det3x3, & + math_I3, & + math_Mandel3333to66, & + math_Mandel33to6 + use mesh, only: mesh_init, & + mesh_FEasCP, & + mesh_NcpElems, & + mesh_maxNips + use lattice, only: lattice_init + use material, only: material_init, & + homogenization_maxNgrains + use constitutive, only: constitutive_init,& + constitutive_state0,constitutive_state + use crystallite, only: crystallite_init, & + crystallite_F0, & + crystallite_partionedF, & + crystallite_Fp0, & + crystallite_Fp, & + crystallite_Lp0, & + crystallite_Lp + use homogenization, only: homogenization_init, & + homogenization_sizeState, & + homogenization_state, & + homogenization_state0, & + materialpoint_F, & + materialpoint_F0, & + materialpoint_P, & + materialpoint_dPdF, & + materialpoint_Temperature, & + materialpoint_stressAndItsTangent, & + materialpoint_postResults + + implicit none + + !*** input variables ***! + integer(pInt), intent(in) :: element, & ! FE element number + IP, & ! FE integration point number + ngens ! size of stress strain law + real(pReal), intent(in) :: Temperature, & ! temperature + dt ! time increment + real(pReal), dimension (3,3), intent(in) :: ffn, & ! deformation gradient for t=t0 + ffn1 ! deformation gradient for t=t1 + integer(pInt), intent(in) :: mode ! computation mode 1: regular computation with aged results + ! 2: regular computation + ! 3: collection of FEM data + ! 4: recycling of former results (MARC speciality) + ! 5: record tangent from former converged inc + ! 6: restore tangent from former converged inc + + !*** output variables ***! + real(pReal), dimension(ngens), intent(out) :: cauchyStress ! stress vector in Mandel notation + real(pReal), dimension(ngens,ngens), intent(out) :: jacobian ! jacobian in Mandel notation + + !*** local variables ***! + real(pReal) J_inverse ! inverse of Jacobian + real(pReal), dimension (3,3) :: Kirchhoff + real(pReal), dimension (3,3,3,3) :: H, & + H_sym + integer(pInt) cp_en, & ! crystal plasticity element number + i, & + j, & + k, & + l, & + m, & + n + logical updateJaco ! flag indicating if JAcobian has to be updated + + !*** global variables ***! + ! CPFEM_cs, & + ! CPFEM_dcsdE, & + ! CPFEM_dcsdE_knownGood, & + ! CPFEM_init_done, & + ! CPFEM_calc_done, & + ! CPFEM_odd_stress, & + ! CPFEM_odd_jacobian + - if (.not. CPFEM_init_done) then ! initialization step (three dimensional stress state check missing?) - call math_init() - call FE_init() - call mesh_init() + ! initialization step (three dimensional stress state check missing?) + if (.not. CPFEM_init_done) then + call numerics_init() + call debug_init() + call math_init() + call FE_init() + call mesh_init() + call lattice_init() + call material_init() + call constitutive_init() + call crystallite_init() + call homogenization_init() + call CPFEM_init() + CPFEM_init_done = .true. + endif - FEsolving_execElem = (/1,mesh_NcpElems/) - allocate(FEsolving_execIP(2,mesh_NcpElems)); FEsolving_execIP = 1_pInt - forall (e = 1:mesh_NcpElems) FEsolving_execIP(2,e) = FE_Nips(mesh_element(2,e)) - - call lattice_init() - call material_init() - call constitutive_init() - call crystallite_init() - call homogenization_init() - call CPFEM_init() - CPFEM_init_done = .true. - endif - - cp_en = mesh_FEasCP('elem',CPFEM_en) - if (cp_en == 1 .and. CPFEM_in == 1) then + cp_en = mesh_FEasCP('elem',element) + + if (cp_en == 1 .and. IP == 1) then write(6,*) '#####################################' write(6,'(a10,1x,f8.4,1x,a10,1x,i4,1x,a10,1x,i3,1x,a10,1x,i2,x,a10,1x,i2)') & 'theTime',theTime,'theInc',theInc,'theCycle',theCycle,'theLovl',theLovl,& - 'mode',CPFEM_mode + 'mode',mode write(6,*) '#####################################' - endif + endif - select case (CPFEM_mode) - case (1,2) ! regular computation (with aging of results if mode == 1) - if (CPFEM_mode == 1) then ! age results at start of new increment - crystallite_F0 = crystallite_partionedF ! crystallite deformation (_subF is perturbed...) - crystallite_Fp0 = crystallite_Fp ! crystallite plastic deformation - crystallite_Lp0 = crystallite_Lp ! crystallite plastic velocity - forall (i = 1:homogenization_maxNgrains,& + ! according to our "mode" we decide what to do + select case (mode) + + ! --+>> REGULAR COMPUTATION (WITH AGING OF RESULTS IF MODE == 1) <<+-- + case (1,2) + ! age results if mode == 1 + if (mode == 1) then + crystallite_F0 = crystallite_partionedF ! crystallite deformation (_subF is perturbed...) + crystallite_Fp0 = crystallite_Fp ! crystallite plastic deformation + crystallite_Lp0 = crystallite_Lp ! crystallite plastic velocity + forall ( i = 1:homogenization_maxNgrains, & j = 1:mesh_maxNips, & - k = 1:mesh_NcpElems) & - constitutive_state0(i,j,k)%p = constitutive_state(i,j,k)%p ! microstructure of crystallites - write(6,'(a10,/,4(3(f10.3,x),/))') 'aged state',constitutive_state(1,1,1)%p/1e6 - do j = 1,mesh_maxNips - do k = 1,mesh_NcpElems - if (homogenization_sizeState(j,k) > 0_pInt) & - homogenization_state0(j,k)%p = homogenization_state(j,k)%p ! internal state of homogenization scheme - enddo - enddo - endif - - if (outdatedFFN1 .or. any(abs(ffn1 - materialpoint_F(:,:,CPFEM_in,cp_en)) > relevantStrain)) then - if (.not. outdatedFFN1) write(6,'(a11,x,i5,x,i2,x,a10,/,3(3(f10.3,x),/))') 'outdated at',cp_en,CPFEM_in,'FFN1 now:',ffn1(:,1),ffn1(:,2),ffn1(:,3) - outdatedFFN1 = .true. - CPFEM_cs(1:CPFEM_ngens,CPFEM_in,cp_en) = CPFEM_odd_stress - CPFEM_dcsde(1:CPFEM_ngens,1:CPFEM_ngens,CPFEM_in,cp_en) = CPFEM_odd_jacobian*math_identity2nd(CPFEM_ngens) - else - if (.not. parallelExecution) then - FEsolving_execElem(1) = cp_en - FEsolving_execElem(2) = cp_en - FEsolving_execIP(1,cp_en) = CPFEM_in - FEsolving_execIP(2,cp_en) = CPFEM_in - call materialpoint_stressAndItsTangent(CPFEM_updateJaco, CPFEM_dt) - call materialpoint_postResults(CPFEM_dt) - elseif (.not. CPFEM_calc_done) then - call materialpoint_stressAndItsTangent(CPFEM_updateJaco, CPFEM_dt) ! parallel execution inside - call materialpoint_postResults(CPFEM_dt) - CPFEM_calc_done = .true. - endif - -! translate from P and dP/dF to CS and dCS/dE - Kirchhoff = math_mul33x33(materialpoint_P(:,:,CPFEM_in, cp_en),transpose(materialpoint_F(:,:,CPFEM_in, cp_en))) - J_inverse = 1.0_pReal/math_det3x3(materialpoint_F(:,:,CPFEM_in, cp_en)) - CPFEM_cs(1:CPFEM_ngens,CPFEM_in,cp_en) = math_Mandel33to6(J_inverse*Kirchhoff) - - H = 0.0_pReal - forall(i=1:3,j=1:3,k=1:3,l=1:3,m=1:3,n=1:3) & - H(i,j,k,l) = H(i,j,k,l) + & - materialpoint_F(j,m,CPFEM_in,cp_en) * & - materialpoint_F(l,n,CPFEM_in,cp_en) * & - materialpoint_dPdF(i,m,k,n,CPFEM_in,cp_en) - & - math_I3(j,l)*materialpoint_F(i,m,CPFEM_in,cp_en)*materialpoint_P(k,m,CPFEM_in,cp_en) + & - 0.5_pReal*(math_I3(i,k)*Kirchhoff(j,l) + math_I3(j,l)*Kirchhoff(i,k) + & - math_I3(i,l)*Kirchhoff(j,k) + math_I3(j,k)*Kirchhoff(i,l)) - forall(i=1:3,j=1:3,k=1:3,l=1:3) & - H_sym(i,j,k,l)= 0.25_pReal*(H(i,j,k,l)+H(j,i,k,l)+H(i,j,l,k)+H(j,i,l,k)) ! where to use the symmetric version?? - CPFEM_dcsde(1:CPFEM_ngens,1:CPFEM_ngens,CPFEM_in,cp_en) = math_Mandel3333to66(J_inverse*H) + k = 1:mesh_NcpElems ) & + constitutive_state0(i,j,k)%p = constitutive_state(i,j,k)%p ! microstructure of crystallites + write(6,'(a10,/,4(3(f10.3,x),/))') 'aged state',constitutive_state(1,1,1)%p/1e6 + do j = 1,mesh_maxNips + do k = 1,mesh_NcpElems + if (homogenization_sizeState(j,k) > 0_pInt) & + homogenization_state0(j,k)%p = homogenization_state(j,k)%p ! internal state of homogenization scheme + enddo + enddo endif - case (3) ! collect and return odd result - materialpoint_Temperature(CPFEM_in,cp_en) = Temperature - materialpoint_F0(:,:,CPFEM_in,cp_en) = ffn - materialpoint_F(:,:,CPFEM_in,cp_en) = ffn1 - CPFEM_cs(1:CPFEM_ngens,CPFEM_in,cp_en) = CPFEM_odd_stress - CPFEM_dcsde(1:CPFEM_ngens,1:CPFEM_ngens,CPFEM_in,cp_en) = CPFEM_odd_jacobian*math_identity2nd(CPFEM_ngens) - CPFEM_calc_done = .false. - case (4) ! do nothing since we can recycle the former results (MARC specialty) - case (5) ! record consistent tangent at beginning of new FE increment (while recycling) + ! deformation gradient outdated or any actual deformation gradient differs more than relevantStrain from the stored one + if (outdatedFFN1 .or. any(abs(ffn1 - materialpoint_F(:,:,IP,cp_en)) > relevantStrain)) then + if (.not. outdatedFFN1) & + write(6,'(a11,x,i5,x,i2,x,a10,/,3(3(f10.3,x),/))') 'outdated at',cp_en,IP,'FFN1 now:',ffn1(:,1),ffn1(:,2),ffn1(:,3) + outdatedFFN1 = .true. + CPFEM_cs(1:ngens,IP,cp_en) = CPFEM_odd_stress + CPFEM_dcsde(1:ngens,1:ngens,IP,cp_en) = CPFEM_odd_jacobian*math_identity2nd(ngens) + + ! deformation gradient is not outdated + else + ! set flag for Jacobian update + updateJaco = (mod(cycleCounter-4,4_pInt*iJacoStiffness)==0) + + ! no parallel computation + if (.not. parallelExecution) then + ! we just take one single element and IP + FEsolving_execElem(1) = cp_en + FEsolving_execElem(2) = cp_en + FEsolving_execIP(1,cp_en) = IP + FEsolving_execIP(2,cp_en) = IP + call materialpoint_stressAndItsTangent(updateJaco, dt) ! calculate stress and its tangent + call materialpoint_postResults(dt) ! post results + + ! parallel computation and calulation not yet done + 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 + CPFEM_calc_done = .true. + endif + + ! translate from P to CS + Kirchhoff = math_mul33x33(materialpoint_P(:,:,IP, cp_en),transpose(materialpoint_F(:,:,IP, cp_en))) + J_inverse = 1.0_pReal/math_det3x3(materialpoint_F(:,:,IP, cp_en)) + CPFEM_cs(1:ngens,IP,cp_en) = math_Mandel33to6(J_inverse*Kirchhoff) + + ! translate from dP/dF to dCS/dE + H = 0.0_pReal + forall(i=1:3,j=1:3,k=1:3,l=1:3,m=1:3,n=1:3) & + H(i,j,k,l) = H(i,j,k,l) + & + materialpoint_F(j,m,IP,cp_en) * & + materialpoint_F(l,n,IP,cp_en) * & + materialpoint_dPdF(i,m,k,n,IP,cp_en) - & + math_I3(j,l)*materialpoint_F(i,m,IP,cp_en)*materialpoint_P(k,m,IP,cp_en) + & + 0.5_pReal*(math_I3(i,k)*Kirchhoff(j,l) + math_I3(j,l)*Kirchhoff(i,k) + & + math_I3(i,l)*Kirchhoff(j,k) + math_I3(j,k)*Kirchhoff(i,l)) + forall(i=1:3,j=1:3,k=1:3,l=1:3) & + H_sym(i,j,k,l)= 0.25_pReal*(H(i,j,k,l)+H(j,i,k,l)+H(i,j,l,k)+H(j,i,l,k)) ! where to use the symmetric version?? + CPFEM_dcsde(1:ngens,1:ngens,IP,cp_en) = math_Mandel3333to66(J_inverse*H) + + endif + + ! --+>> COLLECTION OF FEM DATA AND RETURN OF ODD STRESS AND JACOBIAN <<+-- + case (3) + materialpoint_Temperature(IP,cp_en) = Temperature + materialpoint_F0(:,:,IP,cp_en) = ffn + materialpoint_F(:,:,IP,cp_en) = ffn1 + CPFEM_cs(1:ngens,IP,cp_en) = CPFEM_odd_stress + CPFEM_dcsde(1:ngens,1:ngens,IP,cp_en) = CPFEM_odd_jacobian*math_identity2nd(ngens) + CPFEM_calc_done = .false. + + ! --+>> RECYCLING OF FORMER RESULTS (MARC SPECIALTY) <<+-- + case (4) + ! do nothing + + ! --+>> RECORD JACOBIAN FROM FORMER CONVERGED INC <<+-- + case (5) CPFEM_dcsde_knownGood = CPFEM_dcsde - case (6) ! restore consistent tangent after FE cutback + + ! --+>> RESTORE CONSISTENT JACOBIAN FROM FORMER CONVERGED INC <<+-- + case (6) CPFEM_dcsde = CPFEM_dcsde_knownGood - end select + + end select -! return the local stress and the jacobian from storage - CPFEM_stress(1:CPFEM_ngens) = CPFEM_cs(1:CPFEM_ngens,CPFEM_in,cp_en) - CPFEM_jaco(1:CPFEM_ngens,1:CPFEM_ngens) = CPFEM_dcsdE(1:CPFEM_ngens,1:CPFEM_ngens,CPFEM_in,cp_en) + ! return the local stress and the jacobian from storage + cauchyStress(1:ngens) = CPFEM_cs(1:ngens,IP,cp_en) + jacobian(1:ngens,1:ngens) = CPFEM_dcsdE(1:ngens,1:ngens,IP,cp_en) - return + return end subroutine - END MODULE -!############################################################## \ No newline at end of file + END MODULE CPFEM + \ No newline at end of file diff --git a/trunk/IO.f90 b/trunk/IO.f90 index 88521970f..35858b26d 100644 --- a/trunk/IO.f90 +++ b/trunk/IO.f90 @@ -30,7 +30,7 @@ ! open existing file to given unit ! path to file is relative to working directory !******************************************************************** - logical FUNCTION IO_open_file(unit,relPath) + logical function IO_open_file(unit,relPath) use prec, only: pInt implicit none @@ -47,13 +47,13 @@ 100 IO_open_file = .false. return - END FUNCTION + endfunction !******************************************************************** ! open FEM inputfile to given unit !******************************************************************** - logical FUNCTION IO_open_inputFile(unit) + logical function IO_open_inputFile(unit) use prec, only: pReal, pInt implicit none @@ -68,20 +68,20 @@ ext='dat' ! MARC else ext='inp' ! ABAQUS - end if + endif open(unit,status='old',err=100,file=outName(1:extPos-1)//ext) IO_open_inputFile = .true. return 100 IO_open_inputFile = .false. return - END FUNCTION + endfunction !******************************************************************** ! hybrid IA repetition counter !******************************************************************** - FUNCTION hybridIA_reps(dV_V,steps,C) + function hybridIA_reps(dV_V,steps,C) use prec, only: pReal, pInt implicit none @@ -96,18 +96,18 @@ do Phi =1,steps(2) do phi2=1,steps(3) hybridIA_reps = hybridIA_reps+nint(C*dV_V(phi2,Phi,phi1), pInt) - end do - end do - end do + enddo + enddo + enddo return - END FUNCTION + endfunction !******************************************************************** ! hybrid IA sampling of ODFfile !******************************************************************** - FUNCTION IO_hybridIA(Nast,ODFfileName) + function IO_hybridIA(Nast,ODFfileName) use prec, only: pReal, pInt use math, only: inRad @@ -135,7 +135,7 @@ if (pos(1).ne.3) goto 100 do i=1,3 limits(i) = IO_intValue(line,pos,i)*inRad - end do + enddo !--- deltas in phi1, Phi, phi2 --- read(999,fmt=fileFormat,end=100) line @@ -143,7 +143,7 @@ if (pos(1).ne.3) goto 100 do i=1,3 deltas(i) = IO_intValue(line,pos,i)*inRad - end do + enddo steps = nint(limits/deltas,pInt) allocate(dV_V(steps(3),steps(2),steps(1))) @@ -153,7 +153,7 @@ center = 0.5_pReal else center = 0.0_pReal - end if + endif !--- skip blank line --- read(999,fmt=fileFormat,end=100) line @@ -172,11 +172,11 @@ sum_dV_V = sum_dV_V+prob else prob = 0.0_pReal - end if + endif dV_V(phi2,Phi,phi1) = prob*dg_0*sin((Phi-1.0_pReal+center)*deltas(2)) - end do - end do - end do + enddo + enddo + enddo dV_V = dV_V/sum_dV_V ! normalize to 1 @@ -188,7 +188,7 @@ do while (hybridIA_reps(dV_V,steps,upperC) < Nset) lowerC = upperC upperC = upperC*2.0_pReal - end do + enddo !--- binary search for best C --- do C = (upperC+lowerC)/2.0_pReal @@ -203,8 +203,8 @@ upperC = C else exit - end if - end do + endif + enddo allocate(binSet(Nreps)) bin = 0 ! bin counter @@ -216,9 +216,9 @@ binSet(i:i+reps-1) = bin bin = bin+1 ! advance bin i = i+reps ! advance set - end do - end do - end do + enddo + enddo + enddo do i=1,Nast if (i < Nast) then @@ -226,13 +226,13 @@ j = nint(rnd*(Nast-i)+i+0.5_pReal,pInt) else j = i - end if + endif bin = binSet(j) IO_hybridIA(1,i) = deltas(1)*(mod(bin/(steps(3)*steps(2)),steps(1))+center) ! phi1 IO_hybridIA(2,i) = deltas(2)*(mod(bin/ steps(3) ,steps(2))+center) ! Phi IO_hybridIA(3,i) = deltas(3)*(mod(bin ,steps(3))+center) ! phi2 binSet(j) = binSet(i) - end do + enddo close(999) return @@ -241,13 +241,13 @@ close(999) return - END FUNCTION + endfunction !******************************************************************** ! identifies lines without content !******************************************************************** - PURE FUNCTION IO_isBlank (line) + pure function IO_isBlank (line) use prec, only: pInt implicit none @@ -264,12 +264,12 @@ return - END FUNCTION + endfunction !******************************************************************** ! get tagged content of line !******************************************************************** - PURE FUNCTION IO_getTag (line,openChar,closechar) + pure function IO_getTag (line,openChar,closechar) use prec, only: pInt implicit none @@ -288,11 +288,11 @@ return - END FUNCTION + endfunction !********************************************************************* - FUNCTION IO_countSections(file,part) + function IO_countSections(file,part) !********************************************************************* use prec, only: pInt implicit none @@ -321,13 +321,13 @@ 100 return - END FUNCTION + endfunction !********************************************************************* ! return array of myTag counts within for at most N[sections] !********************************************************************* - FUNCTION IO_countTagInPart(file,part,myTag,Nsections) + function IO_countTagInPart(file,part,myTag,Nsections) use prec, only: pInt implicit none @@ -367,13 +367,13 @@ 100 IO_countTagInPart = counter return -END FUNCTION +endfunction !********************************************************************* ! return array of myTag presence within for at most N[sections] !********************************************************************* - FUNCTION IO_spotTagInPart(file,part,myTag,Nsections) + function IO_spotTagInPart(file,part,myTag,Nsections) use prec, only: pInt implicit none @@ -412,7 +412,7 @@ END FUNCTION 100 return -END FUNCTION +endfunction !******************************************************************** @@ -420,7 +420,7 @@ END FUNCTION ! return array containing number of parts found and ! their left/right positions to be used by IO_xxxVal !******************************************************************** - PURE FUNCTION IO_stringPos (line,N) + pure function IO_stringPos (line,N) use prec, only: pReal,pInt implicit none @@ -438,17 +438,17 @@ END FUNCTION IO_stringPos(part*2) = IO_stringPos(part*2-1)+verify(line(IO_stringPos(part*2-1)+1:),sep) IO_stringPos(part*2+1) = IO_stringPos(part*2)+scan(line(IO_stringPos(part*2):),sep)-2 part = part+1 - end do + enddo IO_stringPos(1) = part-1 return - END FUNCTION + endfunction !******************************************************************** ! read string value at pos from line !******************************************************************** - PURE FUNCTION IO_stringValue (line,positions,pos) + pure function IO_stringValue (line,positions,pos) use prec, only: pReal,pInt implicit none @@ -464,13 +464,13 @@ END FUNCTION endif return - END FUNCTION + endfunction !******************************************************************** ! read string value at pos from fixed format line !******************************************************************** - PURE FUNCTION IO_fixedStringValue (line,ends,pos) + pure function IO_fixedStringValue (line,ends,pos) use prec, only: pReal,pInt implicit none @@ -482,13 +482,13 @@ END FUNCTION IO_fixedStringValue = line(ends(pos)+1:ends(pos+1)) return - END FUNCTION + endfunction !******************************************************************** ! read float value at pos from line !******************************************************************** - PURE FUNCTION IO_floatValue (line,positions,pos) + pure function IO_floatValue (line,positions,pos) use prec, only: pReal,pInt implicit none @@ -504,13 +504,13 @@ END FUNCTION 100 IO_floatValue = huge(1.0_pReal) return - END FUNCTION + endfunction !******************************************************************** ! read float value at pos from fixed format line !******************************************************************** - PURE FUNCTION IO_fixedFloatValue (line,ends,pos) + pure function IO_fixedFloatValue (line,ends,pos) use prec, only: pReal,pInt implicit none @@ -524,13 +524,13 @@ END FUNCTION 100 IO_fixedFloatValue = huge(1.0_pReal) return - END FUNCTION + endfunction !******************************************************************** ! read float x.y+z value at pos from format line line !******************************************************************** - PURE FUNCTION IO_fixedNoEFloatValue (line,ends,pos) + pure function IO_fixedNoEFloatValue (line,ends,pos) use prec, only: pReal,pInt implicit none @@ -553,13 +553,13 @@ END FUNCTION 100 IO_fixedNoEFloatValue = huge(1.0_pReal) return - END FUNCTION + endfunction !******************************************************************** ! read int value at pos from line !******************************************************************** - PURE FUNCTION IO_intValue (line,positions,pos) + pure function IO_intValue (line,positions,pos) use prec, only: pReal,pInt implicit none @@ -575,13 +575,13 @@ END FUNCTION 100 IO_intValue = huge(1_pInt) return - END FUNCTION + endfunction !******************************************************************** ! read int value at pos from fixed format line !******************************************************************** - PURE FUNCTION IO_fixedIntValue (line,ends,pos) + pure function IO_fixedIntValue (line,ends,pos) use prec, only: pReal,pInt implicit none @@ -595,13 +595,13 @@ END FUNCTION 100 IO_fixedIntValue = huge(1_pInt) return - END FUNCTION + endfunction !******************************************************************** ! change character in line to lower case !******************************************************************** - PURE FUNCTION IO_lc (line) + pure function IO_lc (line) use prec, only: pInt implicit none @@ -616,13 +616,13 @@ END FUNCTION enddo return - END FUNCTION + endfunction !******************************************************************** ! in place change of character in line to lower case !******************************************************************** - SUBROUTINE IO_lcInplace (line) + subroutine IO_lcInplace (line) use prec, only: pInt implicit none @@ -638,13 +638,13 @@ END FUNCTION line = IO_lc return - END SUBROUTINE + endsubroutine !******************************************************************** ! read on in file to skip (at least) N chunks (may be over multiple lines) !******************************************************************** - SUBROUTINE IO_skipChunks (unit,N) + subroutine IO_skipChunks (unit,N) use prec, only: pReal,pInt implicit none @@ -659,17 +659,17 @@ END FUNCTION read(unit,'(A300)',end=100) line pos = IO_stringPos(line,maxNchunks) remainingChunks = remainingChunks - pos(1) - end do + enddo 100 return - END SUBROUTINE + endsubroutine !******************************************************************** ! count items in consecutive lines of ints concatenated by "c" ! as last char or range of values a "to" b !******************************************************************** - FUNCTION IO_countContinousIntValues (unit) + function IO_countContinousIntValues (unit) use prec, only: pReal,pInt implicit none @@ -695,13 +695,13 @@ END FUNCTION enddo 100 return - END FUNCTION + endfunction !********************************************************************* ! read consecutive lines of ints concatenated by "c" as last char ! or range of values a "to" b !********************************************************************* - FUNCTION IO_continousIntValues (unit,maxN,lookupName,lookupMap,lookupMaxN) + function IO_continousIntValues (unit,maxN,lookupName,lookupMap,lookupMaxN) use prec, only: pReal,pInt implicit none @@ -746,7 +746,7 @@ END FUNCTION enddo 100 return - END FUNCTION + endfunction !******************************************************************** @@ -754,11 +754,9 @@ END FUNCTION ! and terminate the Marc run with exit #9xxx ! in ABAQUS either time step is reduced or execution terminated !******************************************************************** - SUBROUTINE IO_error(ID,e,i,g,ext_msg) + subroutine IO_error(ID,e,i,g,ext_msg) use prec, only: pInt - - use debug implicit none integer(pInt), intent(in) :: ID @@ -819,6 +817,38 @@ END FUNCTION msg = 'Negative diffusion constant' case (240) msg = 'Non-positive Taylor factor' + case (260) + msg = 'Non-positive relevant strain' + case (261) + msg = 'Frequency for Stiffness update smaller than zero' + case (262) + msg = 'Frequency for Jacobian update of Lp residuum smaller than zero' + case (263) + msg = 'Non-positive perturbation value' + case (264) + msg = 'Limit for homogenization loop too small' + case (265) + msg = 'Limit for crystallite loop too small' + case (266) + msg = 'Limit for state loop too small' + case (267) + msg = 'Limit for stress loop too small' + case (268) + msg = 'Non-positive minimum substep size' + case (269) + msg = 'Non-positive relative tolerance for state' + case (270) + msg = 'Non-positive relative tolerance for stress' + case (271) + msg = 'Non-positive absolute tolerance for stress' + case (272) + msg = 'Non-positive relative tolerance of residual in GIA iteration' + case (273) + msg = 'Non-positive absolute tolerance of residual in GIA iteration' + case (274) + msg = 'Non-positive relative maximum value (upper bound) for GIA residual' + case (275) + msg = 'Limit for GIA iteration too small' case (300) msg = 'This material can only be used with elements with three direct stress components' case (500) @@ -850,8 +880,7 @@ END FUNCTION endif endif write(6,'(a38)') '+------------------------------------+' - - call debug_info() + call flush(6) call quit(9000+ID) !$OMP END CRITICAL (write2out) @@ -859,16 +888,15 @@ END FUNCTION ! ABAQUS returns in some cases return - END SUBROUTINE + endsubroutine !******************************************************************** ! write warning statements to standard out !******************************************************************** - SUBROUTINE IO_warning(ID,e,i,g,ext_msg) + subroutine IO_warning(ID,e,i,g,ext_msg) use prec, only: pInt - use debug implicit none integer(pInt), intent(in) :: ID @@ -901,6 +929,6 @@ END FUNCTION endif write(6,'(a38)') '+------------------------------------+' - END SUBROUTINE + endsubroutine END MODULE IO diff --git a/trunk/constitutive_j2.f90 b/trunk/constitutive_j2.f90 index 092f3ffec..1dd57dc30 100644 --- a/trunk/constitutive_j2.f90 +++ b/trunk/constitutive_j2.f90 @@ -176,7 +176,7 @@ subroutine constitutive_j2_init(file) return -end subroutine +endsubroutine function constitutive_j2_stateInit(myInstance) @@ -193,7 +193,7 @@ function constitutive_j2_stateInit(myInstance) constitutive_j2_stateInit = constitutive_j2_s0(myInstance) return -end function +endfunction function constitutive_j2_homogenizedC(state,ipc,ip,el) @@ -221,7 +221,7 @@ function constitutive_j2_homogenizedC(state,ipc,ip,el) return -end function +endfunction subroutine constitutive_j2_microstructure(Temperature,state,ipc,ip,el) @@ -245,7 +245,7 @@ subroutine constitutive_j2_microstructure(Temperature,state,ipc,ip,el) matID = phase_constitutionInstance(material_phase(ipc,ip,el)) -end subroutine +endsubroutine subroutine constitutive_j2_LpAndItsTangent(Lp,dLp_dTstar,Tstar_v,Temperature,state,ipc,ip,el) @@ -310,7 +310,7 @@ subroutine constitutive_j2_LpAndItsTangent(Lp,dLp_dTstar,Tstar_v,Temperature,sta return -end subroutine +endsubroutine function constitutive_j2_dotState(Tstar_v,Temperature,state,ipc,ip,el) @@ -351,7 +351,7 @@ function constitutive_j2_dotState(Tstar_v,Temperature,state,ipc,ip,el) return -end function +endfunction pure function constitutive_j2_postResults(Tstar_v,Temperature,dt,state,ipc,ip,el) @@ -401,6 +401,6 @@ pure function constitutive_j2_postResults(Tstar_v,Temperature,dt,state,ipc,ip,el return -end function +endfunction END MODULE diff --git a/trunk/constitutive_phenomenological.f90 b/trunk/constitutive_phenomenological.f90 index e4f131358..28c67fc92 100644 --- a/trunk/constitutive_phenomenological.f90 +++ b/trunk/constitutive_phenomenological.f90 @@ -180,7 +180,6 @@ subroutine constitutive_phenomenological_init(file) 100 do i = 1,maxNinstance constitutive_phenomenological_structure(i) = lattice_initializeStructure(constitutive_phenomenological_structureName(i), & - constitutive_phenomenological_CoverA(i)) ! sanity checks if (constitutive_phenomenological_structure(i) < 1 .or. & constitutive_phenomenological_structure(i) > 3) call IO_error(201) diff --git a/trunk/crystallite.f90 b/trunk/crystallite.f90 index ff136a83f..cd5bc0b60 100644 --- a/trunk/crystallite.f90 +++ b/trunk/crystallite.f90 @@ -229,8 +229,8 @@ subroutine crystallite_stressAndItsTangent(updateJaco) !*** variables and functions from other modules ***! use prec, only: pInt, & - pReal, & - subStepMin, & + pReal + use numerics, only: subStepMin, & pert_Fg, & nState, & nCryst @@ -325,7 +325,7 @@ subroutine crystallite_stressAndItsTangent(updateJaco) write (6,'(a,/,3(3(f12.7,x)/))') 'crystallite_partionedF of 1 1 1',crystallite_partionedF(1:3,:,1,1,1) write (6,'(a,/,3(3(f12.7,x)/))') 'crystallite_partionedLp0 of 1 1 1',crystallite_partionedLp0(1:3,:,1,1,1) - + !$OMP PARALLEL DO do e = FEsolving_execElem(1),FEsolving_execElem(2) ! iterate over elements to be processed myNgrains = homogenization_Ngrains(mesh_element(3,e)) @@ -346,6 +346,7 @@ subroutine crystallite_stressAndItsTangent(updateJaco) enddo enddo !$OMPEND PARALLEL DO + ! --+>> crystallite loop <<+-- @@ -458,7 +459,6 @@ subroutine crystallite_stressAndItsTangent(updateJaco) ! constitutive_state is internally interpolated with .._subState0 ! 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 myNgrains = homogenization_Ngrains(mesh_element(3,e)) @@ -610,7 +610,7 @@ subroutine crystallite_stressAndItsTangent(updateJaco) enddo !$OMPEND PARALLEL DO - if (debugger) write (6,*) 'Stiffness calculation finished' + if (debugger) write (6,*) 'Stiffness calculation finished' endif @@ -631,8 +631,8 @@ endsubroutine !*** variables and functions from other modules ***! use prec, only: pReal, & pInt, & - pLongInt, & - rTol_crystalliteState + pLongInt + use numerics, only: rTol_crystalliteState use constitutive, only: constitutive_dotState, & constitutive_sizeDotState, & constitutive_subState0, & @@ -708,12 +708,12 @@ endsubroutine !*** variables and functions from other modules ***! use prec, only: pReal, & pInt, & - pLongInt, & - nStress, & + pLongInt + use numerics, only: nStress, & aTol_crystalliteStress, & rTol_crystalliteStress, & - relevantStrain, & - iJacoLpresiduum + iJacoLpresiduum, & + relevantStrain use debug, only: debugger, & debug_cumLpCalls, & debug_cumLpTicks, & diff --git a/trunk/debug.f90 b/trunk/debug.f90 index a9fedf546..af4dbce09 100644 --- a/trunk/debug.f90 +++ b/trunk/debug.f90 @@ -4,12 +4,11 @@ !############################################################## use prec - implicit none - integer(pInt), dimension(nStress) :: debug_StressLoopDistribution = 0_pInt - integer(pInt), dimension(nState) :: debug_StateLoopDistribution = 0_pInt - integer(pInt), dimension(nState) :: debug_StiffnessStateLoopDistribution = 0_pInt - integer(pInt), dimension(nCryst) :: debug_CrystalliteLoopDistribution = 0_pInt + integer(pInt), dimension(:), allocatable :: debug_StressLoopDistribution + integer(pInt), dimension(:), allocatable :: debug_StateLoopDistribution + integer(pInt), dimension(:), allocatable :: debug_StiffnessStateLoopDistribution + integer(pInt), dimension(:), allocatable :: debug_CrystalliteLoopDistribution integer(pLongInt) :: debug_cumLpTicks = 0_pInt integer(pLongInt) :: debug_cumDotStateTicks = 0_pInt integer(pInt) :: debug_cumLpCalls = 0_pInt @@ -19,35 +18,53 @@ CONTAINS +subroutine debug_init() + + use prec, only: pInt + use numerics, only: nStress, & + nState, & + nCryst + implicit none + + allocate(debug_StressLoopDistribution(nStress)) ; debug_StressLoopDistribution = 0_pInt + allocate(debug_StateLoopDistribution(nState)) ; debug_StateLoopDistribution = 0_pInt + allocate(debug_StiffnessStateLoopDistribution(nState)) ; debug_StiffnessStateLoopDistribution = 0_pInt + allocate(debug_CrystalliteLoopDistribution(nCryst)) ; debug_CrystalliteLoopDistribution = 0_pInt +endsubroutine + !******************************************************************** ! reset debug distributions !******************************************************************** -SUBROUTINE debug_reset() +subroutine debug_reset() - use prec - implicit none + use prec + implicit none - debug_StressLoopDistribution = 0_pInt ! initialize debugging data - debug_StateLoopDistribution = 0_pInt - debug_StiffnessStateLoopDistribution = 0_pInt - debug_CrystalliteLoopDistribution = 0_pInt - debug_cumLpTicks = 0_pInt - debug_cumDotStateTicks = 0_pInt - debug_cumLpCalls = 0_pInt - debug_cumDotStateCalls = 0_pInt + debug_StressLoopDistribution = 0_pInt ! initialize debugging data + debug_StateLoopDistribution = 0_pInt + debug_StiffnessStateLoopDistribution = 0_pInt + debug_CrystalliteLoopDistribution = 0_pInt + debug_cumLpTicks = 0_pInt + debug_cumDotStateTicks = 0_pInt + debug_cumLpCalls = 0_pInt + debug_cumDotStateCalls = 0_pInt -END SUBROUTINE +endsubroutine !******************************************************************** ! write debug statements to standard out !******************************************************************** - SUBROUTINE debug_info() + subroutine debug_info() use prec + use numerics, only: nStress, & + nState, & + nCryst implicit none integer(pInt) i,integral integer(pLongInt) tickrate + write(6,*) write(6,*) 'DEBUG Info' @@ -112,6 +129,6 @@ END SUBROUTINE write(6,'(a15,i10,i10)') ' total',sum(debug_CrystalliteLoopDistribution),integral write(6,*) - END SUBROUTINE + endsubroutine END MODULE debug diff --git a/trunk/homogenization.f90 b/trunk/homogenization.f90 index 4c86ad96b..23217babc 100644 --- a/trunk/homogenization.f90 +++ b/trunk/homogenization.f90 @@ -155,7 +155,7 @@ subroutine homogenization_init() return -end subroutine +endsubroutine !******************************************************************** @@ -168,11 +168,17 @@ subroutine materialpoint_stressAndItsTangent(& dt & ! time increment ) - use prec, only: pInt,pReal, subStepMin,nHomog - use FEsolving, only: FEsolving_execElem, FEsolving_execIP - use mesh, only: mesh_element - use material, only: homogenization_Ngrains - use constitutive, only: constitutive_state0, constitutive_partionedState0, constitutive_state + use prec, only: pInt, & + pReal + use numerics, only: subStepMin, & + nHomog + use FEsolving, only: FEsolving_execElem, & + FEsolving_execIP + use mesh, only: mesh_element + use material, only: homogenization_Ngrains + use constitutive, only: constitutive_state0, & + constitutive_partionedState0, & + constitutive_state use crystallite implicit none @@ -331,7 +337,7 @@ subroutine materialpoint_stressAndItsTangent(& ! how to deal with stiffness? return -end subroutine +endsubroutine !******************************************************************** @@ -371,7 +377,7 @@ subroutine materialpoint_postResults(dt) enddo !$OMP END PARALLEL DO - end subroutine + endsubroutine !******************************************************************** @@ -400,7 +406,7 @@ subroutine homogenization_partitionDeformation(& homogenization_state(ip,el),ip,el) end select -end subroutine +endsubroutine !******************************************************************** @@ -431,7 +437,7 @@ function homogenization_updateState(& return -end function +endfunction !******************************************************************** @@ -459,7 +465,7 @@ subroutine homogenization_averageStressAndItsTangent(& return -end subroutine +endsubroutine !******************************************************************** @@ -489,6 +495,6 @@ function homogenization_postResults(& return -end function +endfunction END MODULE \ No newline at end of file diff --git a/trunk/homogenization_isostrain.f90 b/trunk/homogenization_isostrain.f90 index 38bd0be9a..fdae1ad97 100644 --- a/trunk/homogenization_isostrain.f90 +++ b/trunk/homogenization_isostrain.f90 @@ -109,7 +109,7 @@ subroutine homogenization_isostrain_init(& return -end subroutine +endsubroutine !********************************************************************* @@ -127,7 +127,7 @@ function homogenization_isostrain_stateInit(myInstance) return -end function +endfunction !******************************************************************** @@ -161,7 +161,7 @@ subroutine homogenization_isostrain_partitionDeformation(& return -end subroutine +endsubroutine !******************************************************************** @@ -195,7 +195,7 @@ function homogenization_isostrain_updateState(& return -end function +endfunction !******************************************************************** @@ -232,7 +232,7 @@ subroutine homogenization_isostrain_averageStressAndItsTangent(& return -end subroutine +endsubroutine !******************************************************************** @@ -270,6 +270,6 @@ pure function homogenization_isostrain_postResults(& return -end function +endfunction END MODULE diff --git a/trunk/lattice.f90 b/trunk/lattice.f90 index 6c4b48f20..17176db56 100644 --- a/trunk/lattice.f90 +++ b/trunk/lattice.f90 @@ -585,38 +585,38 @@ subroutine lattice_init() !************************************** !* Module initialization * !************************************** - use IO, only: IO_open_file,IO_countSections,IO_countTagInPart,IO_error - use material, only: material_configfile,material_partPhase - implicit none + use IO, only: IO_open_file,IO_countSections,IO_countTagInPart,IO_error + use material, only: material_configfile,material_partPhase + implicit none - integer(pInt), parameter :: fileunit = 200 - integer(pInt) i,Nsections + integer(pInt), parameter :: fileunit = 200 + integer(pInt) i,Nsections - if(.not. IO_open_file(fileunit,material_configFile)) call IO_error (100) ! corrupt config file - Nsections = IO_countSections(fileunit,material_partPhase) - lattice_Nstructure = 2_pInt + sum(IO_countTagInPart(fileunit,material_partPhase,'covera_ratio',Nsections)) ! fcc + bcc + all hex - close(fileunit) + if(.not. IO_open_file(fileunit,material_configFile)) call IO_error (100) ! corrupt config file + Nsections = IO_countSections(fileunit,material_partPhase) + lattice_Nstructure = 2_pInt + sum(IO_countTagInPart(fileunit,material_partPhase,'covera_ratio',Nsections)) ! fcc + bcc + all hex + close(fileunit) - allocate(lattice_Sslip(3,3,lattice_maxNslip,lattice_Nstructure)); lattice_Sslip = 0.0_pReal - allocate(lattice_Sslip_v(6,lattice_maxNslip,lattice_Nstructure)); lattice_Sslip_v = 0.0_pReal - allocate(lattice_sd(3,lattice_maxNslip,lattice_Nstructure)); lattice_sd = 0.0_pReal - allocate(lattice_st(3,lattice_maxNslip,lattice_Nstructure)); lattice_st = 0.0_pReal - allocate(lattice_sn(3,lattice_maxNslip,lattice_Nstructure)); lattice_sn = 0.0_pReal + allocate(lattice_Sslip(3,3,lattice_maxNslip,lattice_Nstructure)); lattice_Sslip = 0.0_pReal + allocate(lattice_Sslip_v(6,lattice_maxNslip,lattice_Nstructure)); lattice_Sslip_v = 0.0_pReal + allocate(lattice_sd(3,lattice_maxNslip,lattice_Nstructure)); lattice_sd = 0.0_pReal + allocate(lattice_st(3,lattice_maxNslip,lattice_Nstructure)); lattice_st = 0.0_pReal + allocate(lattice_sn(3,lattice_maxNslip,lattice_Nstructure)); lattice_sn = 0.0_pReal - allocate(lattice_Qtwin(3,3,lattice_maxNtwin,lattice_Nstructure)); lattice_Qtwin = 0.0_pReal - allocate(lattice_Stwin(3,3,lattice_maxNtwin,lattice_Nstructure)); lattice_Stwin = 0.0_pReal - allocate(lattice_Stwin_v(6,lattice_maxNtwin,lattice_Nstructure)); lattice_Stwin_v = 0.0_pReal - allocate(lattice_td(3,lattice_maxNtwin,lattice_Nstructure)); lattice_td = 0.0_pReal - allocate(lattice_tt(3,lattice_maxNtwin,lattice_Nstructure)); lattice_tt = 0.0_pReal - allocate(lattice_tn(3,lattice_maxNtwin,lattice_Nstructure)); lattice_tn = 0.0_pReal - allocate(lattice_shearTwin(lattice_maxNtwin,lattice_Nstructure)); lattice_shearTwin = 0.0_pReal + allocate(lattice_Qtwin(3,3,lattice_maxNtwin,lattice_Nstructure)); lattice_Qtwin = 0.0_pReal + allocate(lattice_Stwin(3,3,lattice_maxNtwin,lattice_Nstructure)); lattice_Stwin = 0.0_pReal + allocate(lattice_Stwin_v(6,lattice_maxNtwin,lattice_Nstructure)); lattice_Stwin_v = 0.0_pReal + allocate(lattice_td(3,lattice_maxNtwin,lattice_Nstructure)); lattice_td = 0.0_pReal + allocate(lattice_tt(3,lattice_maxNtwin,lattice_Nstructure)); lattice_tt = 0.0_pReal + allocate(lattice_tn(3,lattice_maxNtwin,lattice_Nstructure)); lattice_tn = 0.0_pReal + allocate(lattice_shearTwin(lattice_maxNtwin,lattice_Nstructure)); lattice_shearTwin = 0.0_pReal - allocate(lattice_interactionSlipSlip(lattice_maxNslip,lattice_maxNslip,lattice_Nstructure)); lattice_interactionSlipSlip = 0_pInt - allocate(lattice_interactionSlipTwin(lattice_maxNslip,lattice_maxNtwin,lattice_Nstructure)); lattice_interactionSlipTwin = 0_pInt - allocate(lattice_interactionTwinTwin(lattice_maxNtwin,lattice_maxNtwin,lattice_Nstructure)); lattice_interactionTwinTwin = 0_pInt - allocate(lattice_interactionTwinSlip(lattice_maxNtwin,lattice_maxNtwin,lattice_Nstructure)); lattice_interactionTwinSlip = 0_pInt - -end subroutine + allocate(lattice_interactionSlipSlip(lattice_maxNslip,lattice_maxNslip,lattice_Nstructure)); lattice_interactionSlipSlip = 0_pInt + allocate(lattice_interactionSlipTwin(lattice_maxNslip,lattice_maxNtwin,lattice_Nstructure)); lattice_interactionSlipTwin = 0_pInt + allocate(lattice_interactionTwinTwin(lattice_maxNtwin,lattice_maxNtwin,lattice_Nstructure)); lattice_interactionTwinTwin = 0_pInt + allocate(lattice_interactionTwinSlip(lattice_maxNtwin,lattice_maxNtwin,lattice_Nstructure)); lattice_interactionTwinSlip = 0_pInt + +endsubroutine function lattice_initializeStructure(struct,CoverA) @@ -758,7 +758,7 @@ function lattice_initializeStructure(struct,CoverA) lattice_initializeStructure = myStructure -end function +endfunction END MODULE diff --git a/trunk/material.f90 b/trunk/material.f90 index 96baf4ddb..1e9e4c42a 100644 --- a/trunk/material.f90 +++ b/trunk/material.f90 @@ -102,7 +102,7 @@ subroutine material_init() call material_populateGrains() -end subroutine +endsubroutine !********************************************************************* @@ -170,7 +170,7 @@ subroutine material_parseHomogenization(file,myPart) 100 homogenization_maxNgrains = maxval(homogenization_Ngrains) return - end subroutine + endsubroutine !********************************************************************* @@ -240,7 +240,7 @@ subroutine material_parseMicrostructure(file,myPart) 100 return - end subroutine + endsubroutine !********************************************************************* @@ -302,7 +302,7 @@ subroutine material_parsePhase(file,myPart) 100 return - end subroutine + endsubroutine !********************************************************************* @@ -418,7 +418,7 @@ subroutine material_parseTexture(file,myPart) 100 return - end subroutine + endsubroutine !********************************************************************* @@ -499,7 +499,7 @@ subroutine material_populateGrains() endif enddo NgrainsOfConstituent(t) = NgrainsOfConstituent(t) + sgn ! change that by one - end do + enddo ! ---------------------------------------------------------------------------- phaseOfGrain = 0_pInt orientationOfGrain = 0.0_pReal @@ -563,7 +563,7 @@ subroutine material_populateGrains() endif grain = grain + NgrainsOfConstituent(i) ! advance microstructure grain index - end do ! constituent + enddo ! constituent ! ---------------------------------------------------------------------------- do i=1,myNgrains-1 ! walk thru grains @@ -575,7 +575,7 @@ subroutine material_populateGrains() orientation = orientationOfGrain(:,t) orientationOfGrain(:,t) = orientationOfGrain(:,i) orientationOfGrain(:,i) = orientation - end do + enddo !calc fraction after weighing with volumePerGrain !exchange in MC steps to improve result... @@ -606,7 +606,7 @@ subroutine material_populateGrains() return - end subroutine + endsubroutine END MODULE diff --git a/trunk/mesh.f90 b/trunk/mesh.f90 index 64e5a3366..742e967f2 100644 --- a/trunk/mesh.f90 +++ b/trunk/mesh.f90 @@ -181,15 +181,16 @@ !*********************************************************** ! initialization !*********************************************************** - SUBROUTINE mesh_init () + subroutine mesh_init () use prec, only: pInt use IO, only: IO_error,IO_open_InputFile - use FEsolving, only: parallelExecution + use FEsolving, only: parallelExecution, FEsolving_execElem, FEsolving_execIP implicit none integer(pInt), parameter :: fileUnit = 222 + integer(pInt) e mesh_Nelems = 0_pInt mesh_NcpElems = 0_pInt @@ -228,14 +229,18 @@ call IO_error(100) ! cannot open input file endif - END SUBROUTINE + FEsolving_execElem = (/1,mesh_NcpElems/) + allocate(FEsolving_execIP(2,mesh_NcpElems)); FEsolving_execIP = 1_pInt + forall (e = 1:mesh_NcpElems) FEsolving_execIP(2,e) = FE_Nips(mesh_element(2,e)) + + endsubroutine !*********************************************************** ! mapping of FE element types to internal representation !*********************************************************** - FUNCTION FE_mapElemtype(what) + function FE_mapElemtype(what) implicit none @@ -259,9 +264,9 @@ FE_mapElemtype = 7 ! Three-dimensional Arbitrarily Distorted qudratic hexahedral case default FE_mapElemtype = 0 ! unknown element --> should raise an error upstream..! - end select + endselect - END FUNCTION + endfunction @@ -270,7 +275,7 @@ ! ! valid questions are 'elem', 'node' !*********************************************************** - FUNCTION mesh_FEasCP(what,id) + function mesh_FEasCP(what,id) use prec, only: pInt use IO, only: IO_lc @@ -289,7 +294,7 @@ lookupMap => mesh_mapFEtoCPnode case default return - end select + endselect lower = 1_pInt upper = size(lookupMap,2) @@ -313,17 +318,17 @@ else mesh_FEasCP = lookupMap(2,center) exit - end if - end do + endif + enddo return - END FUNCTION + endfunction !*********************************************************** ! find face-matching element of same type !!*********************************************************** - FUNCTION mesh_faceMatch(face,elem) + function mesh_faceMatch(face,elem) use prec, only: pInt implicit none @@ -344,7 +349,7 @@ minN = NsharedElems ! remember min # shared elems lonelyNode = faceNode ! remember most lonely node endif - end do + enddo candidate: do i=1,minN ! iterate over lonelyNode's shared elements mesh_faceMatch = mesh_sharedElem(1+i,nodeMap(lonelyNode)) ! present candidate elem if (mesh_faceMatch == elem) then ! my own element ? @@ -358,13 +363,13 @@ candidate: do i=1,minN ! iterate over lonelyNode's shared elements mesh_faceMatch = 0_pInt ! set to "no match" (so far) cycle candidate ! next candidate elem endif - end do + enddo exit ! surviving candidate - end do candidate + enddo candidate return - END FUNCTION + endfunction !******************************************************************** @@ -373,7 +378,7 @@ candidate: do i=1,minN ! iterate over lonelyNode's shared elements ! assign globals: ! FE_nodesAtIP, FE_ipNeighbor, FE_subNodeParent, FE_subNodeOnIPFace !******************************************************************** - SUBROUTINE mesh_get_FEdata () + subroutine mesh_get_FEdata () use prec, only: pInt implicit none @@ -1007,7 +1012,7 @@ candidate: do i=1,minN ! iterate over lonelyNode's shared elements return - END SUBROUTINE + endsubroutine !******************************************************************** ! get count of elements, nodes, and cp elements in mesh @@ -1016,7 +1021,7 @@ candidate: do i=1,minN ! iterate over lonelyNode's shared elements ! assign globals: ! _Nelems, _Nnodes, _NcpElems !******************************************************************** - SUBROUTINE mesh_get_meshDimensions (unit) + subroutine mesh_get_meshDimensions (unit) use prec, only: pInt use IO @@ -1046,19 +1051,19 @@ candidate: do i=1,minN ! iterate over lonelyNode's shared elements case('element') ! Count the number of encountered element sets mesh_NelemSets=mesh_NelemSets+1 mesh_maxNelemInSet = max(mesh_maxNelemInSet,IO_countContinousIntValues(unit)) - end select + endselect case('hypoelastic') do i=1,3+hypoelasticTableStyle ! Skip 3 or 4 lines read (unit,610,END=620) line - end do + enddo mesh_NcpElems = mesh_NcpElems + IO_countContinousIntValues(unit) - end select + endselect - end do + enddo 620 return - END SUBROUTINE + endsubroutine !!******************************************************************** @@ -1068,7 +1073,7 @@ candidate: do i=1,minN ! iterate over lonelyNode's shared elements ! assign globals: ! _maxNnodes, _maxNips, _maxNipNeighbors, _maxNsharedElems !******************************************************************** -SUBROUTINE mesh_get_nodeElemDimensions (unit) +subroutine mesh_get_nodeElemDimensions (unit) use prec, only: pInt use IO @@ -1107,25 +1112,25 @@ SUBROUTINE mesh_get_nodeElemDimensions (unit) n = mesh_FEasCP('node',IO_IntValue (line,pos,j+2)) if (all(node_seen /= n)) node_count(n) = node_count(n)+1 node_seen(j) = n - end do + enddo call IO_skipChunks(unit,FE_NoriginalNodes(t)-(pos(1)-2)) ! read on if FE_Nnodes exceeds node count present on current line - end if - end do + endif + enddo exit - end if - end do + endif + enddo 630 mesh_maxNsharedElems = maxval(node_count) return - END SUBROUTINE + endsubroutine !******************************************************************** ! Build element set mapping ! ! allocate globals: mesh_nameElemSet, mesh_mapElemSet !******************************************************************** - SUBROUTINE mesh_build_elemSetMapping (unit) + subroutine mesh_build_elemSetMapping (unit) use prec, only: pInt use IO @@ -1151,10 +1156,10 @@ SUBROUTINE mesh_get_nodeElemDimensions (unit) elem_set = elem_set+1 mesh_nameElemSet(elem_set) = IO_stringValue(line,pos,4) mesh_mapElemSet(:,elem_set) = IO_continousIntValues(unit,mesh_maxNelemInSet,mesh_nameElemSet,mesh_mapElemSet,mesh_NelemSets) - end if - end do + endif + enddo 640 return - END SUBROUTINE + endsubroutine !******************************************************************** @@ -1163,7 +1168,7 @@ SUBROUTINE mesh_get_nodeElemDimensions (unit) ! allocate globals: ! _mapFEtoCPnode !******************************************************************** - SUBROUTINE mesh_build_nodeMapping (unit) + subroutine mesh_build_nodeMapping (unit) use prec, only: pInt use math, only: qsort @@ -1191,15 +1196,15 @@ SUBROUTINE mesh_get_nodeElemDimensions (unit) read (unit,610,END=650) line mesh_mapFEtoCPnode(1,i) = IO_fixedIntValue (line,(/0,10/),1) mesh_mapFEtoCPnode(2,i) = i - end do + enddo exit - end if - end do + endif + enddo 650 call qsort(mesh_mapFEtoCPnode,1,size(mesh_mapFEtoCPnode,2)) return - END SUBROUTINE + endsubroutine !******************************************************************** @@ -1208,7 +1213,7 @@ SUBROUTINE mesh_get_nodeElemDimensions (unit) ! allocate globals: ! _mapFEtoCPelem !******************************************************************** - SUBROUTINE mesh_build_elemMapping (unit) + subroutine mesh_build_elemMapping (unit) use prec, only: pInt use math, only: qsort @@ -1234,20 +1239,20 @@ SUBROUTINE mesh_get_nodeElemDimensions (unit) if( IO_lc(IO_stringValue(line,pos,1)) == 'hypoelastic' ) then do i=1,3+hypoelasticTableStyle ! skip three (or four if new table style!) lines read (unit,610,END=660) line - end do + enddo contInts = IO_continousIntValues(unit,mesh_NcpElems,mesh_nameElemSet,mesh_mapElemSet,mesh_NelemSets) do i = 1,contInts(1) CP_elem = CP_elem+1 mesh_mapFEtoCPelem(1,CP_elem) = contInts(1+i) mesh_mapFEtoCPelem(2,CP_elem) = CP_elem enddo - end if - end do + endif + enddo 660 call qsort(mesh_mapFEtoCPelem,1,size(mesh_mapFEtoCPelem,2)) ! should be mesh_NcpElems return - END SUBROUTINE + endsubroutine !******************************************************************** @@ -1256,7 +1261,7 @@ SUBROUTINE mesh_get_nodeElemDimensions (unit) ! allocate globals: ! _node !******************************************************************** - SUBROUTINE mesh_build_nodes (unit) + subroutine mesh_build_nodes (unit) use prec, only: pInt use IO @@ -1283,15 +1288,15 @@ SUBROUTINE mesh_get_nodeElemDimensions (unit) m = mesh_FEasCP('node',IO_fixedIntValue (line,node_ends,1)) do j=1,3 mesh_node(j,m) = IO_fixedNoEFloatValue (line,node_ends,j+1) - end do - end do + enddo + enddo exit - end if - end do + endif + enddo 670 return - END SUBROUTINE + endsubroutine !******************************************************************** @@ -1300,7 +1305,7 @@ SUBROUTINE mesh_get_nodeElemDimensions (unit) ! allocate globals: ! _element !******************************************************************** - SUBROUTINE mesh_build_elements (unit) + subroutine mesh_build_elements (unit) use prec, only: pInt use IO @@ -1334,8 +1339,8 @@ SUBROUTINE mesh_get_nodeElemDimensions (unit) forall (j=1:FE_Nnodes(mesh_element(2,e))) & mesh_element(j+4,e) = IO_IntValue (line,pos,j+2) ! copy FE ids of nodes call IO_skipChunks(unit,FE_NoriginalNodes(mesh_element(2,e))-(pos(1)-2)) ! read on if FE_Nnodes exceeds node count present on current line - end if - end do + endif + enddo exit endif @@ -1378,7 +1383,7 @@ SUBROUTINE mesh_get_nodeElemDimensions (unit) 620 return - END SUBROUTINE + endsubroutine !******************************************************************** @@ -1387,7 +1392,7 @@ SUBROUTINE mesh_get_nodeElemDimensions (unit) ! allocate globals: ! _sharedElem !******************************************************************** - SUBROUTINE mesh_build_sharedElems (unit) + subroutine mesh_build_sharedElems (unit) use prec, only: pInt use IO @@ -1422,19 +1427,19 @@ SUBROUTINE mesh_get_nodeElemDimensions (unit) if (all(node_seen /= n)) then mesh_sharedElem(1,n) = mesh_sharedElem(1,n) + 1 mesh_sharedElem(1+mesh_sharedElem(1,n),n) = e - end if + endif node_seen(j) = n enddo call IO_skipChunks(unit,FE_NoriginalNodes(mesh_element(2,e))-(pos(1)-2)) ! read on if FE_Nnodes exceeds node count present on current line - end if - end do + endif + enddo exit - end if - end do + endif + enddo 620 return - END SUBROUTINE + endsubroutine !*********************************************************** @@ -1443,7 +1448,7 @@ SUBROUTINE mesh_get_nodeElemDimensions (unit) ! allocate globals ! _ipNeighborhood !*********************************************************** - SUBROUTINE mesh_build_ipNeighborhood() + subroutine mesh_build_ipNeighborhood() use prec, only: pInt implicit none @@ -1508,7 +1513,7 @@ SUBROUTINE mesh_get_nodeElemDimensions (unit) return - END SUBROUTINE + endsubroutine @@ -1518,7 +1523,7 @@ SUBROUTINE mesh_get_nodeElemDimensions (unit) ! allocate globals ! _subNodeCoord !*********************************************************** - SUBROUTINE mesh_build_subNodeCoords() + subroutine mesh_build_subNodeCoords() use prec, only: pInt,pReal implicit none @@ -1545,7 +1550,7 @@ SUBROUTINE mesh_get_nodeElemDimensions (unit) return - END SUBROUTINE + endsubroutine !*********************************************************** @@ -1554,7 +1559,7 @@ SUBROUTINE mesh_get_nodeElemDimensions (unit) ! allocate globals ! _ipVolume !*********************************************************** - SUBROUTINE mesh_build_ipVolumes() + subroutine mesh_build_ipVolumes() use prec, only: pInt use math, only: math_volTetrahedron @@ -1578,8 +1583,8 @@ SUBROUTINE mesh_get_nodeElemDimensions (unit) do n = 1,FE_NipFaceNodes ! loop over nodes on interface gravityNode(FE_subNodeOnIPFace(n,f,i,t)) = 1 gravityNodePos(:,FE_subNodeOnIPFace(n,f,i,t)) = mesh_subNodeCoord(:,FE_subNodeOnIPFace(n,f,i,t),e) - end do - end do + enddo + enddo do j = 1,mesh_maxNnodes+mesh_maxNsubNodes-1 ! walk through entire flagList except last if (gravityNode(j) > 0_pInt) then ! valid node index @@ -1588,10 +1593,10 @@ SUBROUTINE mesh_get_nodeElemDimensions (unit) gravityNode(j) = 0_pInt ! delete first instance gravityNodePos(:,j) = 0.0_pReal exit ! continue with next suspect - end if - end do - end if - end do + endif + enddo + endif + enddo centerOfGravity = sum(gravityNodePos,2)/count(gravityNode > 0) do f = 1,FE_NipNeighbors(t) ! loop over interfaces of IP and add tetrahedra which connect to CoG @@ -1602,13 +1607,13 @@ SUBROUTINE mesh_get_nodeElemDimensions (unit) nPos(:,1+mod(n+j-0,FE_NipFaceNodes)), & centerOfGravity) mesh_ipVolume(i,e) = mesh_ipVolume(i,e) + sum(volume) ! add contribution from this interface - end do + enddo mesh_ipVolume(i,e) = mesh_ipVolume(i,e) / FE_NipFaceNodes ! renormalize with interfaceNodeNum due to loop over them - end do - end do + enddo + enddo return - END SUBROUTINE + endsubroutine !*********************************************************** @@ -1617,7 +1622,7 @@ SUBROUTINE mesh_get_nodeElemDimensions (unit) ! allocate globals ! _ipArea, _ipAreaNormal !*********************************************************** - SUBROUTINE mesh_build_ipAreas() + subroutine mesh_build_ipAreas() use prec, only: pInt,pReal use math @@ -1652,7 +1657,7 @@ SUBROUTINE mesh_get_nodeElemDimensions (unit) enddo return - END SUBROUTINE + endsubroutine !*********************************************************** @@ -1660,7 +1665,7 @@ SUBROUTINE mesh_get_nodeElemDimensions (unit) ! to the output file ! !*********************************************************** - SUBROUTINE mesh_tell_statistics() + subroutine mesh_tell_statistics() use prec, only: pInt use math, only: math_range @@ -1707,9 +1712,9 @@ SUBROUTINE mesh_get_nodeElemDimensions (unit) write (6,"(i5,x,i5,x,e15.8)") e,i,mesh_IPvolume(i,e) do f = 1,FE_NipNeighbors(mesh_element(2,e)) write (6,"(i33,x,e15.8,x,3(f6.3,x))") f,mesh_ipArea(f,i,e),mesh_ipAreaNormal(:,f,i,e) - end do - end do - end do + enddo + enddo + enddo !write (6,*) !write (6,*) "Input Parser: SUBNODE COORDINATES" !write (6,*) @@ -1723,10 +1728,10 @@ SUBROUTINE mesh_get_nodeElemDimensions (unit) ! mesh_subNodeCoord(1,FE_subNodeOnIPFace(n,f,i,t),e),& ! mesh_subNodeCoord(2,FE_subNodeOnIPFace(n,f,i,t),e),& ! mesh_subNodeCoord(3,FE_subNodeOnIPFace(n,f,i,t),e) - ! end do - ! end do - ! end do - !end do + ! enddo + ! enddo + ! enddo + !enddo write (6,*) write (6,*) write (6,*) "Input Parser: STATISTICS" @@ -1757,7 +1762,7 @@ SUBROUTINE mesh_get_nodeElemDimensions (unit) return - END SUBROUTINE + endsubroutine END MODULE mesh diff --git a/trunk/mpie_cpfem_marc.f90 b/trunk/mpie_cpfem_marc.f90 index e75870d5a..ac7f9b353 100644 --- a/trunk/mpie_cpfem_marc.f90 +++ b/trunk/mpie_cpfem_marc.f90 @@ -16,6 +16,8 @@ ! - set statevariable 3 to index of microstructure ! - make sure the file "material.config" exists in the working ! directory +! - make sure the file "numerics.config" exists in the working +! directory ! - use nonsymmetric option for solver (e.g. direct ! profile or multifrontal sparse, the latter seems ! to be faster!) @@ -33,9 +35,10 @@ !******************************************************************** ! include "prec.f90" ! uses nothing else - include "debug.f90" ! uses prec include "math.f90" ! uses prec - include "IO.f90" ! uses prec, debug, math + include "IO.f90" ! uses prec, math + include "numerics.f90" ! uses prec, IO + include "debug.f90" ! uses prec, numerics include "FEsolving.f90" ! uses prec, IO include "mesh.f90" ! uses prec, math, IO, FEsolving include "material.f90" ! uses prec, math, IO, mesh @@ -44,10 +47,10 @@ include "constitutive_j2.f90" ! uses prec, math, IO, lattice, material, debug include "constitutive_dislobased.f90" ! uses prec, math, IO, lattice, material, debug include "constitutive.f90" ! uses prec, IO, math, lattice, mesh, debug - include "crystallite.f90" ! uses - include "homogenization_isostrain.f90" ! uses - include "homogenization.f90" ! uses - include "CPFEM.f90" ! uses prec, math, mesh, constitutive, FEsolving, debug, lattice, IO, crystallite + include "crystallite.f90" ! uses prec, math, IO, numerics + include "homogenization_isostrain.f90" ! uses prec, math, IO, + include "homogenization.f90" ! uses prec, math, IO, numerics + include "CPFEM.f90" ! uses prec, math, IO, numerics, debug, FEsolving, mesh, lattice, constitutive, crystallite !******************************************************************** @@ -123,11 +126,21 @@ subroutine hypela2(& ifu & ! set to 1 if stretch has been calculated ) - use prec, only: pReal,pInt, iJacoStiffness - use FEsolving - use CPFEM, only: CPFEM_general - use math, only: invnrmMandel - use debug + use prec, only: pReal, & + pInt + use FEsolving, only: cycleCounter, & + theInc, & + theCycle, & + theLovl, & + theTime, & + lastIncConverged, & + outdatedByNewInc, & + outdatedFFN1, & + symmetricSolver + use CPFEM, only: CPFEM_general + use math, only: invnrmMandel + use debug, only: debug_info, & + debug_reset implicit none ! ** Start of generated type statements ** @@ -147,7 +160,7 @@ subroutine hypela2(& include "concom%%MARCVERSION%%" ! concom is needed for inc, subinc, ncycle, lovl include "creeps%%MARCVERSION%%" ! creeps is needed for timinc (time increment) - integer(pInt) computationMode,i + integer(pInt) computationMode, i if (inc == 0) then cycleCounter = 4 @@ -186,7 +199,7 @@ subroutine hypela2(& theCycle = ncycle ! record current cycle count theLovl = lovl ! record current lovl - call CPFEM_general(computationMode,ffn,ffn1,t(1),timinc,n(1),nn,s,mod(cycleCounter-4,4_pInt*iJacoStiffness)==0,d,ngens) + call CPFEM_general(computationMode,ffn,ffn1,t(1),timinc,n(1),nn,s,d,ngens) ! Mandel: 11, 22, 33, SQRT(2)*12, SQRT(2)*23, SQRT(2)*13 ! Marc: 11, 22, 33, 12, 23, 13 diff --git a/trunk/numerics.config b/trunk/numerics.config new file mode 100644 index 000000000..3891604bb --- /dev/null +++ b/trunk/numerics.config @@ -0,0 +1,18 @@ + +relevantStrain 1.0e-7 +iJacoStiffness 1 +iJacoLpresiduum 1 +pert_Fg 1.0e-6 +nHomog 10 +nCryst 20 +nState 10 +nStress 40 +subStepMin 1.0e-3 +rTol_crystalliteState 1.0e-6 +rTol_crystalliteStress 1.0e-6 +aTol_crystalliteStress 1.0e-8 + +resToler 1.0e-4 +resAbsol 1.0e+2 +resBound 1.0e+1 +NRiterMax 24 \ No newline at end of file diff --git a/trunk/numerics.f90 b/trunk/numerics.f90 new file mode 100644 index 000000000..c60140e70 --- /dev/null +++ b/trunk/numerics.f90 @@ -0,0 +1,169 @@ +!############################################################## +MODULE numerics +!############################################################## + +use prec, only: pInt, pReal +implicit none + +character(len=64), parameter :: numerics_configFile = 'numerics.config' ! name of configuration file +integer(pInt) iJacoStiffness, & ! frequency of stiffness update + iJacoLpresiduum, & ! frequency of Jacobian update of residuum in Lp + nHomog, & ! homogenization loop limit + nCryst, & ! crystallite loop limit (only for debugging info, real loop limit is "subStepMin") + nState, & ! state loop limit + nStress, & ! stress loop limit + NRiterMax ! maximum number of GIA iteration +real(pReal) relevantStrain, & ! strain increment considered significant + pert_Fg, & ! strain perturbation for FEM Jacobi + subStepMin, & ! minimum (relative) size of sub-step allowed during cutback in crystallite + rTol_crystalliteState, & ! relative tolerance in crystallite state loop + rTol_crystalliteStress, & ! relative tolerance in crystallite stress loop + aTol_crystalliteStress, & ! absolute tolerance in crystallite stress loop + resToler, & ! relative tolerance of residual in GIA iteration + resAbsol, & ! absolute tolerance of residual in GIA iteration (corresponds to ~1 Pa) + resBound ! relative maximum value (upper bound) for GIA residual + + +CONTAINS + +!******************************************* +! initialization subroutine +!******************************************* +subroutine numerics_init() + + !*** variables and functions from other modules ***! + use prec, only: pInt, & + pReal + use IO, only: IO_error, & + IO_open_file, & + IO_isBlank, & + IO_stringPos, & + IO_stringValue, & + IO_lc, & + IO_floatValue, & + IO_intValue + + implicit none + + !*** input variables ***! + + !*** output variables ***! + + !*** local variables ***! + integer(pInt), parameter :: fileunit = 300 + integer(pInt), parameter :: maxNchunks = 2 + integer(pInt), dimension(1+2*maxNchunks) :: positions + character(len=64) tag + character(len=1024) line + + !*** global variables ***! + ! relevantStrain + ! iJacoStiffness + ! iJacoLpresiduum + ! pert_Fg + ! nHomog + ! nCryst + ! nState + ! nStress + ! subStepMin + ! rTol_crystalliteState + ! rTol_crystalliteStress + ! aTol_crystalliteStress + ! resToler + ! resAbsol + ! resBound + ! NRiterMax + + ! initialize all values to zero + relevantStrain = 0.0_pReal + iJacoStiffness = 0_pInt + iJacoLpresiduum = 0_pInt + pert_Fg = 0.0_pReal + nHomog = 0_pInt + nCryst = 0_pInt + nState = 0_pInt + nStress = 0_pInt + subStepMin = 0.0_pReal + rTol_crystalliteState = 0.0_pReal + rTol_crystalliteStress = 0.0_pReal + aTol_crystalliteStress = 0.0_pReal + resToler = 0.0_pReal + resAbsol = 0.0_pReal + resBound = 0.0_pReal + NRiterMax = 0_pInt + + ! try to open the config file and call error if corrupt + if(.not. IO_open_file(fileunit,numerics_configFile)) call IO_error (100) + + line = '' + ! read variables from config file + do + read(fileunit,'(a1024)',END=100) line + if (IO_isBlank(line)) cycle ! skip empty lines + positions = IO_stringPos(line,maxNchunks) + tag = IO_lc(IO_stringValue(line,positions,1)) ! extract key + select case(tag) + case ('relevantstrain') + relevantStrain = IO_floatValue(line,positions,2) + case ('ijacostiffness') + iJacoStiffness = IO_intValue(line,positions,2) + case ('ijacolpresiduum') + iJacoLpresiduum = IO_intValue(line,positions,2) + case ('pert_fg') + pert_Fg = IO_floatValue(line,positions,2) + case ('nhomog') + nHomog = IO_intValue(line,positions,2) + case ('ncryst') + nCryst = IO_intValue(line,positions,2) + case ('nstate') + nState = IO_intValue(line,positions,2) + case ('nstress') + nStress = IO_intValue(line,positions,2) + case ('substepmin') + subStepMin = IO_floatValue(line,positions,2) + case ('rtol_crystallitestate') + rTol_crystalliteState = IO_floatValue(line,positions,2) + case ('rtol_crystallitestress') + rTol_crystalliteStress = IO_floatValue(line,positions,2) + case ('atol_crystallitestress') + aTol_crystalliteStress = IO_floatValue(line,positions,2) + case ('restoler') + resToler = IO_floatValue(line,positions,2) + case ('resabsol') + resAbsol = IO_floatValue(line,positions,2) + case ('resbound') + resBound = IO_floatValue(line,positions,2) + case ('nritermax') + NRiterMax = IO_intValue(line,positions,2) + end select + enddo + + 100 write(6,*) + ! sanity check + if (relevantStrain <= 0.0_pReal) call IO_error(260) + if (iJacoStiffness < 1_pInt) call IO_error(261) + if (iJacoLpresiduum < 1_pInt) call IO_error(262) + if (pert_Fg <= 0.0_pReal) call IO_error(263) + if (nHomog < 1_pInt) call IO_error(264) + if (nCryst < 1_pInt) call IO_error(265) + if (nState < 1_pInt) call IO_error(266) + if (nStress < 1_pInt) call IO_error(267) + if (subStepMin <= 0.0_pReal) call IO_error(268) + if (rTol_crystalliteState <= 0.0_pReal) call IO_error(269) + if (rTol_crystalliteStress <= 0.0_pReal) call IO_error(270) + if (aTol_crystalliteStress <= 0.0_pReal) call IO_error(271) + if (resToler <= 0.0_pReal) call IO_error(272) + if (resAbsol <= 0.0_pReal) call IO_error(273) + if (resBound <= 0.0_pReal) call IO_error(274) + if (NRiterMax < 1_pInt) call IO_error(275) + + close(fileunit) + + write(6,*) + write(6,*) '<<<+- numerics init -+>>>' + write(6,*) '...done' + write(6,*) + +endsubroutine + +END MODULE numerics \ No newline at end of file diff --git a/trunk/prec.f90 b/trunk/prec.f90 index 6b812f86d..33d1f4d41 100644 --- a/trunk/prec.f90 +++ b/trunk/prec.f90 @@ -4,36 +4,14 @@ !############################################################## implicit none - + ! *** Precision of real and integer variables *** integer, parameter :: pReal = selected_real_kind(15,300) ! 15 significant digits, up to 1e+-300 integer, parameter :: pInt = selected_int_kind(9) ! up to +- 1e9 integer, parameter :: pLongInt = 8 ! should be 64bit - type :: p_vec real(pReal), dimension(:), pointer :: p end type p_vec - -! *** Strain increment considered significant *** - real(pReal), parameter :: relevantStrain = 1.0e-7_pReal - -! *** Numerical parameters *** - integer(pInt), parameter :: iJacoStiffness = 1_pInt ! frequency of stiffness update - integer(pInt), parameter :: iJacoLpresiduum = 6_pInt ! frequency of Jacobian update of residuum in Lp - real(pReal), parameter :: pert_Fg = 1.0e-6_pReal ! strain perturbation for FEM Jacobi - integer(pInt), parameter :: nHomog = 10_pInt ! homogenization loop limit - integer(pInt), parameter :: nCryst = 20_pInt ! crystallite loop limit (only for debugging info, real loop limit is "subStepMin") - integer(pInt), parameter :: nState = 10_pInt ! state loop limit - integer(pInt), parameter :: nStress = 40_pInt ! stress loop limit - real(pReal), parameter :: rTol_crystalliteState = 1.0e-5_pReal ! relative tolerance in crystallite state loop - real(pReal), parameter :: rTol_crystalliteStress = 1.0e-6_pReal ! relative tolerance in crystallite stress loop - real(pReal), parameter :: aTol_crystalliteStress = 1.0e-8_pReal ! absolute tolerance in crystallite stress loop - real(pReal), parameter :: subStepMin = 1.0e-3_pReal ! minimum (relative) size of sub-step allowed during cutback in crystallite -! - real(pReal), parameter :: resToler = 1.0e-4_pReal ! relative tolerance of residual in GIA iteration - real(pReal), parameter :: resAbsol = 1.0e+2_pReal ! absolute tolerance of residual in GIA iteration (corresponds to ~1 Pa) - real(pReal), parameter :: resBound = 1.0e+1_pReal ! relative maximum value (upper bound) for GIA residual - integer(pInt), parameter :: NRiterMax = 24_pInt ! maximum number of GIA iteration END MODULE prec