removed hard-coded tolerances from mpie_spectral and put them to numerics/IO

removed storage of full cauchy stres field from mpie_spectral.f90, only average is stored now
added cauchy stress and von mises equivalent calculation to spectral post.

code/IO.f90

@@ -1146,7 +1146,7 @@ endfunction
  case (47)
    msg = 'mask consistency violated in spectral loadcase'
  case (48)
-   msg = 'Non-positive relative tolerance for defGrad correction'
+   msg = 'Non-positive relative parameter for spectral method'
  case (50)
    msg = 'Error writing constitutive output description'
  case (100)

code/mpie_spectral.f90

@@ -29,7 +29,7 @@ program mpie_spectral
  use IO
  use math
  use CPFEM, only: CPFEM_general, CPFEM_initAll 
- use numerics, only: relevantStrain, rTol_crystalliteStress, mpieNumThreadsInt        !ToDo: change to really needed variables
+ use numerics, only: mpieNumThreadsInt, err_div_tol, err_stress_tol, err_stress_tolrel, err_defgrad_tol, itmax       
  use homogenization, only: materialpoint_sizeResults, materialpoint_results
 !$ use OMP_LIB                                                          ! the openMP function library
 
@@ -72,7 +72,7 @@ program mpie_spectral
  real(pReal), dimension(6) ::                           cstress                                ! cauchy stress in Mandel notation
  real(pReal), dimension(6,6) ::                         dsde, c066, s066                         
  real(pReal), dimension(:,:,:), allocatable ::          ddefgrad                  
- real(pReal), dimension(:,:,:,:,:), allocatable ::      pstress_field, defgrad, defgradold, cstress_field
+ real(pReal), dimension(:,:,:,:,:), allocatable ::      pstress_field, defgrad, defgradold
  
 ! variables storing information for spectral method
  complex(pReal), dimension(:,:,:,:,:), allocatable ::   workfft
@@ -84,9 +84,8 @@ program mpie_spectral
  integer*8, dimension(2,3,3) ::                         plan_fft
 
 ! convergence etc.
- real(pReal) err_div, err_stress, err_defgrad, err_div_temp
+ real(pReal) err_div, err_stress, err_defgrad, sigma0
- real(pReal) err_div_tol, err_stress_tol, err_stress_tolrel, err_defgrad_tol, sigma0
+ integer(pInt) ierr
- integer(pInt) itmax, ierr
  logical errmatinv
 
 ! loop variables etc.
@@ -105,18 +104,11 @@ program mpie_spectral
  
  N_l = 0_pInt; N_s = 0_pInt
  N_t = 0_pInt; N_n = 0_pInt
-
+ gotResolution =.false.; gotDimension =.false.; gotHomogenization = .false.
  resolution = 1_pInt; geomdimension = 0.0_pReal
  
- err_div_tol = 1.0e-4
- itmax = 250_pInt
- err_stress_tolrel=0.01
- err_defgrad_tol=1.0e-5
- err_defgrad_tol=1.0e-3 ! for test
  temperature = 300.0_pReal
 
- gotResolution =.false.; gotDimension =.false.; gotHomogenization = .false.
-
  if (IargC() /= 2) call IO_error(102)                     ! check for correct number of given arguments
 
 ! Reading the loadcase file and assign variables
@@ -171,16 +163,16 @@ program mpie_spectral
          do k = 1,9
            if (bc_maskvector(k)) valuevector(k) = IO_floatValue(line,posInput,j+k) ! assign values for the velocity gradient matrix
          enddo
-         bc_mask(:,:,1,i) = reshape(bc_maskvector,(/3,3/))
+         bc_mask(:,:,1,i) = transpose(reshape(bc_maskvector,(/3,3/)))
-         bc_velocityGrad(:,:,i) = reshape(valuevector,(/3,3/))
+         bc_velocityGrad(:,:,i) = math_transpose3x3(reshape(valuevector,(/3,3/)))
        case('s','stress')
          valuevector = 0.0_pReal
          forall (k = 1:9) bc_maskvector(k) = IO_stringValue(line,posInput,j+k) /= '#'
          do k = 1,9
            if (bc_maskvector(k)) valuevector(k) = IO_floatValue(line,posInput,j+k)  ! assign values for the bc_stress matrix
          enddo
