OpenMP seems to work now

watch out for critical sections
write and matmul statements need to be marked accordingly!!!

trunk/CPFEM_GIA8.f90

@@ -78,6 +78,7 @@
  enddo
 !
 !    *** Output to MARC output file ***
+!$OMP CRITICAL (write2out)
  write(6,*)
  write(6,*) 'CPFEM Initialization'
  write(6,*)
@@ -98,6 +99,7 @@
  write(6,*) 'GIA_bNorm:           ', shape(GIA_bNorm)
  write(6,*)
  call flush(6)
+!$OMP END CRITICAL (write2out)
  return
 !
  END SUBROUTINE
@@ -156,21 +158,28 @@
  endif
 !
  cp_en = mesh_FEasCP('elem',CPFEM_en)
-  if (cp_en == 1 .and. CPFEM_in == 1) &
+  if (cp_en == 1 .and. CPFEM_in == 1) then
+!$OMP CRITICAL (write2out)
     write(6,'(a6,x,i4,x,a4,x,i4,x,a10,x,f8.4,x,a10,x,i2,x,a10,x,i2,x,a10,x,i2,x,a10,x,i2)') &
     'elem',cp_en,'IP',CPFEM_in,&
     'theTime',theTime,'theInc',theInc,'theCycle',theCycle,'theLovl',theLovl,&
     'mode',CPFEM_mode
+!$OMP END CRITICAL (write2out)
+  endif
 !
  select case (CPFEM_mode)
     case (2,1)     ! regular computation (with aging of results)
        if (.not. CPFEM_calc_done) then                ! puuh, me needs doing all the work...
+!$OMP CRITICAL (write2out)
            write (6,*) 'puuh me needs doing all the work', cp_en
+!$OMP END CRITICAL (write2out)
            if (CPFEM_mode == 1) then                  ! age results at start of new increment
              CPFEM_Fp_old            = CPFEM_Fp_new
              constitutive_state_old  = constitutive_state_new
              GIA_rVect_old           = GIA_rVect_new
+!$OMP CRITICAL (write2out)
              write (6,*) '#### aged results'
+!$OMP END CRITICAL (write2out)
            endif
            debug_cutbackDistribution = 0_pInt         ! initialize debugging data
            debug_InnerLoopDistribution = 0_pInt
@@ -186,7 +195,9 @@
            CPFEM_calc_done = .true.                   ! now calc is done
          endif    
 !       translate from P and dP/dF to CS and dCS/dE
+!$OMP CRITICAL (evilmatmul)
        Kirchhoff_bar = matmul(CPFEM_PK1_bar(:,:,CPFEM_in, cp_en),transpose(CPFEM_ffn1_bar(:,:,CPFEM_in, cp_en)))
+!$OMP END CRITICAL (evilmatmul)
        J_inverse  = 1.0_pReal/math_det3x3(CPFEM_ffn1_bar(:,:,CPFEM_in, cp_en))
        CPFEM_stress_bar(1:CPFEM_ngens,CPFEM_in,cp_en) = math_Mandel33to6(J_inverse*Kirchhoff_bar)
 !
@@ -306,7 +317,9 @@
                       dTime,cp_en,CPFEM_in,grain,.true.,&
                       CPFEM_Temperature(CPFEM_in,cp_en),F1(:,:,grain),F0(:,:,grain),Fp0(:,:,grain),state0(:,grain))
      if (msg /= 'ok') then                    ! solution not reached --> exit NRIteration
+!$OMP CRITICAL (write2out)
        write(6,*) 'GIA: grain loop failed to converge @ EL:',cp_en,' IP:',CPFEM_in
+!$OMP END CRITICAL (write2out)
        NRconvergent = .false.
        exit NRiteration
      endif
@@ -363,7 +376,10 @@
    enddo
    resNorm = sqrt(resNorm)
 !
