diff --git a/trunk/CPFEM_GIA8.f90 b/trunk/CPFEM_GIA8.f90 index f59a46f58..891e29024 100644 --- a/trunk/CPFEM_GIA8.f90 +++ b/trunk/CPFEM_GIA8.f90 @@ -136,7 +136,7 @@ ! integer(pInt) CPFEM_en, CPFEM_in, cp_en, CPFEM_ngens, i,j,k,l,m,n, e real(pReal), dimension (3,3) :: ffn,ffn1,Kirchhoff_bar - real(pReal), dimension (3,3,3,3) :: H_bar + real(pReal), dimension (3,3,3,3) :: H_bar, H_bar_sym real(pReal), dimension(CPFEM_ngens) :: CPFEM_stress real(pReal), dimension(CPFEM_ngens,CPFEM_ngens) :: CPFEM_jaco real(pReal) Temperature,CPFEM_dt,J_inverse @@ -208,6 +208,8 @@ math_I3(j,l)*CPFEM_ffn1_bar(i,m,CPFEM_in,cp_en)*CPFEM_PK1_bar(k,m,CPFEM_in,cp_en)) + & 0.5_pReal*(math_I3(i,k)*Kirchhoff_bar(j,l) + math_I3(j,l)*Kirchhoff_bar(i,k) + & math_I3(i,l)*Kirchhoff_bar(j,k) + math_I3(j,k)*Kirchhoff_bar(i,l)) + forall(i=1:3,j=1:3,k=1:3,l=1:3) & + H_bar_sym(i,j,k,l)= 0.25_pReal*(H_bar(i,j,k,l) +H_bar(j,i,k,l) +H_bar(i,j,l,k) +H_bar(j,i,l,k)) CPFEM_jaco_bar(1:CPFEM_ngens,1:CPFEM_ngens,CPFEM_in,cp_en) = math_Mandel3333to66(J_inverse*H_bar) ! case (3) ! collect and return odd result @@ -312,7 +314,8 @@ ! -------------- grain loop ----------------- do grain = 1,texture_Ngrains(mesh_element(4,cp_en)) call SingleCrystallite(msg,PK1(:,:,grain),dPdF(:,:,:,:,grain),& - CPFEM_results(5:4+constitutive_Nresults(grain,CPFEM_in,cp_en),grain,CPFEM_in,cp_en),& + CPFEM_results(CPFEM_Nresults+1:CPFEM_Nresults+constitutive_Nresults(grain,CPFEM_in,cp_en),& + grain,CPFEM_in,cp_en),& Fp1(:,:,grain),Fe1(:,:,grain),state1(:,grain),& ! output up to here dTime,cp_en,CPFEM_in,grain,.true.,& CPFEM_Temperature(CPFEM_in,cp_en),F1(:,:,grain),F0(:,:,grain),Fp0(:,:,grain),state0(:,grain)) @@ -500,7 +503,8 @@ call GIA_RelaxedDeformation(F1,F1_bar,rx) do grain = 1,8 call SingleCrystallite(msg,PK1(:,:,grain),dPdF(:,:,:,:,grain),& - CPFEM_results(5:4+constitutive_Nresults(grain,CPFEM_in,cp_en),grain,CPFEM_in,cp_en),& + CPFEM_results(CPFEM_Nresults+1:CPFEM_Nresults+constitutive_Nresults(grain,CPFEM_in,cp_en),& + grain,CPFEM_in,cp_en),& Fp1(:,:,grain),Fe1(:,:,grain),state1(:,grain),& ! output up to here dTime,cp_en,CPFEM_in,grain,.true.,& CPFEM_Temperature(CPFEM_in,cp_en),F1(:,:,grain),F0(:,:,grain),Fp0(:,:,grain),state0(:,grain)) diff --git a/trunk/CPFEM_Taylor.f90 b/trunk/CPFEM_Taylor.f90 index 7988ac6f2..b87fb30ca 100644 --- a/trunk/CPFEM_Taylor.f90 +++ b/trunk/CPFEM_Taylor.f90 @@ -124,7 +124,7 @@ ! integer(pInt) CPFEM_en, CPFEM_in, cp_en, CPFEM_ngens, i,j,k,l,m,n, e real(pReal), dimension (3,3) :: ffn,ffn1,Kirchhoff_bar - real(pReal), dimension (3,3,3,3) :: H_bar + real(pReal), dimension (3,3,3,3) :: H_bar, H_bar_sym real(pReal), dimension(CPFEM_ngens) :: CPFEM_stress real(pReal), dimension(CPFEM_ngens,CPFEM_ngens) :: CPFEM_jaco, odd_jaco real(pReal) Temperature,CPFEM_dt,J_inverse @@ -204,7 +204,8 @@ math_I3(j,l)*CPFEM_ffn1_bar(i,m,CPFEM_in,cp_en)*CPFEM_PK1_bar(k,m,CPFEM_in,cp_en) + & 0.5_pReal*(math_I3(i,k)*Kirchhoff_bar(j,l) + math_I3(j,l)*Kirchhoff_bar(i,k) + & math_I3(i,l)*Kirchhoff_bar(j,k) + math_I3(j,k)*Kirchhoff_bar(i,l)) -! Do we have to symmetrize H_bar? + forall(i=1:3,j=1:3,k=1:3,l=1:3) & + H_bar_sym(i,j,k,l)= 0.25_pReal*(H_bar(i,j,k,l) +H_bar(j,i,k,l) +H_bar(i,j,l,k) +H_bar(j,i,l,k)) CPFEM_jaco_bar(1:CPFEM_ngens,1:CPFEM_ngens,CPFEM_in,cp_en) = math_Mandel3333to66(J_inverse*H_bar) ! if (CPFEM_in==8 .and. cp_en==80) then ! do e=1,80 @@ -288,7 +289,8 @@ do grain = 1,texture_Ngrains(mesh_element(4,cp_en)) dPdF = dPdF_bar_old ! preguess consistent tangent of grain with avg call SingleCrystallite(msg,PK1,dPdF,& - CPFEM_results(5:4+constitutive_Nresults(grain,CPFEM_in,cp_en),grain,CPFEM_in,cp_en),& + CPFEM_results(CPFEM_Nresults+1:CPFEM_Nresults+constitutive_Nresults(grain,CPFEM_in,cp_en),& + grain,CPFEM_in,cp_en),& CPFEM_Lp(:,:,grain,CPFEM_in,cp_en),& CPFEM_Fp_new(:,:,grain,CPFEM_in,cp_en),Fe1,constitutive_state_new(:,grain,CPFEM_in,cp_en),& ! output up to here CPFEM_dt,cp_en,CPFEM_in,grain,updateJaco,& diff --git a/trunk/CPFEM_Taylor_sequential.f90 b/trunk/CPFEM_Taylor_sequential.f90 index 25ca68c3b..41d523bf3 100644 --- a/trunk/CPFEM_Taylor_sequential.f90 +++ b/trunk/CPFEM_Taylor_sequential.f90 @@ -125,7 +125,7 @@ ! integer(pInt) CPFEM_en, CPFEM_in, cp_en, CPFEM_ngens, i,j,k,l,m,n, e real(pReal), dimension (3,3) :: ffn,ffn1,Kirchhoff_bar - real(pReal), dimension (3,3,3,3) :: H_bar + real(pReal), dimension (3,3,3,3) :: H_bar, H_bar_sym real(pReal), dimension(CPFEM_ngens) :: CPFEM_stress real(pReal), dimension(CPFEM_ngens,CPFEM_ngens) :: CPFEM_jaco, odd_jaco real(pReal) Temperature,CPFEM_dt,J_inverse @@ -186,7 +186,9 @@ math_I3(j,l)*CPFEM_ffn1_bar(i,m,CPFEM_in,cp_en)*CPFEM_PK1_bar(k,m,CPFEM_in,cp_en) + & 0.5_pReal*(math_I3(i,k)*Kirchhoff_bar(j,l) + math_I3(j,l)*Kirchhoff_bar(i,k) + & math_I3(i,l)*Kirchhoff_bar(j,k) + math_I3(j,k)*Kirchhoff_bar(i,l)) - CPFEM_jaco_bar(1:CPFEM_ngens,1:CPFEM_ngens,CPFEM_in,cp_en) = math_Mandel3333to66(J_inverse*H_bar) + forall(i=1:3,j=1:3,k=1:3,l=1:3) & + H_bar_sym(i,j,k,l)= 0.25_pReal*(H_bar(i,j,k,l) +H_bar(j,i,k,l) +H_bar(i,j,l,k) +H_bar(j,i,l,k)) + CPFEM_jaco_bar(1:CPFEM_ngens,1:CPFEM_ngens,CPFEM_in,cp_en) = math_Mandel3333to66(J_inverse*H_bar) case (3) ! *** collect and return odd result *** CPFEM_Temperature(CPFEM_in,cp_en) = Temperature @@ -264,7 +266,8 @@ endif call SingleCrystallite(msg,PK1,dPdF,& - CPFEM_results(5:4+constitutive_Nresults(grain,CPFEM_in,cp_en),grain,CPFEM_in,cp_en),& + CPFEM_results(CPFEM_Nresults+1:CPFEM_Nresults+constitutive_Nresults(grain,CPFEM_in,cp_en),& + grain,CPFEM_in,cp_en),& CPFEM_Lp(:,:,grain,CPFEM_in,cp_en),& CPFEM_Fp_new(:,:,grain,CPFEM_in,cp_en),Fe1,constitutive_state_new(:,grain,CPFEM_in,cp_en),& ! output up to here CPFEM_dt,cp_en,CPFEM_in,grain,updateJaco,& diff --git a/trunk/FEsolving.f90 b/trunk/FEsolving.f90 index 5287b5456..d97b55ace 100644 --- a/trunk/FEsolving.f90 +++ b/trunk/FEsolving.f90 @@ -10,5 +10,42 @@ integer(pInt) theInc,theCycle,theLovl real(pReal) theTime logical :: lastIncConverged = .false.,outdatedByNewInc = .false., outdatedFFN1 = .false. + logical :: symmetricSolver = .false. + + CONTAINS + +!*********************************************************** +! determine wether a symmetric solver is used +!*********************************************************** + subroutine FE_get_solverSymmetry(unit) + + use prec, only: pInt + use IO + implicit none + + integer(pInt) unit + integer(pInt), dimension (133) :: pos + character*300 line + +610 FORMAT(A300) + + rewind(unit) + do + read (unit,610,END=630) line + pos = IO_stringPos(line,1) + if( IO_lc(IO_stringValue(line,pos,1)) == 'solver' ) then + read (unit,610,END=630) line ! Garbage line + pos = IO_stringPos(line,2) ! limit to 64 nodes max (plus ID, type) + if(IO_intValue(line,pos,2) /= 1_pInt) symmetricSolver = .true. +!$OMP CRITICAL (write2out) + write (6,*) + write (6,*) 'Symmetric solver detected. d-Matrix will be symmetrized!' +!$OMP END CRITICAL (write2out) + endif + enddo + +630 return + + end subroutine END MODULE FEsolving diff --git a/trunk/constitutive_pheno.f90 b/trunk/constitutive_pheno.f90 index 3e0d88227..b83028303 100644 --- a/trunk/constitutive_pheno.f90 +++ b/trunk/constitutive_pheno.f90 @@ -701,7 +701,7 @@ do e=1,mesh_NcpElems constitutive_MatVolFrac(g,i,e) = 1.0_pReal ! singular material (so far) constitutive_TexVolFrac(g,i,e) = texVolfrac(s)/multiplicity(texID)/Nsym(texID) constitutive_Nstatevars(g,i,e) = material_Nslip(matID) ! number of state variables (i.e. tau_c of each slip system) - constitutive_Nresults(g,i,e) = 0 ! number of constitutive results +! constitutive_Nresults(g,i,e) = 2*material_Nslip(matID) ! number of constitutive results (shears in this case) constitutive_EulerAngles(:,g,i,e) = Euler(:,s) ! store initial orientation forall (l=1:constitutive_Nstatevars(g,i,e)) ! initialize state variables constitutive_state_old(l,g,i,e) = material_s0_slip(matID) @@ -900,7 +900,7 @@ implicit none !* Definition of variables integer(pInt) ipc,ip,el integer(pInt) matID,i -real(pReal) dt,Temperature,tau_slip +real(pReal) dt,Temperature,tau_slip, active_rate real(pReal), dimension(6) :: Tstar_v real(pReal), dimension(constitutive_Nstatevars(ipc,ip,el)) :: state real(pReal), dimension(constitutive_Nresults(ipc,ip,el)) :: constitutive_post_results @@ -910,11 +910,18 @@ matID = constitutive_matID(ipc,ip,el) if(constitutive_Nresults(ipc,ip,el)==0) return +constitutive_post_results=0 do i=1,material_Nslip(matID) - constitutive_post_results(i) = state(i) +!do i=1,constitutive_Nresults(ipc,ip,el) +! constitutive_post_results(i) = state(i) tau_slip=dot_product(Tstar_v,lattice_Sslip_v(:,i,material_CrystalStructure(matID))) - constitutive_post_results(i+material_Nslip(matID)) = & - dt*material_gdot0_slip(matID)*(abs(tau_slip)/state(i))**material_n_slip(matID)*sign(1.0_pReal,tau_slip) +! constitutive_post_results(i+material_Nslip(matID)) = & + constitutive_post_results(i) = & + material_gdot0_slip(matID)*(abs(tau_slip)/state(i))**material_n_slip(matID)*sign(1.0_pReal,tau_slip) +enddo +active_rate = 0.1_pReal*MAXVAL(abs(constitutive_post_results)) +do i=1,material_Nslip(matID) + if(abs(constitutive_post_results(i)) > active_rate) constitutive_post_results(i+material_Nslip(matID))=1.0_pReal enddo return diff --git a/trunk/mesh.f90 b/trunk/mesh.f90 index b9cb604cd..95dc3985e 100644 --- a/trunk/mesh.f90 +++ b/trunk/mesh.f90 @@ -220,7 +220,8 @@ SUBROUTINE mesh_init () use prec, only: pInt - use IO, only: IO_error,IO_open_InputFile + use IO, only: IO_error,IO_open_InputFile + use FEsolving, only: FE_get_solverSymmetry implicit none integer(pInt), parameter :: fileUnit = 222 @@ -237,7 +238,8 @@ ! call to various subroutines to parse the stuff from the input file... - if (IO_open_inputFile(fileUnit)) then + if (IO_open_inputFile(fileUnit)) then + call FE_get_solverSymmetry(fileUnit) call mesh_get_meshDimensions(fileUnit) call mesh_build_nodeMapping(fileUnit) call mesh_build_elemMapping(fileUnit) diff --git a/trunk/mpie_cpfem_marc2005r3.f90 b/trunk/mpie_cpfem_marc2005r3.f90 index ac6b9c545..6514cd0b7 100644 --- a/trunk/mpie_cpfem_marc2005r3.f90 +++ b/trunk/mpie_cpfem_marc2005r3.f90 @@ -27,16 +27,16 @@ ! - creeps: timinc !******************************************************************** ! - include "prec.f90" - include "FEsolving.f90" - include "debug.f90" - include "math.f90" - include "IO.f90" - include "mesh.f90" - include "lattice.f90" - include "constitutive.f90" - include "crystallite.f90" - include "CPFEM.f90" + include "prec.f90" ! uses nothing else + include "debug.f90" ! uses prec + include "math.f90" ! uses prec + include "IO.f90" ! uses prec, debug, math + include "FEsolving.f90" ! uses prec, IO + include "mesh.f90" ! uses prec, IO, math, FEsolving + include "lattice.f90" ! uses prec, math + include "constitutive.f90" ! uses prec, IO, math, lattice, mesh, debug + include "crystallite.f90" ! uses prec, debug, constitutive, mesh, math, IO + include "CPFEM.f90" ! uses prec, math, mesh, constitutive, FEsolving, debug, lattice, IO, crystallite ! SUBROUTINE hypela2(d,g,e,de,s,t,dt,ngens,n,nn,kcus,matus,ndi,& @@ -236,7 +236,7 @@ ! Marc: 11, 22, 33, 12, 23, 13 forall(i=1:ngens) d(1:ngens,i) = invnrmMandel(i)*d(1:ngens,i)*invnrmMandel(1:ngens) s(1:ngens) = s(1:ngens)*invnrmMandel(1:ngens) - + if(symmetricSolver) d(1:ngens,1:ngens) = 0.5_pReal*(d(1:ngens,1:ngens)+transpose(d(1:ngens,1:ngens))) return END SUBROUTINE diff --git a/trunk/mpie_cpfem_marc2007r1.f90 b/trunk/mpie_cpfem_marc2007r1.f90 index c04b41314..89e4871d0 100644 --- a/trunk/mpie_cpfem_marc2007r1.f90 +++ b/trunk/mpie_cpfem_marc2007r1.f90 @@ -27,16 +27,16 @@ ! - creeps: timinc !******************************************************************** ! - include "prec.f90" - include "FEsolving.f90" - include "debug.f90" - include "math.f90" - include "IO.f90" - include "mesh.f90" - include "lattice.f90" - include "constitutive.f90" - include "crystallite.f90" - include "CPFEM.f90" + include "prec.f90" ! uses nothing else + include "debug.f90" ! uses prec + include "math.f90" ! uses prec + include "IO.f90" ! uses prec, debug, math + include "FEsolving.f90" ! uses prec, IO + include "mesh.f90" ! uses prec, IO, math, FEsolving + include "lattice.f90" ! uses prec, math + include "constitutive.f90" ! uses prec, IO, math, lattice, mesh, debug + include "crystallite.f90" ! uses prec, debug, constitutive, mesh, math, IO + include "CPFEM.f90" ! uses prec, math, mesh, constitutive, FEsolving, debug, lattice, IO, crystallite ! SUBROUTINE hypela2(d,g,e,de,s,t,dt,ngens,n,nn,kcus,matus,ndi,& @@ -213,7 +213,7 @@ ! Marc: 11, 22, 33, 12, 23, 13 forall(i=1:ngens) d(1:ngens,i) = invnrmMandel(i)*d(1:ngens,i)*invnrmMandel(1:ngens) s(1:ngens) = s(1:ngens)*invnrmMandel(1:ngens) - + if(symmetricSolver) d(1:ngens,1:ngens) = 0.5_pReal*(d(1:ngens,1:ngens)+transpose(d(1:ngens,1:ngens))) return END SUBROUTINE @@ -253,7 +253,7 @@ ! assign result variable v=CPFEM_results(mod(jpltcd-1_pInt, CPFEM_Nresults+constitutive_maxNresults)+1_pInt,& (jpltcd-1_pInt)/(CPFEM_Nresults+constitutive_maxNresults)+1_pInt,& - nn, mesh_FEasCP('elem', m)) + nn, mesh_FEasCP('elem', m)) return END SUBROUTINE ! diff --git a/trunk/mpie_cpfem_marc2007r1_sequential.f90 b/trunk/mpie_cpfem_marc2007r1_sequential.f90 index 31292607e..313de55df 100644 --- a/trunk/mpie_cpfem_marc2007r1_sequential.f90 +++ b/trunk/mpie_cpfem_marc2007r1_sequential.f90 @@ -27,16 +27,16 @@ ! - creeps: timinc !******************************************************************** ! - include "prec.f90" - include "FEsolving.f90" - include "debug.f90" - include "math.f90" - include "IO.f90" - include "mesh.f90" - include "lattice.f90" - include "constitutive.f90" - include "crystallite.f90" - include "CPFEM_sequential.f90" + include "prec.f90" ! uses nothing else + include "debug.f90" ! uses prec + include "math.f90" ! uses prec + include "IO.f90" ! uses prec, debug, math + include "FEsolving.f90" ! uses prec, IO + include "mesh.f90" ! uses prec, IO, math, FEsolving + include "lattice.f90" ! uses prec, math + include "constitutive.f90" ! uses prec, IO, math, lattice, mesh, debug + include "crystallite.f90" ! uses prec, debug, constitutive, mesh, math, IO + include "CPFEM_sequential.f90" ! uses prec, math, mesh, constitutive, FEsolving, debug, lattice, IO, crystallite ! SUBROUTINE hypela2(d,g,e,de,s,t,dt,ngens,n,nn,kcus,matus,ndi,& @@ -207,7 +207,7 @@ ! Marc: 11, 22, 33, 12, 23, 13 forall(i=1:ngens) d(1:ngens,i) = invnrmMandel(i)*d(1:ngens,i)*invnrmMandel(1:ngens) s(1:ngens) = s(1:ngens)*invnrmMandel(1:ngens) - + if(symmetricSolver) d(1:ngens,1:ngens) = 0.5_pReal*(d(1:ngens,1:ngens)+transpose(d(1:ngens,1:ngens))) return END SUBROUTINE