-         bc_mask(:,:,2,i) = reshape(bc_maskvector,(/3,3/))
+         bc_mask(:,:,2,i) = transpose(reshape(bc_maskvector,(/3,3/)))
-         bc_stress(:,:,i) = reshape(valuevector,(/3,3/))
+         bc_stress(:,:,i) = math_transpose3x3(reshape(valuevector,(/3,3/)))
        case('t','time','delta')                                            ! increment time
            bc_timeIncrement(i) = IO_floatValue(line,posInput,j+1)
        case('n','incs','increments','steps')                               ! bc_steps
@@ -192,10 +184,10 @@ program mpie_spectral
 
  do i = 1, N_Loadcases
    if (any(bc_mask(:,:,1,i) == bc_mask(:,:,2,i))) call IO_error(47,i)     ! bc_mask consistency
-   print '(a,/,3(3(f12.6,x)/))','L',bc_velocityGrad(:,:,i)
+   print '(a,/,3(3(f12.6,x)/))','L',math_transpose3x3(bc_velocityGrad(:,:,i))
-   print '(a,/,3(3(f12.6,x)/))','bc_stress',bc_stress(:,:,i)
+   print '(a,/,3(3(f12.6,x)/))','bc_stress',math_transpose3x3(bc_stress(:,:,i))
-   print '(a,/,3(3(l,x)/))','bc_mask for velocitygrad',bc_mask(:,:,1,i)
+   print '(a,/,3(3(l,x)/))','bc_mask for velocitygrad',transpose(bc_mask(:,:,1,i))
-   print '(a,/,3(3(l,x)/))','bc_mask for stress',bc_mask(:,:,2,i)
+   print '(a,/,3(3(l,x)/))','bc_mask for stress',transpose(bc_mask(:,:,2,i))
    print *,'time',bc_timeIncrement(i)
    print *,'incs',bc_steps(i)
    print *, ''
@@ -253,7 +245,6 @@ program mpie_spectral
  allocate (gamma_hat(resolution(1)/2+1,resolution(2),resolution(3),3,3,3,3));    gamma_hat            = 0.0_pReal
  allocate (xi(resolution(1)/2+1,resolution(2),resolution(3),3));                 xi                   = 0.0_pReal
  allocate (pstress_field(resolution(1),resolution(2),resolution(3),3,3));        pstress_field        = 0.0_pReal
- allocate (cstress_field(resolution(1),resolution(2),resolution(3),3,3));        cstress_field        = 0.0_pReal
  allocate (defgrad(resolution(1),resolution(2),resolution(3),3,3));              defgrad              = 0.0_pReal
  allocate (defgradold(resolution(1),resolution(2),resolution(3),3,3));           defgradold           = 0.0_pReal
  allocate (ddefgrad(resolution(1),resolution(2),resolution(3)));                ddefgrad              = 0.0_pReal
@@ -365,7 +356,11 @@ program mpie_spectral
      enddo; enddo; enddo
 
      guessmode = 1.0_pReal                              ! keep guessing along former trajectory during same loadcase
+     if(all(bc_mask(:,:,1,loadcase))) then
+       calcmode = 1_pInt                                ! if no stress BC is given (calmode 0 is not needed)
+     else
        calcmode = 0_pInt                                ! start calculation of BC fulfillment
+     endif
      CPFEM_mode = 1_pInt                                ! winding forward
      iter = 0_pInt
      err_div= 2_pReal * err_div_tol                     ! go into loop 
@@ -374,12 +369,13 @@ program mpie_spectral
     
 !*************************************************************
 ! convergence loop
