From fec2c14a4e8351758792589284d6f4a260e9fc85 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Mon, 31 Jan 2011 17:07:42 +0000 Subject: [PATCH] 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 | 2 +- code/mpie_spectral.f90 | 74 +++++++++++++----------------- code/numerics.f90 | 42 ++++++++++++----- processing/post/spectral_post.py | 79 ++++++++++++++++++++++++++------ 4 files changed, 128 insertions(+), 69 deletions(-) diff --git a/code/IO.f90 b/code/IO.f90 index 3d1da4a0c..60e57175d 100644 --- a/code/IO.f90 +++ b/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) diff --git a/code/mpie_spectral.f90 b/code/mpie_spectral.f90 index cbbc18549..1b11e6cec 100644 --- a/code/mpie_spectral.f90 +++ b/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_tol, err_stress_tol, err_stress_tolrel, err_defgrad_tol, sigma0 - integer(pInt) itmax, ierr + real(pReal) err_div, err_stress, err_defgrad, sigma0 + integer(pInt) 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_velocityGrad(:,:,i) = reshape(valuevector,(/3,3/)) + bc_mask(:,:,1,i) = transpose(reshape(bc_maskvector,(/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_stress(:,:,i) = reshape(valuevector,(/3,3/)) + bc_mask(:,:,2,i) = transpose(reshape(bc_maskvector,(/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)/))','bc_stress',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 stress',bc_mask(:,:,2,i) + print '(a,/,3(3(f12.6,x)/))','L',math_transpose3x3(bc_velocityGrad(:,:,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',transpose(bc_mask(:,:,1,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 @@ -364,26 +355,31 @@ program mpie_spectral defgradold(i,j,k,:,:) = temp33_Real enddo; enddo; enddo - guessmode = 1.0_pReal ! keep guessing along former trajectory during same loadcase - calcmode = 0_pInt ! start calculation of BC fulfillment - CPFEM_mode = 1_pInt ! winding forward + 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 - defgradAimCorr = 0.0_pReal ! reset damping calculation + err_div= 2_pReal * err_div_tol ! go into loop + defgradAimCorr = 0.0_pReal ! reset damping calculation damper = damper * 0.9_pReal !************************************************************* ! 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 - select case (calcmode) + select case (calcmode) case (0) print *, 'Update Stress Field (constitutive evaluation P(F))' ielem = 0_pInt @@ -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 + pstress_av(m,n) = sum(pstress_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(f10.4,x)/))', ' Cauchy Stress [MPa]: ',cstress_av(1:3,:)/1.e6 + print '(a,/,3(3(f12.7,x)/))', ' Deformation Gradient: ',math_transpose3x3(defgrad_av) + 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 diff --git a/code/numerics.f90 b/code/numerics.f90 index d74a14609..74b542ffa 100644 --- a/code/numerics.f90 +++ b/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 <<>> -integer(pInt) fixedSeed ! fixed seeding for pseudo-random number generator -! OpenMP variable +!* Random seeding parameters + fixedSeed ! fixed seeding for pseudo-random number generator +!* 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 <<>> fixedSeed = 0_pInt @@ -209,7 +216,7 @@ subroutine numerics_init() case ('integratorstiffness') integratorStiffness = IO_intValue(line,positions,2) -!* RGC parameters: added <<>> +!* 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') - rTol_defgradAvg = IO_floatValue(line,positions,2) + case ('err_div_tol') + 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 <<>> +!* 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 diff --git a/processing/post/spectral_post.py b/processing/post/spectral_post.py index 1894674a0..4372c6bee 100644 --- a/processing/post/spectral_post.py +++ b/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 z in range(0,resolution[2]): for y in range(0,resolution[1]): - for z in range(0,resolution[2]): + 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,15 +140,48 @@ 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 z in range(0,resolution[2]): for y in range(0,resolution[1]): - for z in range(0,resolution[2]): + for x in range(0,resolution[0]): currentPosition = currentPosition + distance*8 p.file.seek(currentPosition) for i in range(0,3): for j in range(0,3): 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) - centroids, defgrad_av = centroids(p.resolution,p.dimension,defgrad) - print centroids.shape, defgrad_av.shape - ms = mesh(p.resolution,p.dimension,defgrad_av,centroids) - writeVtkAscii('mesh-%i.vtk'%i,ms,grain,p.resolution) - writeVtkAsciidefgrad_av('box-%i.vtk'%i,p.dimension,defgrad_av) - writeVtkAsciiDots('points-%i.vtk'%i,centroids,grain,p.resolution) + p_stress = readTensor(p.resolution,p.file,p.N_element_scalars,c_pos,52) + #c_stress = calculateCauchyStress(p_stress,defgrad,p.resolution) + #grain = calculateVonMises(c_stress,p.resolution) + centroids_coord, defgrad_av = centroids(p.resolution,p.dimension,defgrad) + ms = mesh(p.resolution,p.dimension,defgrad_av,centroids_coord) + writeVtkAscii(name+'-mesh-%i.vtk'%i,ms,p_stress[:,:,:,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()