some polishing for single precision version.

Now only prec and mpie_spectral are needed in different versions

code/makefile

@@ -14,16 +14,12 @@ precision :=double
 ifeq ($(precision),single)
 cpspectral_single.exe: mpie_spectral_single.o CPFEM.a
 	ifort -openmp -o cpspectral_single.exe mpie_spectral_single.o CPFEM.a libfftw3f_threads.a libfftw3f.a constitutive.a advanced.a basics.a -lpthread
-	rm *.mod
-
 mpie_spectral_single.o: mpie_spectral_single.f90 CPFEM.o 
 	ifort -openmp -c -O3 -heap-arrays 500000000 mpie_spectral_single.f90 
 
 else
 cpspectral_double.exe: mpie_spectral.o CPFEM.a
 	ifort -openmp -o cpspectral_double.exe mpie_spectral.o CPFEM.a libfftw3_threads.a libfftw3.a constitutive.a advanced.a basics.a -lpthread
-	rm *.mod
-
 mpie_spectral.o: mpie_spectral.f90 CPFEM.o 
 	ifort -openmp -c -O3 -heap-arrays 500000000 mpie_spectral.f90 
 endif
@@ -66,31 +62,20 @@ constitutive_phenopowerlaw.o: constitutive_phenopowerlaw.f90 basics.a advanced.a
 	ifort -openmp -c -O3 -heap-arrays 500000000 constitutive_phenopowerlaw.f90
 
 
-ifeq ($(precision),single)
+
-advanced.a: lattice.o
-	ar rc advanced.a FEsolving.o mesh_single.o material.o lattice.o
-else
 advanced.a: lattice.o
 	ar rc advanced.a FEsolving.o mesh.o material.o lattice.o
-endif
+
 
 lattice.o: lattice.f90 material.o
 	ifort -openmp -c -O3 -heap-arrays 500000000 lattice.f90
-ifeq ($(precision),single)
-material.o: material.f90 mesh_single.o
-	ifort -openmp -c -O3 -heap-arrays 500000000 material.f90
-mesh_single.o: mesh_single.f90 FEsolving.o
-	ifort -openmp -c -O3 -heap-arrays 500000000 mesh_single.f90
-else
 material.o: material.f90 mesh.o
 	ifort -openmp -c -O3 -heap-arrays 500000000 material.f90
 mesh.o: mesh.f90 FEsolving.o
 	ifort -openmp -c -O3 -heap-arrays 500000000 mesh.f90
-endif
 FEsolving.o: FEsolving.f90 basics.a
 	ifort -openmp -c -O3 -heap-arrays 500000000 FEsolving.f90
 
-
 ifeq ($(precision),single)
 basics.a: debug.o math.o
 	ar rc basics.a debug.o math.o numerics.o IO.o mpie_spectral_interface.o prec_single.o
@@ -101,7 +86,6 @@ endif
 
 debug.o: debug.f90 numerics.o
 	ifort -openmp -c -O3 debug.f90
-
 math.o: math.f90 numerics.o
 	ifort -openmp -c -O3 math.f90
 
@@ -109,6 +93,7 @@ numerics.o: numerics.f90 IO.o
 	ifort -openmp -c -O3 numerics.f90
 IO.o: IO.f90 mpie_spectral_interface.o
 	ifort -openmp -c -O3 IO.f90
+
 ifeq ($(precision),single)
 mpie_spectral_interface.o: mpie_spectral_interface.f90 prec_single.o
 	ifort -openmp -c -O3 mpie_spectral_interface.f90
@@ -124,6 +109,7 @@ endif
 
 clean: 
 	rm -rf *.o 
+	rm -rf *.mod
 	rm -rf *.exe 
 	rm -rf basics.a
 	rm -rf advanced.a

code/math.f90

@@ -118,7 +118,7 @@ real(pReal), dimension(4,36), parameter :: math_symOperations = &
 !**************************************************************************
  SUBROUTINE math_init ()
 