-     do while( iter <= itmax .and. &     
+     do while(iter < itmax .and. &     
              (err_div > err_div_tol .or. &
               err_stress > err_stress_tol .or. &
               err_defgrad > err_defgrad_tol))    
        iter = iter + 1
        print '(3(A,I5.5,tr2))', ' Loadcase = ',loadcase, ' Step = ',steps,'Iteration = ',iter
+       cstress_av = 0.0_pReal
 !*************************************************************
        
 ! adjust defgrad to fulfill BCs 
@@ -403,15 +399,12 @@ program mpie_spectral
                               cstress,dsde, pstress, dPdF)
            CPFEM_mode = 2_pInt
            pstress_field(i,j,k,:,:) = pstress
-           cstress_field(i,j,k,:,:) = math_mandel6to33(cstress)          
+           cstress_av = cstress_av + math_mandel6to33(cstress)          
          enddo; enddo; enddo
-   c066 = c066 * wgt
-   c0 = math_mandel66to3333(c066)
-
          
+         cstress_av = cstress_av * wgt
          do m = 1,3; do n = 1,3
            pstress_av(m,n) = sum(pstress_field(:,:,:,m,n)) * wgt   
-           cstress_av(m,n) = sum(cstress_field(:,:,:,m,n)) * wgt
            defgrad_av(m,n) = sum(defgrad(:,:,:,m,n)) * wgt
          enddo; enddo
          
@@ -437,12 +430,12 @@ program mpie_spectral
          enddo; enddo
          err_div = 2 * err_div_tol 
          err_defgrad = maxval(abs(mask_defgrad * (defgrad_av - defgradAim)))         
-         print '(a,/,3(3(f12.7,x)/))', ' Deformation Gradient:   ',defgrad_av(1:3,:)                           
+         print '(a,/,3(3(f12.7,x)/))', ' Deformation Gradient:   ',math_transpose3x3(defgrad_av)                           
-         print '(a,/,3(3(f10.4,x)/))', ' Cauchy Stress [MPa]: ',cstress_av(1:3,:)/1.e6        
+         print '(a,/,3(3(f10.4,x)/))', ' Cauchy Stress [MPa]: ',math_transpose3x3(cstress_av)/1.e6        
          print '(2(a,E8.2))', ' error stress               ',err_stress,'  Tol. = ', err_stress_tol 
          print '(2(a,E8.2))', ' error deformation gradient ',err_defgrad,'  Tol. = ', err_defgrad_tol*0.8  
          if(err_stress < err_stress_tol*0.8) then 
-           calcmode = 1
+           calcmode = 1_pInt
          endif
              
 ! Using the spectral method to calculate the change of deformation gradient, check divergence of stress field in fourier space        
@@ -464,7 +457,7 @@ program mpie_spectral
                               temperature,timeinc,ielem,1_pInt,&
                               cstress,dsde, pstress, dPdF)
            pstress_field(i,j,k,:,:) = pstress 
-           cstress_field(i,j,k,:,:) = math_mandel6to33(cstress)          
+           cstress_av = cstress_av + math_mandel6to33(cstress)          
          enddo; enddo; enddo       
          
          print *, 'Calculating equilibrium using spectral method'
@@ -485,14 +478,13 @@ program mpie_spectral
          enddo; enddo; enddo
          workfft(1,1,1,:,:) = defgrad_av - math_I3
          
-        ! err_div = err_div/real((prodnn/resolution(1)*(resolution(1)/2+1)), pReal)/sigma0       !weighting of error
          err_div = err_div/sigma0       !weighting of error
          
+         cstress_av = cstress_av * wgt
          do m = 1,3; do n = 1,3
             call dfftw_execute_dft_c2r(plan_fft(2,m,n), workfft(:,:,:,m,n),ddefgrad(:,:,:))
             defgrad(:,:,:,m,n) = defgrad(:,:,:,m,n) + ddefgrad * wgt
             pstress_av(m,n) = sum(pstress_field(:,:,:,m,n))*wgt
-            cstress_av(m,n) = sum(cstress_field(:,:,:,m,n))*wgt
             defgrad_av(m,n) = sum(defgrad(:,:,:,m,n))*wgt
             defgrad(:,:,:,m,n) = defgrad(:,:,:,m,n) + (defgradAim(m,n) - defgrad_av(m,n)) !anticipated target minus current state
          enddo; enddo

code/numerics.f90