-   if (debugger) write(6,'(x,a,i3,a,i3,a,i3,a,e10.4)')'EL:',cp_en,' IP:',CPFEM_in,' Iter:',NRiter,' RNorm:',resNorm
+   if (debugger) then
+!$OMP CRITICAL (write2out)
+    write(6,'(x,a,i3,a,i3,a,i3,a,e10.4)')'EL:',cp_en,' IP:',CPFEM_in,' Iter:',NRiter,' RNorm:',resNorm
+!$OMP END CRITICAL (write2out)
    if (NRiter == 1_pInt) resMax = resNorm
    if ((resNorm < resToler*resMax) .or. (resNorm < resAbsol)) then     ! resNorm < tolerance ===> convergent
      NRconvergent = .true.
@@ -403,7 +419,9 @@
      dvardres = 0.0_pReal
      call math_invert(36,dresdvar,dvardres,dummy,failed)
      if (failed) then
+!$OMP CRITICAL (write2out)
        write(6,*) 'GIA: failed to invert the Jacobian @ EL:',cp_en,' IP:',CPFEM_in
+!$OMP END CRITICAL (write2out)
        NRconvergent = .false.
        exit NRiteration
      endif
@@ -431,7 +449,9 @@
 !
 ! return to the general subroutine when convergence is not reached
  if (.not. NRconvergent) then
+!$OMP CRITICAL (write2out)
    write(6,'(x,a)') 'GIA: convergence is not reached @ EL:',cp_en,' IP:',CPFEM_in
+!$OMP END CRITICAL (write2out)
    call IO_error(600)
    return
  endif
@@ -450,11 +470,13 @@
 !   update results plotted in MENTAT
    call math_pDecomposition(Fe1(:,:,grain),U,R,error) ! polar decomposition
    if (error) then
+!$OMP CRITICAL (write2out)
      write(6,*) Fe1(:,:,grain)
      write(6,*) 'polar decomposition'
      write(6,*) 'Grain:             ',grain
      write(6,*) 'Integration point: ',CPFEM_in
      write(6,*) 'Element:           ',mesh_element(1,cp_en)
+!$OMP END CRITICAL (write2out)
      call IO_error(650)
      return
    endif
@@ -483,7 +505,9 @@
                       dTime,cp_en,CPFEM_in,grain,.true.,&
                       CPFEM_Temperature(CPFEM_in,cp_en),F1(:,:,grain),F0(:,:,grain),Fp0(:,:,grain),state0(:,grain))
          if (msg /= 'ok') then                    ! solution not reached --> exit NRIteration
+!$OMP CRITICAL (write2out)
            write(6,*) 'GIA: perturbation grain loop failed to converge within allowable step-size'
+!$OMP END CRITICAL (write2out)
            NRconvergent = .false.
            exit NRPerturbation
          endif
@@ -535,7 +559,11 @@
        enddo
        resNorm = sqrt(resNorm)
 !
-!       if (debugger) write(6,'(x,a,i3,a,i3,a,i3,a,i3,a,e10.4)')'EL = ',cp_en,':IP = ',CPFEM_in,':pert = ',3*(ip-1)+jp,':Iter = ',NRiter,':RNorm = ',resNorm
+!       if (debugger) then
+!!$OMP CRITICAL (write2out)
+!            write(6,'(x,a,i3,a,i3,a,i3,a,i3,a,e10.4)')'EL = ',cp_en,':IP = ',CPFEM_in,':pert = ',3*(ip-1)+jp,':Iter = ',NRiter,':RNorm = ',resNorm
+!!$OMP END CRITICAL (write2out)
+!       endif
        if (NRiter == 1_pInt) resMax = resNorm
        if ((resNorm < resToler*resMax) .or. (resNorm < resAbsol)) then     ! resNorm < tolerance ===> convergent
          NRconvergent = .true.
@@ -575,7 +603,9 @@
          dvardres = 0.0_pReal
          call math_invert(36,dresdvar,dvardres,dummy,failed)
          if (failed) then
+!$OMP CRITICAL (write2out)
            write(6,*) 'GIA: perturbation failed to invert the Jacobian'
+!$OMP END CRITICAL (write2out)
            NRconvergent = .false.
            exit NRPerturbation
          endif

trunk/CPFEM_Taylor.f90

@@ -171,10 +171,12 @@
            enddo
 !$OMP END PARALLEL DO
            call debug_info()                          ! output of debugging/performance statistics