- use prec,     only: pReal,pInt
+ use prec,     only: pReal,pInt,tol_math_check
  use numerics, only: fixedSeed
  use IO,       only: IO_error
  implicit none
@@ -155,28 +155,28 @@ real(pReal), dimension(4,36), parameter :: math_symOperations = &
  q = math_qRnd();
  axisangle = math_QuaternionToAxisAngle(q);
  q2 = math_AxisAngleToQuaternion(axisangle(1:3),axisangle(4))
- if ( any(abs( q-q2) > 1.0e-8 ) .and. &
+ if ( any(abs( q-q2) > tol_math_check ) .and. &
-      any(abs(-q-q2) > 1.0e-8 ) ) &
+      any(abs(-q-q2) > tol_math_check ) ) &
    call IO_error(670)
  
  ! +++ q -> R -> q  +++
  R = math_QuaternionToR(q);
  q2 = math_RToQuaternion(R)
- if ( any(abs( q-q2) > 1.0e-8 ) .and. &
+ if ( any(abs( q-q2) > tol_math_check ) .and. &
-      any(abs(-q-q2) > 1.0e-8 ) ) &
+      any(abs(-q-q2) > tol_math_check ) ) &
    call IO_error(671)
  
  ! +++ q -> euler -> q  +++
  Eulers = math_QuaternionToEuler(q);
  q2 = math_EulerToQuaternion(Eulers)
- if ( any(abs( q-q2) > 1.0e-8 ) .and. &
+ if ( any(abs( q-q2) > tol_math_check ) .and. &
-      any(abs(-q-q2) > 1.0e-8 ) ) &
+      any(abs(-q-q2) > tol_math_check ) ) &
    call IO_error(672)
  
  ! +++ R -> euler -> R  +++
  Eulers = math_RToEuler(R);
  R2 = math_EulerToR(Eulers)
- if ( any(abs( R-R2) > 1.0e-8 ) ) &
+ if ( any(abs( R-R2) > tol_math_check ) ) &
    call IO_error(673)
  
  

code/prec.f90

@@ -9,6 +9,8 @@
  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
+ real(pReal), parameter :: tol_math_check = 1.0e-8_pReal
+ real(pReal), parameter :: tol_gravityNodePos = 1.0e-100_pReal
  
  type :: p_vec
      real(pReal), dimension(:), pointer :: p
@@ -17,6 +19,8 @@
 CONTAINS
 
 subroutine prec_init
+ implicit none
+
  write(6,*)
  write(6,*) '<<<+-  prec init  -+>>>'
  write(6,*) '$Id$'

code/prec_single.f90

@@ -9,6 +9,8 @@
  integer, parameter :: pReal = selected_real_kind(6,37)       ! 6 significant digits, up to 1e+-37
  integer, parameter :: pInt  = selected_int_kind(9)           ! up to +- 1e9
  integer, parameter :: pLongInt  = 8                          ! should be 64bit
+ real(pReal), parameter :: tol_math_check = 1.0e-6_pReal
+ real(pReal), parameter :: tol_gravityNodePos = 1.0e-36_pReal
  
  type :: p_vec
      real(pReal), dimension(:), pointer :: p

processing/post/spectral_post.py

@@ -324,45 +324,15 @@ for i in range(p.N_element_scalars):
   c_pos = p.dataOffset + i*8.0 + 4.0
   p.file.seek(c_pos)
   print(i, struct.unpack('d',p.file.read(8)))
-for i in range(1,2):
-  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,7)
-  defgrad_av = postprocessingMath.tensor_avg(res_x,res_y,res_z,defgrad)
-  centroids_coord = postprocessingMath.deformed(res_x,res_y,res_z,p.dimension,defgrad,defgrad_av)
-  undeformed = postprocessingMath.mesh(p.resolution[0],p.resolution[1],p.resolution[2],p.dimension,defgrad_av,centroids_coord)
-  #writeVtkAscii(name+'-mesh-undeformed.vtk',undeformed,defgrad[:,:,:,1,2],p.resolution)
 