@@ -46,11 +46,15 @@ real(pReal)                     relevantStrain, &                       ! strain
                                 volDiscrMod_RGC, &                      ! stiffness of RGC volume discrepancy (zero = without volume discrepancy constraint)
                                 volDiscrPow_RGC, &                      ! powerlaw penalty for volume discrepancy
 !* spectral parameters:
-                                rTol_defgradAvg                         ! relative tolerance for correction to deformation gradient aim
+                                err_div_tol, &                          ! error of divergence in fourier space
+                                err_stress_tol, &                       ! absolut stress error, will be computed from err_stress_tolrel (dont prescribe a value)
+                                err_stress_tolrel, &                    ! factor to multiply with highest stress to get err_stress_tol
+                                err_defgrad_tol                         ! tolerance for error of defgrad compared to prescribed defgrad
+integer(pInt)                   itmax , &                               ! maximum number of iterations
 
-                                !* Random seeding parameters: added <<<updated 27.08.2009>>>
+!* Random seeding parameters
-integer(pInt)                   fixedSeed                               ! fixed seeding for pseudo-random number generator
+                                fixedSeed                               ! fixed seeding for pseudo-random number generator
-! OpenMP variable
+!* OpenMP variable
 !$ integer(pInt)                mpieNumThreadsInt                       ! value stored in environment variable MPIE_NUM_THREADS
 
 
@@ -136,7 +140,10 @@ subroutine numerics_init()
   volDiscrPow_RGC         = 5.0
 
 !* spectral parameters:  
-  rTol_defgradAvg         = 1.0e-6
+  err_div_tol             = 1.0e-4  
+  err_defgrad_tol         = 1.0e-3
+  err_stress_tolrel       = 0.01
+  itmax                   = 20_pInt
 
 !* Random seeding parameters: added <<<updated 27.08.2009>>>
   fixedSeed               = 0_pInt
@@ -209,7 +216,7 @@ subroutine numerics_init()
         case ('integratorstiffness')
               integratorStiffness = IO_intValue(line,positions,2)
 
-!* RGC parameters: added <<<updated 17.12.2009>>>
+!* RGC parameters: 
         case ('atol_rgc')
               absTol_RGC = IO_floatValue(line,positions,2)
         case ('rtol_rgc')
@@ -238,10 +245,16 @@ subroutine numerics_init()
               volDiscrPow_RGC = IO_floatValue(line,positions,2)
 
 !* spectral parameters
-        case ('rTol_defgradAvg')
+        case ('err_div_tol')
-              rTol_defgradAvg = IO_floatValue(line,positions,2)
+              err_div_tol = IO_floatValue(line,positions,2)
+        case ('err_defgrad_tol')
+              err_defgrad_tol = IO_floatValue(line,positions,2)
+        case ('err_stress_tolrel')
+              err_stress_tolrel = IO_floatValue(line,positions,2)
+        case ('itmax')
+              itmax = IO_intValue(line,positions,2)
 
-!* Random seeding parameters: added <<<updated 27.08.2009>>>
+!* Random seeding parameters
         case ('fixed_seed')
               fixedSeed = IO_floatValue(line,positions,2)
       endselect
@@ -300,8 +313,10 @@ subroutine numerics_init()
   write(6,*)
 
 !* spectral parameters
-  write(6,'(a24,x,e8.1)') 'rTol_defgradAvg:        ',rTol_defgradAvg
+  write(6,'(a24,x,e8.1)') 'err_div_tol:             ',err_div_tol
-
+  write(6,'(a24,x,e8.1)') 'err_defgrad_tol:         ',err_defgrad_tol
+  write(6,'(a24,x,e8.1)') 'err_stress_tolrel:       ',err_stress_tolrel
+  write(6,'(a24,x,i8)')   'itmax:                   ',itmax
   write(6,*)
 
 !* Random seeding parameters
@@ -356,7 +371,10 @@ subroutine numerics_init()
   if (volDiscrPow_RGC <= 0.0_pReal)         call IO_error(289)
 
 !* spectral parameters