-           CPFEM_calc_done = .true.                   ! now calc is done
+           CPFEM_calc_done = .true.                   ! now calc is done
-         endif    
+       endif    
 !       translate from P and dP/dF to CS and dCS/dE
+!$OMP CRITICAL (evilmatmul)
        Kirchhoff_bar = matmul(CPFEM_PK1_bar(:,:,CPFEM_in, cp_en),transpose(CPFEM_ffn1_bar(:,:,CPFEM_in, cp_en)))
+!$OMP END CRITICAL (evilmatmul)
        J_inverse  = 1.0_pReal/math_det3x3(CPFEM_ffn1_bar(:,:,CPFEM_in, cp_en))
        CPFEM_stress_bar(1:CPFEM_ngens,CPFEM_in,cp_en) = math_Mandel33to6(J_inverse*Kirchhoff_bar)
 !
@@ -187,7 +189,18 @@
                           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)
+       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
+!       do i=1,8
+!       write(6,*)
+!       write(6,*) e, i
+!       write(6,*) CPFEM_stress_bar(1:CPFEM_ngens,i,e)
+!       write(6,*)
+!       write(6,*) CPFEM_jaco_bar(1:CPFEM_ngens,1:CPFEM_ngens,i,e)
+!       enddo
+!   enddo
+! endif
 !
     case (3)    ! collect and return odd result
        CPFEM_Temperature(CPFEM_in,cp_en)  = Temperature
@@ -196,6 +209,15 @@
        CPFEM_stress_bar(1:CPFEM_ngens,CPFEM_in,cp_en) = CPFEM_odd_stress
        CPFEM_jaco_bar(1:CPFEM_ngens,1:CPFEM_ngens,CPFEM_in,cp_en) = CPFEM_odd_jacobian*math_identity2nd(CPFEM_ngens)
        CPFEM_calc_done = .false.
+! if (CPFEM_in==8 .and. cp_en==80) then
+!   do e=1,80
+!       do i=1,8
+!       write(6,*)
+!       write(6,*) e, i
+!       write(6,*) ffn1
+!       enddo
+!   enddo
+! endif
 
     case (4)    ! do nothing since we can recycle the former results (MARC specialty)
     case (5)    ! record consistent tangent at beginning of new increment
@@ -207,7 +229,7 @@
 ! return the local stress and the jacobian from storage
  CPFEM_stress(1:CPFEM_ngens) = CPFEM_stress_bar(1:CPFEM_ngens,CPFEM_in,cp_en)
  CPFEM_jaco(1:CPFEM_ngens,1:CPFEM_ngens) = CPFEM_jaco_bar(1:CPFEM_ngens,1:CPFEM_ngens,CPFEM_in,cp_en)
-! 
+!
  return
 !
  END SUBROUTINE
@@ -248,7 +270,7 @@
  endif
 
  do grain = 1,texture_Ngrains(mesh_element(4,cp_en))
-   dPdF = dPdF_bar_old                                                 ! preguess consistent tangent of grain with avg
+   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_Lp(:,:,grain,CPFEM_in,cp_en),&
@@ -259,7 +281,9 @@
                       CPFEM_Fp_old(:,:,grain,CPFEM_in,cp_en),constitutive_state_old(:,grain,CPFEM_in,cp_en))
 
    if (msg /= 'ok') then                                               ! solution not reached --> exit
+!$OMP CRITICAL (write2out)
        write(6,*) 'grain loop failed to converge @ EL:',cp_en,' IP:',CPFEM_in
+!$OMP END CRITICAL (write2out)
        call IO_error(600)
        return 
    endif

trunk/constitutive_dislo.f90

@@ -315,6 +315,7 @@ do while(.true.)
       section=section+1
    else
       if (section>0) then
+!$OMP CRITICAL (write2out)
          select case(tag)
 	     case ('lattice_structure')
               material_LatticeStructure(section)=IO_intValue(line,positions,2)
@@ -400,6 +401,7 @@ do while(.true.)
               material_c9(section)=IO_floatValue(line,positions,2)
 			  write(6,*) 'c9', material_c9(section)
          end select