-for i in range(240,241):
+for i in range(200,201): # define here the steps
   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,7)
+  defgrad = readTensor(p.resolution,p.file,p.N_element_scalars,c_pos,7) #define here the position of the deformation gradient
+  rotation = readTensor(p.resolution,p.file,p.N_element_scalars,c_pos,7) #define here the position of the tensor
   defgrad_av = postprocessingMath.tensor_avg(res_x,res_y,res_z,defgrad)
-  #defgrad = numpy.zeros([p.resolution[0],p.resolution[1],p.resolution[2],3,3], 'd')
+  centroids_coord = postprocessingMath.deformed_fft(res_x,res_y,res_z,p.dimension,defgrad,defgrad_av,1.0)
-  #for z in range(p.resolution[2]):
-  #   defgrad[:,:,z,1,2] = (2.0*z)/(p.resolution[2]-1.0)+ 3.8*math.sin(z*20.0/(p.resolution[2]-1.0)*2*math.pi)
-  #   defgrad[:,:,z,0,2] = (2.0*z)/(p.resolution[2]-1.0)+ 5.0*math.cos(z/(p.resolution[2]-1.0)*2*math.pi)
-  #defgrad[:,:,:,0,0] = 1.0
-  #defgrad[:,:,:,1,1] = 1.0
-  #defgrad[:,:,:,2,2] = 1.0
-  #logstrain = postprocessingMath.logstrain_mat(res_x,res_y,res_z,defgrad)  
-  #logstrain2 = postprocessingMath.logstrain_spat(res_x,res_y,res_z,defgrad)
-  #p_stress = readTensor(p.resolution,p.file,p.N_element_scalars,c_pos,52)         
-  #c_stress = postprocessingMath.calculate_cauchy(res_x,res_y,res_z,defgrad,p_stress)
-  #vm = postprocessingMath.calculate_mises(res_x,res_y,res_z,c_stress)
-  #defgrad_av = postprocessingMath.tensor_avg(res_x,res_y,res_z,defgrad)
-  #subroutine inverse_reconstruction(res_x,res_y,res_z,reference_configuration,current_configuration,defgrad)
-  centroids_coord = postprocessingMath.deformed(res_x,res_y,res_z,p.dimension,defgrad,defgrad_av)
-  deformed = postprocessingMath.mesh(p.resolution[0],p.resolution[1],p.resolution[2],p.dimension,defgrad_av,centroids_coord)
-  #writeVtkAscii(name+'-mesh-deformed.vtk',deformed,defgrad[:,:,:,1,2],p.resolution)
-  defgrad = postprocessingMath.inverse_reconstruction(res_x,res_y,res_z,undeformed,deformed)
-  defgrad_av = postprocessingMath.tensor_avg(res_x,res_y,res_z,defgrad)
-  centroids_coord = postprocessingMath.deformed(res_x,res_y,res_z,p.dimension,defgrad,defgrad_av)
-  centroids_coord2 = postprocessingMath.deformed_fft(res_x,res_y,res_z,p.dimension,defgrad,defgrad_av,1.0)
   ms = postprocessingMath.mesh(p.resolution[0],p.resolution[1],p.resolution[2],p.dimension,defgrad_av,centroids_coord)
-  ms2 = postprocessingMath.mesh(p.resolution[0],p.resolution[1],p.resolution[2],p.dimension,defgrad_av,centroids_coord2)
+  writeVtkAscii(name+'-mesh-fft-%s.vtk'%i,ms,defgrad[:,:,:,1,2],p.resolution)
-  writeVtkAscii(name+'-mesh-usual-%s.vtk'%i,ms,defgrad[:,:,:,1,2],p.resolution)
-  writeVtkAscii(name+'-mesh-fft-%s.vtk'%i,ms2,defgrad[:,:,:,1,2],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()