-  if (rTol_defgradAvg <= 0.0_pReal)         call IO_error(48)
+  if (err_div_tol <= 0.0_pReal)             call IO_error(48)
+  if (err_defgrad_tol <= 0.0_pReal)         call IO_error(48)
+  if (err_stress_tolrel <= 0.0_pReal)       call IO_error(48)
+  if (itmax <= 1.0_pInt)                    call IO_error(48)
   
   if (fixedSeed <= 0_pInt)                  write(6,'(a)') 'Random is random!'
 endsubroutine

processing/post/spectral_post.py

@@ -63,6 +63,8 @@ class MPIEspectral_result:
     self.N_increments =  self._keyedInt('increments')
     self.N_element_scalars = self._keyedInt('materialpoint_sizeResults')
     self.resolution = self._keyedPackedArray('resolution',3,'i')
+    #print self.resolution
+    #self.resolution = numpy.array([10,10,10],'i')
     self.N_nodes = (self.resolution[0]+1)*(self.resolution[1]+1)*(self.resolution[2]+1)
     self.N_elements = self.resolution[0]*self.resolution[1]*self.resolution[2]
     self.dimension = self._keyedPackedArray('dimension',3,'d')
@@ -127,9 +129,9 @@ class MPIEspectral_result:
 def readScalar(resolution,file,distance,startingPosition,offset):
   currentPosition = startingPosition+offset*8+4 - distance*8 # we add distance later on
   field = numpy.zeros([resolution[0],resolution[1],resolution[2]], 'd')
-  for x in range(0,resolution[0]):
-    for y in range(0,resolution[1]):
   for z in range(0,resolution[2]):
+    for y in range(0,resolution[1]):
+      for x in range(0,resolution[0]):
         currentPosition = currentPosition + distance*8
         p.file.seek(currentPosition)
         field[x][y][z]=struct.unpack('d',p.file.read(8))[0]
@@ -138,9 +140,9 @@ def readScalar(resolution,file,distance,startingPosition,offset):
 def readTensor(resolution,file,distance,startingPosition,offset):
   currentPosition = startingPosition+offset*8+4 - distance*8 # we add distance later on
   field = numpy.zeros([resolution[0],resolution[1],resolution[2],3,3], 'd')
-  for x in range(0,resolution[0]):
-    for y in range(0,resolution[1]):
   for z in range(0,resolution[2]):
+    for y in range(0,resolution[1]):
+      for x in range(0,resolution[0]):
         currentPosition = currentPosition + distance*8
         p.file.seek(currentPosition)
         for i in range(0,3):
@@ -148,6 +150,39 @@ def readTensor(resolution,file,distance,startingPosition,offset):
             field[x][y][z][i][j]=struct.unpack('d',p.file.read(8))[0]
   return field
 
+#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++   
+def calculateCauchyStress(p_stress,defgrad,res):
+#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 
+  c_stress = numpy.zeros([res[0],res[1],res[2],3,3],'d')
+  for z in range(res[2]):
+    for y in range(res[1]):
+      for x in range(res[0]):
+        jacobi =  numpy.linalg.det(defgrad[x,y,z])
+        c_stress[x,y,z] = numpy.dot(p_stress[x,y,z],numpy.transpose(defgrad[x,y,z]))/jacobi
+  return c_stress
+  
+#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++   
+def calculateVonMises(tensor,res):
+#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 
+  vonMises = numpy.zeros([res[0],res[1],res[2]],'d')
+  deviator = numpy.zeros([3,3],'d')
+  delta = numpy.zeros([3,3],'d')
+  delta[0,0] = 1.0
+  delta[1,1] = 1.0
+  delta[2,2] = 1.0
+  for z in range(res[2]):
+    for y in range(res[1]):
+      for x in range(res[0]): 
+       deviator = tensor[x,y,z] - 1.0/3.0*tensor[x,y,z,0,0]*tensor[x,y,z,1,1]*tensor[x,y,z,2,2]*delta
+       J_2 = deviator[0,0]*deviator[1,1]\
+           + deviator[1,1]*deviator[2,2]\
+           + deviator[0,0]*deviator[2,2]\
+           - (deviator[0,1])**2\
+           - (deviator[1,2])**2\
+           - (deviator[0,2])**2
+       vonMises[x,y,z] = numpy.sqrt(3*J_2)
+  return vonMises
+  
 #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 
 def mesh(res,geomdim,defgrad_av,centroids):
 #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 