+!$OMP END CRITICAL (write2out)
       endif
    endif
 enddo
@@ -697,7 +699,9 @@ do texID=1,texture_maxN
    multiplicity(texID) = max(1_pInt,texture_Ngrains(texID)/Ncomponents(texID)/Nsym(texID))
    if (mod(texture_Ngrains(texID),Ncomponents(texID)*Nsym(texID)) /= 0_pInt) then
       texture_Ngrains(texID) = multiplicity(texID)*Ncomponents(texID)*Nsym(texID)
+!$OMP CRITICAL (write2out)
       write (6,*) 'changed Ngrains to',texture_Ngrains(texID),' for texture',texID
+!$OMP END CRITICAL (write2out)
    endif
 enddo
 
@@ -930,8 +934,10 @@ matID = constitutive_matID(ipc,ip,el)
 startIdxTwin = material_Nslip(matID)
 
 !* Quantities derived from state - slip
+!$OMP CRITICAL (evilmatmul)
 constitutive_rho_f=matmul(constitutive_Pforest  (1:material_Nslip(matID),1:material_Nslip(matID),matID),state)
 constitutive_rho_p=matmul(constitutive_Pparallel(1:material_Nslip(matID),1:material_Nslip(matID),matID),state)	
+!$OMP END CRITICAL (evilmatmul)
 do i=1,material_Nslip(matID)
    constitutive_passing_stress(i) = material_tau0(matID)+material_c1(matID)*material_Gmod(matID)*material_bg(matID)*&
                                     sqrt(constitutive_rho_p(i))

trunk/constitutive_pheno.f90

@@ -432,7 +432,7 @@ fileunit=200
 !* First reading: number of materials and textures
 !-----------------------------
 !* determine material_maxN and texture_maxN from last respective parts
-if(IO_open_file(fileunit,filename)==.false.) goto 100
+if(.not. IO_open_file(fileunit,filename)) call IO_error (200) ! corrupt mattex file
 part = '_dummy_'
 do while (part/='')
    formerPart = part
@@ -561,7 +561,6 @@ enddo
 ! MISSING some consistency checks may be..?
 ! if ODFfile present then set NGauss NFiber =0
 return
-100 call IO_error(200) ! corrupt materials_textures file
 end subroutine
 
 
@@ -613,7 +612,9 @@ do texID=1,texture_maxN
    multiplicity(texID) = max(1_pInt,texture_Ngrains(texID)/Ncomponents(texID)/Nsym(texID))
    if (mod(texture_Ngrains(texID),Ncomponents(texID)*Nsym(texID)) /= 0_pInt) then
       texture_Ngrains(texID) = multiplicity(texID)*Ncomponents(texID)*Nsym(texID)
+!$OMP CRITICAL (write2out)
       write (6,*) 'changed Ngrains to',texture_Ngrains(texID),' for texture',texID
+!$OMP END CRITICAL (write2out)
    endif
 enddo
 
@@ -869,8 +870,10 @@ do i=1,constitutive_Nstatevars(ipc,ip,el)
 enddo
 
 !* Hardening for all systems
+!$OMP CRITICAL (evilmatmul)
 constitutive_dotState=matmul(constitutive_HardeningMatrix(1:material_Nslip(matID),1:material_Nslip(matID),&
                       matID),self_hardening)
+!$OMP END CRITICAL (evilmatmul)
 return
 end function
 
@@ -913,4 +916,4 @@ return
 
 end function
 
-END MODULE
+END MODULE

trunk/crystallite.f90

@@ -23,7 +23,6 @@ CONTAINS
      Fp_new,&     ! new plastic deformation gradient
      Fe_new,&     ! new "elastic" deformation gradient
      state_new,&  ! new state variable array
-!
      dt,&         ! time increment
      cp_en,&      ! element number
      ip,&         ! integration point number
@@ -64,15 +63,13 @@ CONTAINS
  nCutbacks = 0_pInt
 !
  Fg_aim = Fg_old                                                   ! make "new", "aim" a synonym for "old"
-!$OMP CRITICAL (fpnew)
  Fp_new = Fp_old
