From 298cecbfecdcd7cead276d41399c6fa7064cb5ea Mon Sep 17 00:00:00 2001 From: Philip Eisenlohr Date: Tue, 3 Aug 2010 23:47:00 +0000 Subject: [PATCH] corrected syntax errors (long lines, line continuation by \) and logical mistake in mpie_cpfem_marc / abq_std which went unnoticed in ifort. Using SunStudio f90 surfaced those... --- code/FEsolving.f90 | 4 +- code/IO.f90 | 2 +- code/constitutive_dislotwin.f90 | 70 ++++++++++++++++++--------------- code/constitutive_nonlocal.f90 | 8 ++-- code/constitutive_titanmod.f90 | 17 ++++---- code/crystallite.f90 | 38 +++++++++--------- code/homogenization_RGC.f90 | 4 +- code/mesh.f90 | 20 +++++----- code/mpie_cpfem_abaqus_std.f | 4 +- code/mpie_cpfem_marc.f90 | 8 ++-- code/mpie_spectral.f90 | 9 +++-- 11 files changed, 98 insertions(+), 86 deletions(-) diff --git a/code/FEsolving.f90 b/code/FEsolving.f90 index c08efadc4..b739deaa3 100644 --- a/code/FEsolving.f90 +++ b/code/FEsolving.f90 @@ -43,10 +43,10 @@ rewind(fileunit) do read (fileunit,'(a1024)',END=100) line - positions = IO_stringPos(line,1) + positions(1:1+2*1) = IO_stringPos(line,1) if( IO_lc(IO_stringValue(line,positions,1)) == 'solver' ) then read (fileunit,'(a1024)',END=100) line ! Garbage line - positions = IO_stringPos(line,2) + positions(1:1+2*2) = IO_stringPos(line,2) symmetricSolver = (IO_intValue(line,positions,2) /= 1_pInt) exit endif diff --git a/code/IO.f90 b/code/IO.f90 index 85b4014ef..ac4ea4ef2 100644 --- a/code/IO.f90 +++ b/code/IO.f90 @@ -135,7 +135,7 @@ end function use mpie_interface implicit none - character(len=*), parameter :: pathSep = achar(47)//achar(92) ! /, \ + character(len=*), parameter :: pathSep = achar(47)//achar(92) ! forward and backward slash character(len=*) relPath integer(pInt) unit diff --git a/code/constitutive_dislotwin.f90 b/code/constitutive_dislotwin.f90 index 875e660cd..9f23af579 100644 --- a/code/constitutive_dislotwin.f90 +++ b/code/constitutive_dislotwin.f90 @@ -13,13 +13,13 @@ implicit none character(len=*), parameter :: constitutive_dislotwin_label = 'dislotwin' character(len=18), dimension(2), parameter:: constitutive_dislotwin_listBasicSlipStates = (/'rhoEdge ', & 'rhoEdgeDip'/) -character(len=18), dimension(1), parameter:: constitutive_dislotwin_listBasicTwinStates = (/'twinFraction'/) +character(len=18), dimension(1), parameter:: constitutive_dislotwin_listBasicTwinStates = (/'twinFraction'/) character(len=18), dimension(4), parameter:: constitutive_dislotwin_listDependentSlipStates =(/'invLambdaSlip ', & 'invLambdaSlipTwin', & 'meanFreePathSlip ', & 'tauSlipThreshold '/) character(len=18), dimension(4), parameter:: constitutive_dislotwin_listDependentTwinStates =(/'invLambdaTwin ', & - 'meanFreePathTwin', & + 'meanFreePathTwin', & 'tauTwinThreshold', & 'twinVolume '/) real(pReal), parameter :: kB = 1.38e-23_pReal ! Boltzmann constant in J/Kelvin @@ -453,7 +453,7 @@ do i = 1,maxNinstance 'resolved_stress_slip', & 'threshold_stress_slip', & 'edge_dipole_distance', & - 'stress_exponent' & + 'stress_exponent' & ) mySize = constitutive_dislotwin_totalNslip(i) case('twin_fraction', & @@ -817,7 +817,7 @@ forall (s = 1:ns) & state(g,ip,el)%p(5*ns+3*nt+s) = & constitutive_dislotwin_Gmod(myInstance)*constitutive_dislotwin_burgersPerSlipSystem(s,myInstance)*& sqrt(dot_product((state(g,ip,el)%p(1:ns)+state(g,ip,el)%p(ns+1:2*ns)),& - constitutive_dislotwin_interactionMatrixSlipSlip(1:ns,s,myInstance))) + constitutive_dislotwin_interactionMatrixSlipSlip(1:ns,s,myInstance))) !* threshold stress for growing twin forall (t = 1:nt) & @@ -920,7 +920,8 @@ do f = 1,lattice_maxNslipFamily ! loop over all constitutive_dislotwin_v0PerSlipSystem(f,myInstance) !* Shear rates due to slip - gdot_slip(j) = DotGamma0*exp(-BoltzmannRatio*(1-StressRatio_p)**constitutive_dislotwin_q(myInstance))*sign(1.0_pReal,tau_slip(j)) + gdot_slip(j) = DotGamma0*exp(-BoltzmannRatio*(1-StressRatio_p)**constitutive_dislotwin_q(myInstance))*& + sign(1.0_pReal,tau_slip(j)) !* Derivatives of shear rates dgdot_dtauslip(j) = & @@ -953,10 +954,10 @@ do f = 1,lattice_maxNtwinFamily ! loop over all !* Resolved shear stress on twin system tau_twin(j) = dot_product(Tstar_v,lattice_Stwin_v(:,index_myFamily+i,myStructure)) - !* Stress ratios + !* Stress ratios StressRatio_r = (state(g,ip,el)%p(6*ns+3*nt+j)/tau_twin(j))**constitutive_dislotwin_r(myInstance) - !* Shear rates and their derivatives due to twin + !* Shear rates and their derivatives due to twin if ( tau_twin(j) > 0.0_pReal ) then gdot_twin(j) = & (constitutive_dislotwin_MaxTwinFraction(myInstance)-sumf)*lattice_shearTwin(index_myFamily+i,myStructure)*& @@ -1008,10 +1009,13 @@ function constitutive_dislotwin_dotState(Tstar_v,Temperature,state,g,ip,el) use prec, only: pReal,pInt,p_vec use math, only: pi -use mesh, only: mesh_NcpElems,mesh_maxNips -use material, only: homogenization_maxNgrains,material_phase, phase_constitutionInstance -use lattice, only: lattice_Sslip,lattice_Sslip_v,lattice_Stwin,lattice_Stwin_v,lattice_maxNslipFamily,lattice_maxNtwinFamily, & - lattice_NslipSystem,lattice_NtwinSystem,lattice_shearTwin +use mesh, only: mesh_NcpElems, mesh_maxNips +use material, only: homogenization_maxNgrains, material_phase, phase_constitutionInstance +use lattice, only: lattice_Sslip,lattice_Sslip_v, & + lattice_Stwin,lattice_Stwin_v, & + lattice_maxNslipFamily,lattice_maxNtwinFamily, & + lattice_NslipSystem,lattice_NtwinSystem, & + lattice_shearTwin implicit none !* Input-Output variables @@ -1029,7 +1033,8 @@ real(pReal), dimension(constitutive_dislotwin_totalNslip(phase_constitutionInsta gdot_slip,tau_slip,DotRhoMultiplication,EdgeDipDistance,DotRhoEdgeEdgeAnnihilation,DotRhoEdgeDipAnnihilation,& ClimbVelocity,DotRhoEdgeDipClimb,DotRhoDipFormation -real(pReal), dimension(constitutive_dislotwin_totalNtwin(phase_constitutionInstance(material_phase(g,ip,el)))) :: gdot_twin,tau_twin +real(pReal), dimension(constitutive_dislotwin_totalNtwin(phase_constitutionInstance(material_phase(g,ip,el)))) :: & + gdot_twin,tau_twin !* Shortened notation myInstance = phase_constitutionInstance(material_phase(g,ip,el)) @@ -1064,7 +1069,8 @@ do f = 1,lattice_maxNslipFamily ! loop over all constitutive_dislotwin_v0PerSlipSystem(f,myInstance) !* Shear rates due to slip - gdot_slip(j) = DotGamma0*exp(-BoltzmannRatio*(1-StressRatio_p)**constitutive_dislotwin_q(myInstance))*sign(1.0_pReal,tau_slip(j)) + gdot_slip(j) = DotGamma0*exp(-BoltzmannRatio*(1-StressRatio_p)**constitutive_dislotwin_q(myInstance))*& + sign(1.0_pReal,tau_slip(j)) !* Multiplication DotRhoMultiplication(j) = abs(gdot_slip(j))/& @@ -1079,8 +1085,8 @@ do f = 1,lattice_maxNslipFamily ! loop over all EdgeDipDistance(j) = & (3.0_pReal*constitutive_dislotwin_Gmod(myInstance)*constitutive_dislotwin_burgersPerSlipSystem(f,myInstance))/& (16.0_pReal*pi*abs(tau_slip(j))) - if (EdgeDipDistance(j)>state(g,ip,el)%p(4*ns+2*nt+j)) EdgeDipDistance(j)=state(g,ip,el)%p(4*ns+2*nt+j) - if (EdgeDipDistance(j)state(g,ip,el)%p(4*ns+2*nt+j)) EdgeDipDistance(j)=state(g,ip,el)%p(4*ns+2*nt+j) + if (EdgeDipDistance(j) 0_pInt) then j = 0_pInt do f = 1,lattice_maxNtwinFamily ! loop over all slip families - index_myFamily = sum(lattice_NtwinSystem(1:f-1,myStructure)) ! at which index starts my family - do i = 1,constitutive_dislotwin_Ntwin(f,myInstance) ! process each (active) slip system in family + index_myFamily = sum(lattice_NtwinSystem(1:f-1,myStructure)) ! at which index starts my family + do i = 1,constitutive_dislotwin_Ntwin(f,myInstance) ! process each (active) slip system in family j = j + 1_pInt constitutive_dislotwin_postResults(c+j) = dot_product(Tstar_v,lattice_Stwin_v(:,index_myFamily+i,myStructure)) enddo; enddo @@ -1366,11 +1372,11 @@ do o = 1,phase_Noutput(material_phase(g,ip,el)) StressRatio_pminus1*(1-StressRatio_p)**(constitutive_dislotwin_q(myInstance)-1.0_pReal) !* Stress exponent - if (gdot_slip==0.0_pReal) then - constitutive_dislotwin_postResults(c+j) = 0.0_pReal - else - constitutive_dislotwin_postResults(c+j) = (tau/gdot_slip)*dgdot_dtauslip - endif + if (gdot_slip==0.0_pReal) then + constitutive_dislotwin_postResults(c+j) = 0.0_pReal + else + constitutive_dislotwin_postResults(c+j) = (tau/gdot_slip)*dgdot_dtauslip + endif enddo ; enddo c = c + ns diff --git a/code/constitutive_nonlocal.f90 b/code/constitutive_nonlocal.f90 index 8828d2a27..9c117eafc 100644 --- a/code/constitutive_nonlocal.f90 +++ b/code/constitutive_nonlocal.f90 @@ -583,11 +583,11 @@ do i = 1,maxNinstance constitutive_nonlocal_forestProjectionEdge(s1, s2, i) & = abs(math_mul3x3(lattice_sn(:, constitutive_nonlocal_slipSystemLattice(s1,i), myStructure), & - lattice_st(:, constitutive_nonlocal_slipSystemLattice(s2,i), myStructure))) ! forest projection of edge dislocations is the projection of (t = b x n) onto the slip normal of the respective slip plane + lattice_st(:, constitutive_nonlocal_slipSystemLattice(s2,i), myStructure))) ! forest projection of edge dislocations is the projection of (t = b x n) onto the slip normal of the respective slip plane constitutive_nonlocal_forestProjectionScrew(s1, s2, i) & = abs(math_mul3x3(lattice_sn(:, constitutive_nonlocal_slipSystemLattice(s1,i), myStructure), & - lattice_sd(:, constitutive_nonlocal_slipSystemLattice(s2,i), myStructure))) ! forest projection of screw dislocations is the projection of b onto the slip normal of the respective splip plane + lattice_sd(:, constitutive_nonlocal_slipSystemLattice(s2,i), myStructure))) ! forest projection of screw dislocations is the projection of b onto the slip normal of the respective splip plane !*** calculation of interaction matrices @@ -1128,7 +1128,7 @@ ns = constitutive_nonlocal_totalNslip(myInstance) forall (t = 1:8) & rhoSgl(:,t) = state(g,ip,el)%p((t-1)*ns+1:t*ns) -forall (s = 1:ns, t = 5:8, rhoSgl(s,t) * constitutive_nonlocal_v(s,t-4,g,ip,el) < 0.0_pReal) & ! contribution of used rho for changing sign of v +forall (s = 1:ns, t = 5:8, rhoSgl(s,t) * constitutive_nonlocal_v(s,t-4,g,ip,el) < 0.0_pReal) & ! contribution of used rho for changing sign of v rhoSgl(s,t-4) = rhoSgl(s,t-4) + abs(rhoSgl(s,t)) tauThreshold = state(g,ip,el)%p(11*ns+1:12*ns) @@ -1473,7 +1473,7 @@ do n = 1,FE_NipNeighbors(mesh_element(2,el)) neighboring_fluxdensity = fluxdensity ! ... then copy flux density to neighbor to ensure zero gradient in fluxdensity endif else ! if no neighbor existent... - if ( all(periodicSurfaceFlux(maxloc(abs(mesh_ipAreaNormal(:,n,ip,el))))) ) then ! ... and we want periodic fluxes at surface... + if ( all(periodicSurfaceFlux(maxloc(abs(mesh_ipAreaNormal(:,n,ip,el))))) ) then ! ... and we want periodic fluxes at surface... forall (t = 1:4) & ! ... then mirror fluxes neighboring_fluxdensity(:,t) = fluxdensity(:,t-1_pInt+2_pInt*mod(t,2_pInt)) else ! ... and we have a free surface... diff --git a/code/constitutive_titanmod.f90 b/code/constitutive_titanmod.f90 index b78ac2681..d452d9bf5 100644 --- a/code/constitutive_titanmod.f90 +++ b/code/constitutive_titanmod.f90 @@ -13,15 +13,15 @@ implicit none character(len=*), parameter :: constitutive_titanmod_label = 'titanmod' character(len=18), dimension(2), parameter:: constitutive_titanmod_listBasicSlipStates = (/'rho_edge ', & 'rho_screw'/) -character(len=18), dimension(1), parameter:: constitutive_titanmod_listBasicTwinStates = (/'twinFraction'/) +character(len=18), dimension(1), parameter:: constitutive_titanmod_listBasicTwinStates = (/'twinFraction'/) character(len=18), dimension(5), parameter:: constitutive_titanmod_listDependentSlipStates =(/'invLambdaSlipe', & - 'invLambdaSlips', & - 'etauSlipThreshold', & + 'invLambdaSlips', & + 'etauSlipThreshold', & 'stauSlipThreshold', & 'invLambdaSlipTwin'/) character(len=18), dimension(4), parameter:: constitutive_titanmod_listDependentTwinStates =(/'invLambdaTwin ', & - 'meanFreePathTwin', & + 'meanFreePathTwin', & 'tauTwinThreshold', & 'twinVolume '/) real(pReal), parameter :: kB = 1.38e-23_pReal ! Boltzmann constant in J/Kelvin @@ -1500,7 +1500,8 @@ do o = 1,phase_Noutput(material_phase(g,ip,el)) tau = dot_product(Tstar_v,lattice_Sslip_v(:,index_myFamily+i,myStructure)) !* Stress ratios StressRatio_p = (abs(tau)/state(g,ip,el)%p(5*ns+3*nt+j))**constitutive_titanmod_pe_PerSlipSystem(j,myInstance) - StressRatio_pminus1 = (abs(tau)/state(g,ip,el)%p(5*ns+3*nt+j))**(constitutive_titanmod_pe_PerSlipSystem(j,myInstance)-1.0_pReal) + StressRatio_pminus1 = (abs(tau)/state(g,ip,el)%p(5*ns+3*nt+j))**& + (constitutive_titanmod_pe_PerSlipSystem(j,myInstance)-1.0_pReal) !* Boltzmann ratio BoltzmannRatio = constitutive_titanmod_f0_PerSlipSystem(j,myInstance)/(kB*Temperature) !* Initial shear rates @@ -1511,15 +1512,15 @@ do o = 1,phase_Noutput(material_phase(g,ip,el)) !* Shear rates due to slip constitutive_titanmod_postResults(c+j) = & DotGamma0*exp(-BoltzmannRatio*(1-StressRatio_p)**constitutive_titanmod_qe_PerSlipSystem(j,myInstance))* & - sign(1.0_pReal,tau) + sign(1.0_pReal,tau) enddo ; enddo - + ! invLambdaSlipe', & ! 'invLambdaSlips', & ! 'etauSlipThreshold', & ! 'stauSlipThreshold', & ! 'invLambdaSlipTwin - + c = c + ns case ('edgesegment') constitutive_titanmod_postResults(c+1:c+ns) = state(g,ip,el)%p((4*ns+2*nt+1):(5*ns+2*nt)) diff --git a/code/crystallite.f90 b/code/crystallite.f90 index 04dc3b2b8..b584b2adc 100644 --- a/code/crystallite.f90 +++ b/code/crystallite.f90 @@ -712,7 +712,7 @@ subroutine crystallite_stressAndItsTangent(updateJaco) if (crystallite_todo(g,i,e)) then ! all undone crystallites call constitutive_collectDotState(crystallite_Tstar_v(:,g,i,e), crystallite_subTstar0_v(:,g,i,e), & crystallite_Fe, crystallite_Fp, crystallite_Temperature(g,i,e), & - crystallite_disorientation(:,:,g,i,e), crystallite_subdt(g,i,e), g, i, e) + crystallite_disorientation(:,:,g,i,e), crystallite_subdt(g,i,e), g, i, e) delta_dotState1 = constitutive_dotState(g,i,e)%p - constitutive_previousDotState(g,i,e)%p delta_dotState2 = constitutive_previousDotState(g,i,e)%p - constitutive_previousDotState2(g,i,e)%p dot_prod12 = dot_product(delta_dotState1, delta_dotState2) @@ -941,10 +941,10 @@ subroutine crystallite_stressAndItsTangent(updateJaco) do g = 1,myNgrains selectiveDebugger = (e == debug_e .and. i == debug_i .and. g == debug_g) if (crystallite_todo(g,i,e)) then - crystallite_todo(g,i,e) = crystallite_integrateStress(g,i,e) ! stress integration + crystallite_todo(g,i,e) = crystallite_integrateStress(g,i,e) ! stress integration if ( .not. crystallite_localConstitution(g,i,e) & - .and. .not. crystallite_todo(g,i,e)) & ! if broken non-local... - crystallite_todo = crystallite_todo .and. crystallite_localConstitution ! ...all non-locals skipped + .and. .not. crystallite_todo(g,i,e)) & ! if broken non-local... + crystallite_todo = crystallite_todo .and. crystallite_localConstitution ! ...all non-locals skipped endif enddo; enddo; enddo !$OMPEND PARALLEL DO @@ -954,7 +954,7 @@ subroutine crystallite_stressAndItsTangent(updateJaco) do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) do g = 1,myNgrains if (crystallite_todo(g,i,e)) & - constitutive_dotState(g,i,e)%p = 0.0_pReal ! zero out dotState + constitutive_dotState(g,i,e)%p = 0.0_pReal ! zero out dotState enddo; enddo; enddo crystallite_statedamper = 1.0_pReal @@ -991,21 +991,21 @@ subroutine crystallite_stressAndItsTangent(updateJaco) do g = 1,myNgrains selectiveDebugger = (e == debug_e .and. i == debug_i .and. g == debug_g) if (crystallite_todo(g,i,e)) then - crystallite_stateConverged(g,i,e) = crystallite_updateState(g,i,e) ! update state - crystallite_temperatureConverged(g,i,e) = crystallite_updateTemperature(g,i,e) ! update temperature + crystallite_stateConverged(g,i,e) = crystallite_updateState(g,i,e) ! update state + crystallite_temperatureConverged(g,i,e) = crystallite_updateTemperature(g,i,e) ! update temperature crystallite_converged(g,i,e) = crystallite_stateConverged(g,i,e) & .and. crystallite_temperatureConverged(g,i,e) if ( .not. crystallite_localConstitution(g,i,e) & - .and. .not. crystallite_todo(g,i,e)) & ! if updateState signals broken non-local... - crystallite_todo = crystallite_todo .and. crystallite_localConstitution ! ...all non-locals skipped + .and. .not. crystallite_todo(g,i,e)) & ! if updateState signals broken non-local... + crystallite_todo = crystallite_todo .and. crystallite_localConstitution ! ...all non-locals skipped endif enddo; enddo; enddo !$OMPEND PARALLEL DO - if (any(.not. crystallite_converged .and. .not. crystallite_localConstitution)) & ! any non-local not yet converged? - crystallite_converged = crystallite_converged .and. crystallite_localConstitution ! all non-local not converged + if (any(.not. crystallite_converged .and. .not. crystallite_localConstitution)) & ! any non-local not yet converged? + crystallite_converged = crystallite_converged .and. crystallite_localConstitution ! all non-local not converged - crystallite_todo = crystallite_todo .and. .not. crystallite_converged ! skip all converged + crystallite_todo = crystallite_todo .and. .not. crystallite_converged ! skip all converged enddo @@ -1013,10 +1013,10 @@ subroutine crystallite_stressAndItsTangent(updateJaco) myNgrains = homogenization_Ngrains(mesh_element(3,e)) do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) do g = 1,myNgrains - if (crystallite_converged(g,i,e)) then ! if stiffness calculation converged... - crystallite_dPdF(:,:,k,l,g,i,e) = (crystallite_P(:,:,g,i,e) - storedP(:,:,g,i,e))/pert_Fg ! ... use tangent dP_ij/dFg_kl - elseif (.not. storedConvergenceFlag(g,i,e)) then ! if crystallite didnŐt converge before... - crystallite_dPdF(:,:,:,:,g,i,e) = crystallite_fallbackdPdF(:,:,:,:,g,i,e) ! ... use (elastic) fallback + if (crystallite_converged(g,i,e)) then ! if stiffness calculation converged... + crystallite_dPdF(:,:,k,l,g,i,e) = (crystallite_P(:,:,g,i,e) - storedP(:,:,g,i,e))/pert_Fg ! ... use tangent dP_ij/dFg_kl + elseif (.not. storedConvergenceFlag(g,i,e)) then ! if crystallite didnŐt converge before... + crystallite_dPdF(:,:,:,:,g,i,e) = crystallite_fallbackdPdF(:,:,:,:,g,i,e) ! ... use (elastic) fallback endif enddo; enddo; enddo @@ -1436,7 +1436,8 @@ LpLoop: do if (error) then if (verboseDebugger .and. selectiveDebugger) then !$OMP CRITICAL (write2out) - write(6,'(a,i3,x,i2,x,i5,x,a,x,i3)') '::: integrateStress failed on dR/dLp inversion at ',g,i,e,' ; iteration ', NiterationStress + write(6,'(a,i3,x,i2,x,i5,x,a,x,i3)') '::: integrateStress failed on dR/dLp inversion at ',g,i,e, & + ' ; iteration ', NiterationStress write(6,*) write(6,'(a,/,9(9(e15.3,x)/))') 'dRdLp',dRdLp write(6,'(a,/,9(9(e15.3,x)/))') 'dLpdT_constitutive',dLpdT_constitutive @@ -1469,7 +1470,8 @@ LpLoop: do if (error) then if (verboseDebugger .and. selectiveDebugger) then !$OMP CRITICAL (write2out) - write(6,'(a,i3,x,i2,x,i5,x,a,x,i3)') '::: integrateStress failed on invFp_new inversion at ',g,i,e,' ; iteration ', NiterationStress + write(6,'(a,i3,x,i2,x,i5,x,a,x,i3)') '::: integrateStress failed on invFp_new inversion at ',g,i,e, & + ' ; iteration ', NiterationStress write(6,*) write(6,'(a11,3(i3,x),/,3(3(f12.7,x)/))') 'invFp_new at ',g,i,e,invFp_new !$OMPEND CRITICAL (write2out) diff --git a/code/homogenization_RGC.f90 b/code/homogenization_RGC.f90 index e5fa58a41..8d0e67416 100644 --- a/code/homogenization_RGC.f90 +++ b/code/homogenization_RGC.f90 @@ -476,7 +476,9 @@ function homogenization_RGC_updateState(& state%p(3*nIntFaceTot+8) = maxval(abs(drelax))/dt if (el == 1 .and. ip == 1) then write(6,'(x,a30,x,e14.8)')'Constitutive work: ',constitutiveWork - write(6,'(x,a30,3(x,e14.8))')'Magnitude mismatch: ',sum(NN(1,:))/dble(nGrain),sum(NN(2,:))/dble(nGrain),sum(NN(3,:))/dble(nGrain) + write(6,'(x,a30,3(x,e14.8))')'Magnitude mismatch: ',sum(NN(1,:))/dble(nGrain), & + sum(NN(2,:))/dble(nGrain), & + sum(NN(3,:))/dble(nGrain) write(6,'(x,a30,x,e14.8)')'Penalty energy: ',penaltyEnergy write(6,'(x,a30,x,e14.8)')'Volume discrepancy: ',volDiscrep write(6,*)'' diff --git a/code/mesh.f90 b/code/mesh.f90 index 88572043c..a74e62947 100644 --- a/code/mesh.f90 +++ b/code/mesh.f90 @@ -2006,7 +2006,7 @@ subroutine mesh_marc_count_cpSizes (unit) rewind(unit) do read (unit,610,END=630) line - pos = IO_stringPos(line,1) + pos = IO_stringPos(line,maxNchunks) if( IO_lc(IO_stringValue(line,pos,1)) == 'connectivity' ) then read (unit,610,END=630) line ! Garbage line do i=1,mesh_Nelems ! read all elements @@ -2060,7 +2060,7 @@ subroutine mesh_marc_count_cpSizes (unit) rewind(unit) do read (unit,610,END=620) line - pos = IO_stringPos(line,2) + pos = IO_stringPos(line,maxNchunks) if ( IO_lc(IO_stringValue(line,pos,1)) == '*part' ) inPart = .true. if ( IO_lc(IO_stringValue(line,pos,1)) == '*end' .and. & IO_lc(IO_stringValue(line,pos,2)) == 'part' ) inPart = .false. @@ -2348,7 +2348,7 @@ subroutine mesh_marc_count_cpSizes (unit) do while (e < mesh_NcpElems) read(unit,'(a1024)',END=110) line if (IO_isBlank(line)) cycle ! skip empty lines - pos = IO_stringPos(line,1) + pos(1:1+2*1) = IO_stringPos(line,1) e = e+1 ! valid element entry mesh_element ( 1,e) = e ! FE id @@ -2399,7 +2399,7 @@ subroutine mesh_marc_count_cpSizes (unit) rewind(unit) do read (unit,610,END=620) line - pos = IO_stringPos(line,1) + pos(1:1+2*1) = IO_stringPos(line,1) if( IO_lc(IO_stringValue(line,pos,1)) == 'connectivity' ) then read (unit,610,END=620) line ! Garbage line do i = 1,mesh_Nelems @@ -2421,16 +2421,16 @@ subroutine mesh_marc_count_cpSizes (unit) 620 rewind(unit) ! just in case "initial state" apears before "connectivity" read (unit,610,END=620) line do - pos = IO_stringPos(line,2) + pos(1:1+2*2) = IO_stringPos(line,2) if( (IO_lc(IO_stringValue(line,pos,1)) == 'initial') .and. & (IO_lc(IO_stringValue(line,pos,2)) == 'state') ) then if (initialcondTableStyle == 2) read (unit,610,END=620) line ! read extra line for new style read (unit,610,END=630) line ! read line with index of state var - pos = IO_stringPos(line,1) + pos(1:1+2*1) = IO_stringPos(line,1) sv = IO_IntValue(line,pos,1) ! figure state variable index if( (sv == 2).or.(sv == 3) ) then ! only state vars 2 and 3 of interest read (unit,610,END=620) line ! read line with value of state var - pos = IO_stringPos(line,1) + pos(1:1+2*1) = IO_stringPos(line,1) do while (scan(IO_stringValue(line,pos,1),'+-',back=.true.)>1) ! is noEfloat value? val = NINT(IO_fixedNoEFloatValue(line,(/0,20/),1)) ! state var's value mesh_maxValStateVar(sv-1) = max(val,mesh_maxValStateVar(sv-1)) ! remember max val of homogenization and microstructure index @@ -2445,7 +2445,7 @@ subroutine mesh_marc_count_cpSizes (unit) enddo if (initialcondTableStyle == 0) read (unit,610,END=620) line ! ignore IP range for old table style read (unit,610,END=630) line - pos = IO_stringPos(line,1) + pos(1:1+2*1) = IO_stringPos(line,1) enddo endif else @@ -2487,7 +2487,7 @@ subroutine mesh_marc_count_cpSizes (unit) rewind(unit) do read (unit,610,END=620) line - pos = IO_stringPos(line,2) + pos(1:1+2*2) = IO_stringPos(line,2) if ( IO_lc(IO_stringValue(line,pos,1)) == '*part' ) inPart = .true. if ( IO_lc(IO_stringValue(line,pos,1)) == '*end' .and. & IO_lc(IO_stringValue(line,pos,2)) == 'part' ) inPart = .false. @@ -2539,7 +2539,7 @@ subroutine mesh_marc_count_cpSizes (unit) if (i <= mesh_Nmaterials) then ! found one? elemSetName = mesh_mapMaterial(i) ! take corresponding elemSet read (unit,610,END=630) line ! read homogenization and microstructure - pos = IO_stringPos(line,2) + pos(1:1+2*2) = IO_stringPos(line,2) homog = NINT(IO_floatValue(line,pos,1)) micro = NINT(IO_floatValue(line,pos,2)) do i = 1,mesh_NelemSets ! look thru all elemSet definitions diff --git a/code/mpie_cpfem_abaqus_std.f b/code/mpie_cpfem_abaqus_std.f index 9ab89b09f..7361c284a 100644 --- a/code/mpie_cpfem_abaqus_std.f +++ b/code/mpie_cpfem_abaqus_std.f @@ -174,7 +174,7 @@ subroutine UMAT(STRESS,STATEV,DDSDDE,SSE,SPD,SCD,& calcMode(npt,cp_en) = .not. calcMode(npt,cp_en) ! ping pong (calc <--> collect) if ( calcMode(npt,cp_en) ) then ! now calc - if ( lastMode .ne. calcMode(npt,cp_en) ) then ! first after ping pong + if ( lastMode .neqv. calcMode(npt,cp_en) ) then ! first after ping pong call debug_reset() ! resets debugging outdatedFFN1 = .false. cycleCounter = cycleCounter + 1 @@ -186,7 +186,7 @@ subroutine UMAT(STRESS,STATEV,DDSDDE,SSE,SPD,SCD,& computationMode = 2 ! plain calc endif else ! now collect - if ( lastMode .ne. calcMode(npt,cp_en) ) call debug_info() ! first after ping pong reports debugging + if ( lastMode .neqv. calcMode(npt,cp_en) ) call debug_info() ! first after ping pong reports debugging if ( lastIncConverged ) then lastIncConverged = .false. computationMode = 4 ! collect and backup Jacobian after convergence diff --git a/code/mpie_cpfem_marc.f90 b/code/mpie_cpfem_marc.f90 index aefd619cc..fa311106a 100644 --- a/code/mpie_cpfem_marc.f90 +++ b/code/mpie_cpfem_marc.f90 @@ -57,7 +57,7 @@ end subroutine function getSolverWorkingDirectoryName() implicit none character(1024) getSolverWorkingDirectoryName, outName - character(len=*), parameter :: pathSep = achar(47)//achar(92) ! /, \ + character(len=*), parameter :: pathSep = achar(47)//achar(92) ! forward and backward slash getSolverWorkingDirectoryName='' outName='' @@ -71,7 +71,7 @@ function getSolverJobName() implicit none character(1024) getSolverJobName, outName - character(len=*), parameter :: pathSep = achar(47)//achar(92) ! /, \ + character(len=*), parameter :: pathSep = achar(47)//achar(92) ! forward and backward slash integer(pInt) extPos getSolverJobName='' @@ -262,7 +262,7 @@ subroutine hypela2(& calcMode(nn,cp_en) = .not. calcMode(nn,cp_en) ! ping pong (calc <--> collect) if ( calcMode(nn,cp_en) ) then ! now calc - if ( lastMode .ne. calcMode(nn,cp_en) ) then ! first after ping pong + if ( lastMode .neqv. calcMode(nn,cp_en) ) then ! first after ping pong call debug_reset() ! resets debugging outdatedFFN1 = .false. cycleCounter = cycleCounter + 1 @@ -274,7 +274,7 @@ subroutine hypela2(& computationMode = 2 ! plain calc endif else ! now collect - if ( lastMode /= calcMode(nn,cp_en) ) call debug_info() ! first after ping pong reports debugging + if ( lastMode .neqv. calcMode(nn,cp_en) ) call debug_info() ! first after ping pong reports debugging if ( lastIncConverged ) then lastIncConverged = .false. computationMode = 4 ! collect and backup Jacobian after convergence diff --git a/code/mpie_spectral.f90 b/code/mpie_spectral.f90 index b0675e2b7..80904b8f6 100644 --- a/code/mpie_spectral.f90 +++ b/code/mpie_spectral.f90 @@ -55,7 +55,7 @@ function getSolverWorkingDirectoryName() implicit none character(len=1024) cwd,outname,getSolverWorkingDirectoryName - character(len=*), parameter :: pathSep = achar(47)//achar(92) ! /, \ + character(len=*), parameter :: pathSep = achar(47)//achar(92) ! forward and backward slash call getarg(2,outname) ! path to loadFile @@ -83,7 +83,7 @@ function getSolverJobName() implicit none character(1024) getSolverJobName, outName, cwd - character(len=*), parameter :: pathSep = achar(47)//achar(92) ! /, \ + character(len=*), parameter :: pathSep = achar(47)//achar(92) ! forward and backward slash integer(pInt) posExt,posSep getSolverJobName = '' @@ -119,9 +119,10 @@ function getLoadcaseName() implicit none character(len=1024) getLoadcaseName, outName, cwd - character(len=*), parameter :: pathSep = achar(47)//achar(92) ! /, \ + character(len=*), parameter :: pathSep = achar(47)//achar(92) ! forward and backward slash integer(pInt) posExt,posSep - posExt = 0 !not sure if its needed + + posExt = 0 ! not sure if required call getarg(2,getLoadcaseName) posExt = scan(getLoadcaseName,'.',back=.true.)