@@ -256,8 +291,6 @@ def centroids(res,geomdimension,defgrad):
         parameter_coords=(2.0*numpy.array([i,j,k])-res+fones)/(res-fones)
         pos = (fones + parameter_coords)
         neg = (fones - parameter_coords)
-        if(k<3 and j<3 and i<3):
-          print i,j,k,':',pos, neg,'(',parameter_coords,')'
         centroids[i,j,k] =  ( cornerCoords[0,i,j,k] *neg[0]*neg[1]*neg[2]\
                             + cornerCoords[1,i,j,k] *pos[0]*neg[1]*neg[2]\
                             + cornerCoords[2,i,j,k] *pos[0]*pos[1]*neg[2]\
@@ -414,7 +447,8 @@ print 'Post Processing for Material subroutine for BVP solution using spectral m
 print '*********************************************************************************\n'
 
 #reading in the header of the results file
-p = MPIEspectral_result('32x32x32x100.spectralOut')
+name = 'dipl32'
+p = MPIEspectral_result(name+'.spectralOut')
 p.extrapolation('')
 print p
 
@@ -424,18 +458,33 @@ res_x=p.resolution[0]
 res_y=p.resolution[1]
 res_z=p.resolution[2]
 
+# for i in range(1,3):
+  # print('new step')
+  # c_pos = p.dataOffset + i*(p.N_element_scalars*8*p.N_elements + 8) #8 accounts for header&footer
+  # for j in range(p.N_element_scalars):
+  # #def readScalar(resolution,file,distance,startingPosition,offset):
+  # #currentPosition = startingPosition+offset*8+4 - distance*8 # we add distance later on
+  # #field = numpy.zeros([resolution[0],resolution[1],resolution[2]], 'd')
+  # #for z in range(0,resolution[2]):
+    # #for y in range(0,resolution[1]):
+      # #for x in range(0,resolution[0]):
+    # currentPosition = c_pos + j*8 +4
+    # p.file.seek(currentPosition)
+    # print(struct.unpack('d',p.file.read(8)))
 
     
-for i in range(120,121):
+
+for i in range(40,46):
   c_pos = p.dataOffset + i*(p.N_element_scalars*8*p.N_elements + 8) #8 accounts for header&footer
   defgrad = readTensor(p.resolution,p.file,p.N_element_scalars,c_pos,16)         
-  grain = readScalar(p.resolution,p.file,p.N_element_scalars,c_pos,7)
+  p_stress = readTensor(p.resolution,p.file,p.N_element_scalars,c_pos,52)         
-  centroids, defgrad_av = centroids(p.resolution,p.dimension,defgrad)
+  #c_stress = calculateCauchyStress(p_stress,defgrad,p.resolution)
-  print centroids.shape, defgrad_av.shape
+  #grain = calculateVonMises(c_stress,p.resolution)
-  ms = mesh(p.resolution,p.dimension,defgrad_av,centroids)
+  centroids_coord, defgrad_av = centroids(p.resolution,p.dimension,defgrad)
-  writeVtkAscii('mesh-%i.vtk'%i,ms,grain,p.resolution)
+  ms = mesh(p.resolution,p.dimension,defgrad_av,centroids_coord)
-  writeVtkAsciidefgrad_av('box-%i.vtk'%i,p.dimension,defgrad_av)
+  writeVtkAscii(name+'-mesh-%i.vtk'%i,ms,p_stress[:,:,:,1,2],p.resolution)
-  writeVtkAsciiDots('points-%i.vtk'%i,centroids,grain,p.resolution)
+  writeVtkAsciidefgrad_av(name+'-box-%i.vtk'%i,p.dimension,defgrad_av)
+  #writeVtkAsciiDots(name+'-points-%i.vtk'%i,centroids_coord,grain,p.resolution)
   sys.stdout.flush()