-!$OMP END CRITICAL (fpnew)
  call math_invert3x3(Fp_old,inv,det,error)
+!$OMP CRITICAL (evilmatmul)
  Fe_new = matmul(Fg_old,inv)
-!$OMP CRITICAL (statenew)
+!$OMP END CRITICAL (evilmatmul)
  state_bestguess = state_new                                       ! remember potentially available state guess
  state_new = state_old
-!$OMP END CRITICAL (statenew)
 !
  cuttedBack = .false.
  guessNew = .true.
@@ -92,8 +89,10 @@ CONTAINS
    if (guessNew) &
      state_new = state_bestguess                                   ! try best available guess for state
 
+!!$OMP CRITICAL (timeint)
    call TimeIntegration(msg,Lp,Fp_new,Fe_new,P,state_new,post_results,.true., &   ! def gradients and PK2 at end of time step
                         dt_aim,cp_en,ip,grain,Temperature,Fg_aim,Fp_current,state_current,0_pInt)
+!!$OMP END CRITICAL (timeint)
 !
    if (msg == 'ok') then
      cuttedBack = .false.                 ! no cut back required
@@ -101,14 +100,22 @@ CONTAINS
      subFrac = subFrac + subStep
      subStep = 1.0_pReal - subFrac        ! try one go
 
-	 if (debugger) write (6,*) '--------- one go -----------++##'
+	 if (debugger) then
+!$OMP CRITICAL (write2out)
+		write (6,*) '--------- one go -----------++##'
+!$OMP END CRITICAL (write2out)
+     endif
    else
      nCutbacks = nCutbacks + 1            ! record additional cutback
      cuttedBack = .true.                  ! encountered problems -->
      guessNew = .true.                    ! redo plastic Lp guess
      subStep = subStep / 2.0_pReal        ! cut time step in half
      state_bestguess = state_current      ! current state is then best guess
-	 if (debugger) write (6,*) '>>>>>>>>>>>>>>>>>>>> cutback <<<<<<<<<<<<<<<<<<<<<<'
+	 if (debugger) then
+!$OMP CRITICAL (write2out)
+		write (6,*) '>>>>>>>>>>>>>>>>>>>> cutback <<<<<<<<<<<<<<<<<<<<<<'
+!$OMP END CRITICAL (write2out)
+     endif
 
    endif
  enddo  ! potential substepping
@@ -132,8 +139,30 @@ CONTAINS
        Fg_pert(k,l) = Fg_pert(k,l) + pert_Fg  ! perturb single component
        Lp_pert    = Lp
        state_pert = state_new                 ! initial guess from end of time step
+!!$OMP CRITICAL (timeint)
        call TimeIntegration(msg,Lp_pert,Fp_pert,Fe_pert,P_pert,state_pert,post_results,.false., &   ! def gradients and PK2 at end of time step
                             dt_aim,cp_en,ip,grain,Temperature,Fg_pert,Fp_current,state_current,1_pInt)
+!!$OMP END CRITICAL (timeint)
+!!$OMP CRITICAL (write2out)
+!		if(cp_en==61 .and. ip==1) then
+!		write (6,*) k, l
+!		write (6,*) msg
+!		write (6,*) Lp_pert
+!		write (6,*) Fp_pert
+!		write (6,*) Fe_pert
+!		write (6,*) P_pert
+!		write (6,*) state_pert
+!		write (6,*) post_results
+!		write (6,*) dt_aim
+!		write (6,*) cp_en
+!		write (6,*) ip
+!		write (6,*) grain
+!		write (6,*) Temperature
+!		write (6,*) Fg_pert
+!		write (6,*) Fp_current
+!		write (6,*) state_current
+!		endif
+!!$OMP END CRITICAL (write2out)
        if (msg == 'ok') &
          dPdF(:,:,k,l) = (P_pert-P)/pert_Fg   ! constructing tangent dP_ij/dFg_kl only if valid forward difference
                                               ! otherwise leave component unchanged
@@ -211,12 +240,29 @@ CONTAINS
     return
  endif
 
