From 8cae4d609a874bbe6ff96e05b672573636ffae29 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Fri, 25 Feb 2011 12:41:46 +0000 Subject: [PATCH] some polishing for single precision version. Now only prec and mpie_spectral are needed in different versions --- code/makefile | 22 ++++-------------- code/math.f90 | 16 ++++++------- code/prec.f90 | 6 ++++- code/prec_single.f90 | 4 +++- processing/post/spectral_post.py | 40 ++++---------------------------- 5 files changed, 25 insertions(+), 63 deletions(-) diff --git a/code/makefile b/code/makefile index 01bc539bb..602bb5d94 100644 --- a/code/makefile +++ b/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 diff --git a/code/math.f90 b/code/math.f90 index 401f28b7d..6acc4bfd6 100644 --- a/code/math.f90 +++ b/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. & - any(abs(-q-q2) > 1.0e-8 ) ) & + if ( any(abs( q-q2) > tol_math_check ) .and. & + 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. & - any(abs(-q-q2) > 1.0e-8 ) ) & + if ( any(abs( q-q2) > tol_math_check ) .and. & + 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. & - any(abs(-q-q2) > 1.0e-8 ) ) & + if ( any(abs( q-q2) > tol_math_check ) .and. & + 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) diff --git a/code/prec.f90 b/code/prec.f90 index 7a1e79e5e..a3b21a273 100644 --- a/code/prec.f90 +++ b/code/prec.f90 @@ -9,7 +9,9 @@ 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 end type p_vec @@ -17,6 +19,8 @@ CONTAINS subroutine prec_init + implicit none + write(6,*) write(6,*) '<<<+- prec init -+>>>' write(6,*) '$Id$' diff --git a/code/prec_single.f90 b/code/prec_single.f90 index f3bfb02f3..ae86adf11 100644 --- a/code/prec_single.f90 +++ b/code/prec_single.f90 @@ -9,7 +9,9 @@ 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 end type p_vec diff --git a/processing/post/spectral_post.py b/processing/post/spectral_post.py index ebe755c23..fcbb1c965 100644 --- a/processing/post/spectral_post.py +++ b/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') - #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) + centroids_coord = 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-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) + writeVtkAscii(name+'-mesh-fft-%s.vtk'%i,ms,defgrad[:,:,:,1,2],p.resolution) sys.stdout.flush()