- A = matmul(transpose(invFp_old), matmul(transpose(Fg_new),matmul(Fg_new,invFp_old)))
+!$OMP CRITICAL (evilmatmul)
+ A = matmul(transpose(invFp_old), matmul(transpose(Fg_new),matmul(Fg_new,invFp_old)))
+!$OMP END CRITICAL (evilmatmul)
+!!$OMP CRITICAL (write2out)
+!if(cp_en==61 .and. ip==1) then
+!    write(6,*)
+!    write(6,*) '*************************'
+!    write(6,*) iInner, iOuter
+!    write(6,*) cp_en, ip
+!    write(6,*) 'invFp_old'
+!    write(6,*) invFp_old
+!    write(6,*) 'Fg_new'
+!    write(6,*) Fg_new
+!    write(6,*) 'A'
+!    write(6,*) A
+!    write(6,*) '*************************'
+!    write(6,*)
+!    call flush(6)
+!endif
+!!$OMP END CRITICAL (write2out)
 !
  if (all(state == 0.0_pReal)) then
-!$OMP CRITICAL (statenew)
    state = state_old    ! former state guessed, if none specified
-!$OMP END CRITICAL (statenew)
  endif
  iOuter = 0_pInt                                   ! outer counter
 !
@@ -252,18 +298,24 @@ Inner: do              ! inner iteration: Lp
 !
          B = math_i3 - dt*Lpguess
          BT = transpose(B)
+!$OMP CRITICAL (evilmatmul)
          AB = matmul(A,B)
          BTA = matmul(BT,A)
          Tstar_v = 0.5_pReal*matmul(C_66,math_mandel33to6(matmul(BT,AB)-math_I3))
+!$OMP END CRITICAL (evilmatmul)
          Tstar = math_Mandel6to33(Tstar_v)
          p_hydro=(Tstar_v(1)+Tstar_v(2)+Tstar_v(3))/3.0_pReal
          forall(i=1:3) Tstar_v(i) = Tstar_v(i)-p_hydro                ! subtract hydrostatic pressure
+!!$OMP CRITICAL (calcLp)
          call constitutive_LpAndItsTangent(Lp,dLp, &
                                            Tstar_v,state,Temperature,grain,ip,cp_en)
+!!$OMP END CRITICAL (calcLp)
 !!$OMP CRITICAL (write2out)
+!if(cp_en==61 .and. ip==1) then
 !    write(6,*)
 !    write(6,*) '*************************'
-!    write(6,*) cp_en, ip
+!    write(6,*) iInner, iOuter
+!    write(6,*) cp_en, ip
 !    write(6,*) 'Tstar_v'
 !    write(6,*) Tstar_v
 !    write(6,*) 'A'
@@ -277,10 +329,24 @@ Inner: do              ! inner iteration: Lp
 !    write(6,*) 'C_66'
 !    write(6,*) C_66
 !    write(6,*) '*************************'
-!    write(6,*)
+!    write(6,*)
-!$OMP END CRITICAL (write2out)
+!	call flush(6)
+!endif
+!!$OMP END CRITICAL (write2out)
 !			    
          Rinner = Lpguess - Lp                                        ! update current residuum
+!!$OMP CRITICAL (write2out)
+!    write(6,*)
+!    write(6,*) '*************************'
+!    write(6,*) iInner, iOuter
+!    write(6,*) cp_en, ip
+!    write(6,*) 'Rinner'
+!    write(6,*) Rinner
+!    write(6,*) 'Lp'
+!    write(6,*) Lp
+!    write(6,*) 'Lpguess'
+!    write(6,*) Lpguess
+!!$OMP END CRITICAL (write2out)
 !
          if (.not.(any(Rinner/=Rinner)) .and. &                       ! exclude any NaN in residuum
              ( (maxval(abs(Rinner)) < abstol_Inner) .or. &            ! below abs tol .or.
@@ -312,7 +378,9 @@ Inner: do              ! inner iteration: Lp
              dTdLp(3*(i-1)+j,3*(k-1)+l) = dTdLp(3*(i-1)+j,3*(k-1)+l) + &
              C(i,j,l,n)*AB(k,n)+C(i,j,m,l)*BTA(m,k)
            dTdLp = -0.5_pReal*dt*dTdLp
+!$OMP CRITICAL (evilmatmul)
            dRdLp = eye2 - matmul(dLp,dTdLp)                           ! calc dR/dLp
+!$OMP END CRITICAL (evilmatmul)
            invdRdLp = 0.0_pReal
            call math_invert(9,dRdLp,invdRdLp,dummy,failed)            ! invert dR/dLp --> dLp/dR
            if (failed) then
@@ -345,12 +413,12 @@ Inner: do              ! inner iteration: Lp
 !$OMP CRITICAL (in)
        debug_InnerLoopDistribution(iInner) = debug_InnerLoopDistribution(iInner)+1
 !$OMP END CRITICAL (in)
+!!$OMP CRITICAL (stateupdate)
        ROuter = state - state_old - &
                 dt*constitutive_dotState(Tstar_v,state,Temperature,&
                                          grain,ip,cp_en)          ! residuum from evolution of microstructure
-!$OMP CRITICAL (statenew)
+!!$OMP END CRITICAL (stateupdate)
        state = state - ROuter                                           ! update of microstructure
-!$OMP END CRITICAL (statenew)
        if (maxval(abs(Router/state),state /= 0.0_pReal) < reltol_Outer) exit Outer
      enddo Outer
 !
@@ -358,29 +426,36 @@ Inner: do              ! inner iteration: Lp
  debug_OuterLoopDistribution(iOuter) = debug_OuterLoopDistribution(iOuter)+1
 !$OMP END CRITICAL (out)
 
+!$OMP CRITICAL (evilmatmul)
  invFp_new = matmul(invFp_old,B)
-!$OMP CRITICAL (fpnew)
+!$OMP END CRITICAL (evilmatmul)
  call math_invert3x3(invFp_new,Fp_new,det,failed)
-!$OMP END CRITICAL (fpnew)
  if (failed) then
    msg = 'inversion Fp_new^-1'
    return
  endif
 !
  if (wantsConstitutiveResults) then     ! get the post_results upon request
-!$OMP CRITICAL (res)
    results = 0.0_pReal
    results = constitutive_post_results(Tstar_v,state,Temperature,dt,grain,ip,cp_en)
-!$OMP END CRITICAL (res)
  endif
 !
-!$OMP CRITICAL (fpnew)
  Fp_new = Fp_new*det**(1.0_pReal/3.0_pReal)     ! regularize Fp by det = det(InvFp_new) !!
- Fe_new = matmul(Fg_new,invFp_new)              ! calc resulting Fe
-!$OMP END CRITICAL (fpnew)
  forall (i=1:3) Tstar_v(i) = Tstar_v(i)+p_hydro ! add hydrostatic component back
+!$OMP CRITICAL (evilmatmul)
+ Fe_new = matmul(Fg_new,invFp_new)              ! calc resulting Fe
  P = matmul(Fe_new,matmul(Tstar,transpose(invFp_new)))    ! first PK stress
+!$OMP END CRITICAL (evilmatmul)
 !
+!!$OMP CRITICAL (write2out)
+!    write(6,*)
+!    write(6,*) '*************************'
+!    write(6,*) cp_en, ip
+!    write(6,*) iInner, iOuter
+!    write(6,*) 'P'
+!    write(6,*) P
+!    write(6,*) '*************************'
+!!$OMP END CRITICAL (write2out)
  return
 !
  END SUBROUTINE

trunk/math.f90

@@ -823,7 +823,9 @@
  real(pReal), dimension(3,3) :: r
  real(pReal) math_disorient, tr
 
+!$OMP CRITICAL (evilmatmul)
  r = matmul(math_EulerToR(EulerB),transpose(math_EulerToR(EulerA)))
+!$OMP END CRITICAL (evilmatmul)
  tr = (r(1,1)+r(2,2)+r(3,3)-1.0_pReal)*0.4999999_pReal
  math_disorient = abs(0.5_pReal*pi-asin(tr))
  return
@@ -887,8 +889,10 @@ endif
    disturb(2) = sign(1.0_pReal,rnd(2))*acos(cosScatter+(1.0_pReal-cosScatter)*rnd(4)) ! Phi
    disturb(3) = scatter * rnd(2)                                                      ! phi2
    if (rnd(5) <= exp(-1.0_pReal*(math_disorient(origin,disturb)/scatter)**2)) exit   
- enddo
+ enddo
+!$OMP CRITICAL (evilmatmul)
  math_sampleGaussOri = math_RtoEuler(matmul(math_EulerToR(disturb),math_EulerToR(center)))
+!$OMP END CRITICAL (evilmatmul)
  return
  
  END FUNCTION
@@ -962,7 +966,9 @@ endif
  pRot = math_RodrigtoR(axis,angle)
 
 ! ---# apply the three rotations #---
- math_sampleFiberOri = math_RtoEuler(matmul(pRot,matmul(fRot,oRot))) 
+!$OMP CRITICAL (evilmatmul)
+math_sampleFiberOri = math_RtoEuler(matmul(pRot,matmul(fRot,oRot))) 
+!$OMP END CRITICAL (evilmatmul)
  return
 
  END FUNCTION
@@ -1029,11 +1035,17 @@ endif
  real(pReal) FE(3,3),R(3,3),U(3,3),CE(3,3),EW1,EW2,EW3,EB1(3,3),EB2(3,3),EB3(3,3),UI(3,3),det
 
  error = .false.
+!$OMP CRITICAL (evilmatmul)
  ce=matmul(transpose(fe),fe)
+!$OMP END CRITICAL (evilmatmul)
  CALL math_spectral1(CE,EW1,EW2,EW3,EB1,EB2,EB3)
  U=DSQRT(EW1)*EB1+DSQRT(EW2)*EB2+DSQRT(EW3)*EB3
  call math_invert3x3(U,UI,det,error)
- if (.not. error) R = matmul(fe,ui)
+ if (.not. error) then
+!$OMP CRITICAL (evilmatmul)
+ R = matmul(fe,ui)
+!$OMP END CRITICAL (evilmatmul)
+ endif
  return 
  
  END SUBROUTINE
@@ -1090,7 +1102,9 @@ endif
   D3=1.0_pReal/(EW3-EW1)/(EW3-EW2)
   M1=M-EW1*math_I3
   M2=M-EW2*math_I3
+!$OMP CRITICAL (evilmatmul)
   EB3=MATMUL(M1,M2)*D3
+!$OMP END CRITICAL (evilmatmul)
   EB1=math_I3-EB3
 !  both EB2 and EW2 are set to zero so that they do not
 !  contribute to U in PDECOMPOSITION
@@ -1100,7 +1114,9 @@ endif
   D1=1.0_pReal/(EW1-EW2)/(EW1-EW3)
   M2=M-math_I3*EW2
   M3=M-math_I3*EW3
+!$OMP CRITICAL (evilmatmul)
   EB1=MATMUL(M2,M3)*D1
+!$OMP END CRITICAL (evilmatmul)
   EB2=math_I3-EB1
 !  both EB3 and EW3 are set to zero so that they do not
 !  contribute to U in PDECOMPOSITION
@@ -1110,7 +1126,9 @@ endif
   D2=1.0_pReal/(EW2-EW1)/(EW2-EW3) 
   M1=M-math_I3*EW1
   M3=M-math_I3*EW3
+!$OMP CRITICAL (evilmatmul)
   EB2=MATMUL(M1,M3)*D2
+!$OMP END CRITICAL (evilmatmul)
   EB1=math_I3-EB2
 !  both EB3 and EW3 are set to zero so that they do not
 !  contribute to U in PDECOMPOSITION
@@ -1123,9 +1141,11 @@ endif
   M1=M-EW1*math_I3
   M2=M-EW2*math_I3
   M3=M-EW3*math_I3
+!$OMP CRITICAL (evilmatmul)
   EB1=MATMUL(M2,M3)*D1
   EB2=MATMUL(M1,M3)*D2
   EB3=MATMUL(M1,M2)*D3
+!$OMP END CRITICAL (evilmatmul)
    END IF
  END IF
  RETURN