From 50db944c0c1a26344c9253bc953180fbfa1c3b31 Mon Sep 17 00:00:00 2001 From: Pratheek Shanthraj Date: Wed, 13 Nov 2013 19:21:35 +0000 Subject: [PATCH] reworked phase field interface to damask spectral solvers. now specify 'thermal a b c' or 'fracture a b c' to activate either phase field where a b c are the initial value, diffusion coefficient and mobility respectively. Right now only thermal and fracture phase fields implemented and only in the basic petsc solver --- code/DAMASK_spectral_driver.f90 | 43 ++++- code/DAMASK_spectral_solverBasicPETSc.f90 | 216 ++++++++++++++-------- code/DAMASK_spectral_utilities.f90 | 75 +------- 3 files changed, 189 insertions(+), 145 deletions(-) diff --git a/code/DAMASK_spectral_driver.f90 b/code/DAMASK_spectral_driver.f90 index 668da55fd..ab83105be 100644 --- a/code/DAMASK_spectral_driver.f90 +++ b/code/DAMASK_spectral_driver.f90 @@ -75,8 +75,9 @@ program DAMASK_spectral_Driver geomSize, & tBoundaryCondition, & tSolutionState, & + phaseFieldDataBin, & cutBack, & - utilities_temperatureUpdate + maxPhaseFields use DAMASK_spectral_SolverBasic #ifdef PETSc use DAMASK_spectral_SolverBasicPETSC @@ -89,14 +90,17 @@ program DAMASK_spectral_Driver real(pReal), dimension (3,3) :: rotation = math_I3 !< rotation of BC type(tBoundaryCondition) :: P, & !< stress BC deformation !< deformation BC (Fdot or L) + type(phaseFieldDataBin) :: phaseFieldData(maxPhaseFields) real(pReal) :: time = 0.0_pReal, & !< length of increment - temperature = 300.0_pReal, & !< isothermal starting conditions + temperature = 300.0_pReal, & !< isothermal starting condition density = 0.0_pReal !< density integer(pInt) :: incs = 0_pInt, & !< number of increments outputfrequency = 1_pInt, & !< frequency of result writes restartfrequency = 0_pInt, & !< frequency of restart writes logscale = 0_pInt !< linear/logarithmic time inc flag - logical :: followFormerTrajectory = .true. !< follow trajectory of former loadcase + logical :: followFormerTrajectory = .true., & !< follow trajectory of former loadcase + thermal_active = .false., & !< activate thermal phase field + fracture_active = .false. !< activate fracture phase field end type tLoadCase !-------------------------------------------------------------------------------------------------- @@ -112,7 +116,8 @@ program DAMASK_spectral_Driver integer(pInt) :: & N_t = 0_pInt, & !< # of time indicators found in load case file N_n = 0_pInt, & !< # of increment specifiers found in load case file - N_def = 0_pInt !< # of rate of deformation specifiers found in load case file + N_def = 0_pInt, & !< # of rate of deformation specifiers found in load case file + nActivePhaseFields character(len=65536) :: & line @@ -230,8 +235,22 @@ program DAMASK_spectral_Driver loadCases(currentLoadCase)%P%values = math_plain9to33(temp_valueVector) case('t','time','delta') ! increment time loadCases(currentLoadCase)%time = IO_floatValue(line,positions,i+1_pInt) - case('temp','temperature') ! starting temperature + case('temperature') loadCases(currentLoadCase)%temperature = IO_floatValue(line,positions,i+1_pInt) + case('thermal') ! starting temperature, conductivity and mobility + loadCases(:)%phaseFieldData(1)%label = 'thermal' + loadCases(:)%phaseFieldData(1)%active = .true. + loadCases(:)%phaseFieldData(1)%phaseField0 = 300.0_pReal ! initialize to meaningful value if not defined + loadCases(currentLoadCase)%phaseFieldData(1)%phaseField0 = IO_floatValue(line,positions,i+1_pInt) + loadCases(currentLoadCase)%phaseFieldData(1)%diffusion = IO_floatValue(line,positions,i+2_pInt) + loadCases(currentLoadCase)%phaseFieldData(1)%mobility = IO_floatValue(line,positions,i+3_pInt) + case('fracture') ! starting damage, diffusion and mobility + loadCases(:)%phaseFieldData(2)%label = 'fracture' + loadCases(:)%phaseFieldData(2)%active = .true. + loadCases(:)%phaseFieldData(2)%phaseField0 = 1.0_pReal ! initialize to meaningful value if not defined + loadCases(currentLoadCase)%phaseFieldData(2)%phaseField0 = IO_floatValue(line,positions,i+1_pInt) + loadCases(currentLoadCase)%phaseFieldData(2)%diffusion = IO_floatValue(line,positions,i+2_pInt) + loadCases(currentLoadCase)%phaseFieldData(2)%mobility = IO_floatValue(line,positions,i+3_pInt) case('den','density') ! starting density loadCases(currentLoadCase)%density = IO_floatValue(line,positions,i+1_pInt) case('n','incs','increments','steps') ! number of increments @@ -271,6 +290,14 @@ program DAMASK_spectral_Driver end select enddo; enddo close(myUnit) + ! reorder phase field data to remove redundant non-active fields + nActivePhaseFields = 0_pInt + do i = 1, maxPhaseFields + if (loadCases(1)%phaseFieldData(i)%active) then + nActivePhaseFields = nActivePhaseFields + 1_pInt + loadCases(:)%phaseFieldData(nActivePhaseFields) = loadCases(:)%phaseFieldData(i) + endif + enddo !-------------------------------------------------------------------------------------------------- ! consistency checks and output of load case @@ -336,7 +363,7 @@ program DAMASK_spectral_Driver call basic_init(loadCases(1)%temperature) #ifdef PETSc case (DAMASK_spectral_SolverBasicPETSc_label) - call basicPETSc_init(loadCases(1)%temperature) + call basicPETSc_init(loadCases(1)%temperature,nActivePhaseFields,loadCases(1)%phaseFieldData(1:nActivePhaseFields)) case (DAMASK_spectral_SolverAL_label) if(iand(debug_level(debug_spectral),debug_levelBasic)/= 0) & call IO_warning(42_pInt, ext_msg='debug Divergence') @@ -467,7 +494,9 @@ program DAMASK_spectral_Driver F_BC = loadCases(currentLoadCase)%deformation, & temperature_bc = loadCases(currentLoadCase)%temperature, & rotation_BC = loadCases(currentLoadCase)%rotation, & - density = loadCases(currentLoadCase)%density) + density = loadCases(currentLoadCase)%density, & + nActivePhaseFields = nActivePhaseFields, & + phaseFieldData = loadCases(1)%phaseFieldData(1:nActivePhaseFields)) case (DAMASK_spectral_SolverAL_label) solres = AL_solution (& incInfo,guess,timeinc,timeIncOld,remainingLoadCaseTime, & diff --git a/code/DAMASK_spectral_solverBasicPETSc.f90 b/code/DAMASK_spectral_solverBasicPETSc.f90 index 2af41b462..40f10f965 100644 --- a/code/DAMASK_spectral_solverBasicPETSc.f90 +++ b/code/DAMASK_spectral_solverBasicPETSc.f90 @@ -31,7 +31,9 @@ module DAMASK_spectral_SolverBasicPETSc use math, only: & math_I3 use DAMASK_spectral_Utilities, only: & - tSolutionState + tSolutionState, & + phaseFieldDataBin, & + maxPhaseFields implicit none private @@ -49,6 +51,8 @@ module DAMASK_spectral_SolverBasicPETSc real(pReal) :: timeincOld real(pReal) :: temperature real(pReal) :: density + integer(pInt) :: nActivePhaseFields + type(phaseFieldDataBin) :: phaseFieldData(maxPhaseFields) end type tSolutionParams type(tSolutionParams), private :: params @@ -62,9 +66,10 @@ module DAMASK_spectral_SolverBasicPETSc !-------------------------------------------------------------------------------------------------- ! common pointwise data real(pReal), private, dimension(:,:,:,:,:), allocatable :: F_lastInc, Fdot, F_lastInc2 - real(pReal), private, dimension(:,:,:), allocatable :: & - heatSource_lastInc, & + real(pReal), private, dimension(:,:,:,:), allocatable :: & + phaseFieldRHS_lastInc, & phaseField_lastInc, & + phaseFieldRHS, & phaseFieldDot complex(pReal), private, dimension(:,:,:,:,:), allocatable :: inertiaField_fourier @@ -82,7 +87,8 @@ module DAMASK_spectral_SolverBasicPETSc C_volAvgLastInc = 0.0_pReal, & !< previous volume average stiffness C_minMaxAvg = 0.0_pReal, & !< current (min+max)/2 stiffness S = 0.0_pReal !< current compliance (filled up with zeros) - real(pReal), private :: err_stress, err_div, err_divPrev, err_divDummy, err_phaseField, phaseField_Avg + real(pReal), private :: err_stress, err_div, err_divPrev, err_divDummy + real(pReal), private, dimension(:), allocatable :: err_phaseField, phaseField_Avg logical, private :: ForwardData integer(pInt), private :: & totalIter = 0_pInt !< total iteration in current increment @@ -115,7 +121,7 @@ contains !-------------------------------------------------------------------------------------------------- !> @brief allocates all neccessary fields and fills them with data, potentially from restart info !-------------------------------------------------------------------------------------------------- -subroutine basicPETSc_init(temperature) +subroutine basicPETSc_init(temperature,nActivePhaseFields,phaseFieldData) use, intrinsic :: iso_fortran_env ! to get compiler_version and compiler_options (at least for gfortran >4.6 at the moment) use IO, only: & IO_intOut, & @@ -144,14 +150,14 @@ subroutine basicPETSc_init(temperature) math_invSym3333 implicit none - real(pReal), intent(inout) :: & - temperature + integer(pInt), intent(in) :: nActivePhaseFields + type(phaseFieldDataBin), intent(in) :: phaseFieldData(nActivePhaseFields) + real(pReal), intent(inOut) :: temperature #include #include #include real(pReal), dimension(:,:,:,:,:), allocatable :: P PetscScalar, dimension(:,:,:,:), pointer :: xx_psc, F - PetscScalar, dimension(:,:,:), pointer :: phaseField PetscErrorCode :: ierr PetscObject :: dummy real(pReal), dimension(3,3) :: & @@ -159,23 +165,27 @@ subroutine basicPETSc_init(temperature) real(pReal), dimension(3,3,3,3) :: & temp3333_Real = 0.0_pReal KSP :: ksp + integer(pInt) :: i call Utilities_init() write(6,'(/,a)') ' <<<+- DAMASK_spectral_solverBasicPETSc init -+>>>' write(6,'(a)') ' $Id$' write(6,'(a15,a)') ' Current time: ',IO_timeStamp() #include "compilation_info.f90" - - allocate (P (3,3,grid(1),grid(2),grid(3)),source = 0.0_pReal) + !-------------------------------------------------------------------------------------------------- ! allocate global fields + allocate (P (3,3,grid(1),grid(2),grid(3)),source = 0.0_pReal) allocate (F_lastInc (3,3,grid(1),grid(2),grid(3)),source = 0.0_pReal) allocate (F_lastInc2(3,3,grid(1),grid(2),grid(3)),source = 0.0_pReal) allocate (Fdot (3,3,grid(1),grid(2),grid(3)),source = 0.0_pReal) allocate (inertiaField_fourier (grid1Red,grid(2),grid(3),3,3),source = cmplx(0.0_pReal,0.0_pReal,pReal)) - allocate (heatSource_lastInc (grid(1),grid(2),grid(3)),source = 0.0_pReal) - allocate (phaseField_lastInc(grid(1),grid(2),grid(3)),source = 0.0_pReal) - allocate (phaseFieldDot (grid(1),grid(2),grid(3)),source = 0.0_pReal) + allocate (phaseFieldRHS_lastInc (nActivePhaseFields,grid(1),grid(2),grid(3)),source = 0.0_pReal) + allocate (phaseField_lastInc (nActivePhaseFields,grid(1),grid(2),grid(3)),source = 0.0_pReal) + allocate (phaseFieldDot (nActivePhaseFields,grid(1),grid(2),grid(3)),source = 0.0_pReal) + allocate (phaseFieldRHS (nActivePhaseFields,grid(1),grid(2),grid(3)),source = 0.0_pReal) + allocate (err_phaseField(nActivePhaseFields), source = 0.0_pReal) + allocate (phaseField_Avg(nActivePhaseFields), source = 0.0_pReal) !-------------------------------------------------------------------------------------------------- ! initialize solver specific parts of PETSc @@ -183,31 +193,29 @@ subroutine basicPETSc_init(temperature) call DMDACreate3d(PETSC_COMM_WORLD, & DMDA_BOUNDARY_NONE, DMDA_BOUNDARY_NONE, DMDA_BOUNDARY_NONE, & DMDA_STENCIL_BOX,grid(1),grid(2),grid(3),PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE, & - 10,1,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,da,ierr) + 9+nActivePhaseFields,1,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,da,ierr) CHKERRQ(ierr) call DMCreateGlobalVector(da,solution_vec,ierr); CHKERRQ(ierr) - !call DMDASNESSetFunctionLocal(da,INSERT_VALUES,BasicPETSC_formResidual,dummy,ierr); CHKERRQ(ierr) ! needed for newer versions of petsc - call DMDASetLocalFunction(da,BasicPETSC_formResidual,ierr); CHKERRQ(ierr) + call DMDASNESSetFunctionLocal(da,INSERT_VALUES,BasicPETSC_formResidual,dummy,ierr); CHKERRQ(ierr) ! needed for newer versions of petsc + !call DMDASetLocalFunction(da,BasicPETSC_formResidual,ierr); CHKERRQ(ierr) call SNESSetDM(snes,da,ierr); CHKERRQ(ierr) - call SNESSetConvergenceTest(snes,BasicPETSC_converged,dummy,PETSC_NULL_FUNCTION,ierr) - CHKERRQ(ierr) + call SNESSetConvergenceTest(snes,BasicPETSC_converged,dummy,PETSC_NULL_FUNCTION,ierr); CHKERRQ(ierr) call SNESGetKSP(snes,ksp,ierr); CHKERRQ(ierr) - call KSPSetConvergenceTest(ksp,BasicPETSC_convergedKSP,dummy,PETSC_NULL_FUNCTION,ierr) - CHKERRQ(ierr) + call KSPSetConvergenceTest(ksp,BasicPETSC_convergedKSP,dummy,PETSC_NULL_FUNCTION,ierr); CHKERRQ(ierr) call SNESSetFromOptions(snes,ierr); CHKERRQ(ierr) !-------------------------------------------------------------------------------------------------- ! init fields - call DMDAVecGetArrayF90(da,solution_vec,xx_psc,ierr); CHKERRQ(ierr) ! get the data out of PETSc to work with + call DMDAVecGetArrayF90(da,solution_vec,xx_psc,ierr); CHKERRQ(ierr) ! get the data out of PETSc to work with F => xx_psc(0:8,:,:,:) - phaseField => xx_psc(9,:,:,:) - - if (restartInc == 1_pInt) then ! no deformation (no restart) + if (restartInc == 1_pInt) then ! no deformation (no restart) F_lastInc = spread(spread(spread(math_I3,3,grid(1)),4,grid(2)),5,grid(3)) ! initialize to identity - F = reshape(F_lastInc,[9,grid(1),grid(2),grid(3)]) + xx_psc(0:8,:,:,:) = reshape(F_lastInc,[9,grid(1),grid(2),grid(3)]) F_lastInc2 = F_lastInc - phaseField = temperature - phaseField_lastInc = temperature + do i = 1, nActivePhaseFields + xx_psc(8+i,:,:,:) = phaseFieldData(i)%phaseField0 + phaseField_lastInc(i,:,:,:) = phaseFieldData(i)%phaseField0 + enddo elseif (restartInc > 1_pInt) then ! using old values from file if (iand(debug_level(debug_spectral),debug_spectralRestart)/= 0) & write(6,'(/,a,'//IO_intOut(restartInc-1_pInt)//',a)') & @@ -257,7 +265,8 @@ end subroutine basicPETSc_init !> @brief solution for the Basic PETSC scheme with internal iterations !-------------------------------------------------------------------------------------------------- type(tSolutionState) function basicPETSc_solution( & - incInfoIn,guess,timeinc,timeinc_old,loadCaseTime,P_BC,F_BC,temperature_bc,rotation_BC,density) + incInfoIn,guess,timeinc,timeinc_old,loadCaseTime,P_BC,F_BC,temperature_bc,rotation_BC,density, & + nActivePhaseFields,phaseFieldData) use numerics, only: & update_gamma, & itmax @@ -290,6 +299,8 @@ type(tSolutionState) function basicPETSc_solution( & !-------------------------------------------------------------------------------------------------- ! input data for solution + integer(pInt), intent(in) :: nActivePhaseFields + type(phaseFieldDataBin), intent(in) :: phaseFieldData(nActivePhaseFields) real(pReal), intent(in) :: & timeinc, & !< increment in time for current solution timeinc_old, & !< increment in time of last increment @@ -304,17 +315,17 @@ type(tSolutionState) function basicPETSc_solution( & character(len=*), intent(in) :: & incInfoIn real(pReal), dimension(3,3), intent(in) :: rotation_BC + integer(pInt) :: i !-------------------------------------------------------------------------------------------------- ! PETSc Data - PetscScalar, pointer :: xx_psc(:,:,:,:), F(:,:,:,:), phaseField(:,:,:) + PetscScalar, pointer :: xx_psc(:,:,:,:), F(:,:,:,:) PetscErrorCode :: ierr SNESConvergedReason :: reason incInfo = incInfoIn call DMDAVecGetArrayF90(da,solution_vec,xx_psc,ierr); CHKERRQ(ierr) ! get the data out of PETSc to work with F => xx_psc(0:8,:,:,:) - phaseField => xx_psc(9,:,:,:) !-------------------------------------------------------------------------------------------------- ! restart information for spectral solver if (restartWrite) then @@ -343,9 +354,11 @@ type(tSolutionState) function basicPETSc_solution( & F,[3,3,grid(1),grid(2),grid(3)])),[3,1,product(grid)]) if (cutBack) then F_aim = F_aim_lastInc - F = reshape(F_lastInc,[9,grid(1),grid(2),grid(3)]) + xx_psc(0:8,:,:,:) = reshape(F_lastInc,[9,grid(1),grid(2),grid(3)]) C_volAvg = C_volAvgLastInc - phaseField = phaseField_lastInc + do i = 1, nActivePhaseFields + xx_psc(8+i,:,:,:) = phaseField_lastInc(i,:,:,:) + enddo else C_volAvgLastInc = C_volAvg @@ -368,17 +381,21 @@ type(tSolutionState) function basicPETSc_solution( & F,[3,3,grid(1),grid(2),grid(3)])),[3,1,product(grid)]) Fdot = Utilities_calculateRate(math_rotate_backward33(f_aimDot,params%rotation_BC), & timeinc_old,guess,F_lastInc,reshape(F,[3,3,grid(1),grid(2),grid(3)])) - phaseFieldDot = (phaseField - phaseField_lastInc)/timeinc_old + do i = 1, nActivePhaseFields + phaseFieldDot(i,:,:,:) = (xx_psc(8+i,:,:,:) - phaseField_lastInc(i,:,:,:))/timeinc_old + phaseField_lastInc(i,:,:,:) = xx_psc(8+i,:,:,:) + phaseFieldRHS_lastInc(i,:,:,:) = phaseFieldRHS(i,:,:,:) + enddo F_lastInc2 = F_lastInc - F_lastInc = reshape(F,[3,3,grid(1),grid(2),grid(3)]) - phaseField_lastInc = phaseField - heatSource_lastInc = reshape(materialpoint_heat(1,:),[grid(1),grid(2),grid(3)]) + F_lastInc = reshape(F,[3,3,grid(1),grid(2),grid(3)]) endif F_aim = F_aim + f_aimDot * timeinc F = reshape(Utilities_forwardField(timeinc,F_lastInc,Fdot,math_rotate_backward33(F_aim, & rotation_BC)),[9,grid(1),grid(2),grid(3)]) - phaseField = phaseField_lastInc + phaseFieldDot*timeinc + do i = 1, nActivePhaseFields + xx_psc(8+i,:,:,:) = phaseField_lastInc(i,:,:,:) + phaseFieldDot(i,:,:,:)*timeinc + enddo call DMDAVecRestoreArrayF90(da,solution_vec,xx_psc,ierr); CHKERRQ(ierr) !-------------------------------------------------------------------------------------------------- @@ -397,6 +414,8 @@ type(tSolutionState) function basicPETSc_solution( & params%timeincOld = timeinc_old params%temperature = temperature_BC params%density = density + params%nActivePhaseFields = nActivePhaseFields + params%phaseFieldData(1:nActivePhaseFields) = phaseFieldData(1:nActivePhaseFields) call SNESSolve(snes,PETSC_NULL_OBJECT,solution_vec,ierr); CHKERRQ(ierr) call SNESGetConvergedReason(snes,reason,ierr); CHKERRQ(ierr) @@ -451,34 +470,32 @@ subroutine BasicPETSC_formResidual(in,x_scal,f_scal,dummy,ierr) use crystallite, only: & crystallite_temperature use homogenization, only: & - materialpoint_heat + materialpoint_heat, & + materialpoint_P + use constitutive, only: & + constitutive_damage implicit none DMDALocalInfo, dimension(DMDA_LOCAL_INFO_SIZE) :: & in - PetscScalar, target, dimension(10, & + PetscScalar, target, dimension(9+params%nActivePhaseFields, & XG_RANGE,YG_RANGE,ZG_RANGE) :: & x_scal - PetscScalar, target, dimension(10, & + PetscScalar, target, dimension(9+params%nActivePhaseFields, & X_RANGE,Y_RANGE,Z_RANGE) :: & f_scal PetscScalar, pointer, dimension(:,:,:,:) :: & F, & residual_F - PetscScalar, pointer, dimension(:,:,:) :: & - phaseField, & - residual_phaseField PetscInt :: & PETScIter, & nfuncs PetscObject :: dummy PetscErrorCode :: ierr - real(pReal) :: mobility, diffusivity + integer(pInt) :: i F => x_scal(1:9,1:grid(1),1:grid(2),1:grid(3)) residual_F => f_scal(1:9,1:grid(1),1:grid(2),1:grid(3)) - phaseField => x_scal(10,1:grid(1),1:grid(2),1:grid(3)) - residual_phaseField => f_scal(10,1:grid(1),1:grid(2),1:grid(3)) call SNESGetNumberFunctionEvals(snes,nfuncs,ierr); CHKERRQ(ierr) call SNESGetIterationNumber(snes,PETScIter,ierr); CHKERRQ(ierr) @@ -517,7 +534,15 @@ subroutine BasicPETSC_formResidual(in,x_scal,f_scal,dummy,ierr) !-------------------------------------------------------------------------------------------------- ! evaluate constitutive response - crystallite_temperature(1,1_pInt:product(grid)) = reshape(phaseField,[product(grid)]) + do i = 1, params%nActivePhaseFields + if(params%phaseFieldData(i)%label == 'thermal') & + crystallite_temperature(1,1_pInt:product(grid)) = & + reshape(x_scal(9+i,1:grid(1),1:grid(2),1:grid(3)),[product(grid)]) + if(params%phaseFieldData(i)%label == 'fracture') & + constitutive_damage(1,1,1:product(grid)) = & + reshape(x_scal(9+i,1:grid(1),1:grid(2),1:grid(3)),[product(grid)]) + enddo + call Utilities_constitutiveResponse(F_lastInc,F,params%temperature,params%timeinc, & residual_F,C_volAvg,C_minmaxAvg,P_av,ForwardData,params%rotation_BC) ForwardData = .false. @@ -539,32 +564,78 @@ subroutine BasicPETSC_formResidual(in,x_scal,f_scal,dummy,ierr) call Utilities_fourierConvolution(math_rotate_backward33(F_aim_lastIter-F_aim,params%rotation_BC)) call Utilities_FFTbackward() +!-------------------------------------------------------------------------------------------------- +! constructing phase field residual + do i = 1, params%nActivePhaseFields + select case (params%phaseFieldData(i)%label) + case ('thermal') + phaseField_real = 0.0_pReal + phaseField_real(1:grid(1),1:grid(2),1:grid(3)) = & + phaseField_lastInc(i,1:grid(1),1:grid(2),1:grid(3)) + call utilities_scalarFFTforward() + call utilities_diffusion(params%phaseFieldData(i)%diffusion,params%timeinc) + call utilities_scalarFFTbackward() + f_scal(9+i,1:grid(1),1:grid(2),1:grid(3)) = & + phaseField_real(1:grid(1),1:grid(2),1:grid(3)) + + phaseFieldRHS(i,1:grid(1),1:grid(2),1:grid(3)) = & + reshape(materialpoint_heat(1,1_pInt:product(grid)),[grid(1),grid(2),grid(3)]) + phaseField_real = 0.0_pReal + phaseField_real(1:grid(1),1:grid(2),1:grid(3)) = & + params%timeinc*params%phaseFieldData(i)%mobility* & + (phaseFieldRHS_lastInc(i,1:grid(1),1:grid(2),1:grid(3)) + & + phaseFieldRHS (i,1:grid(1),1:grid(2),1:grid(3)))/2.0_pReal + call utilities_scalarFFTforward() + call utilities_diffusion(params%phaseFieldData(i)%diffusion,params%timeinc/2.0_pReal) + call utilities_scalarFFTbackward() + f_scal(9+i,1:grid(1),1:grid(2),1:grid(3)) = & + x_scal(9+i,1:grid(1),1:grid(2),1:grid(3)) - & + f_scal(9+i,1:grid(1),1:grid(2),1:grid(3)) - & + phaseField_real(1:grid(1),1:grid(2),1:grid(3)) + err_phaseField(i) = maxval(abs(f_scal(9+i,1:grid(1),1:grid(2),1:grid(3)))) + phaseField_Avg(i) = sum(x_scal(9+i,1:grid(1),1:grid(2),1:grid(3)))*wgt + + case ('fracture') + phaseField_real = 0.0_pReal + phaseField_real(1:grid(1),1:grid(2),1:grid(3)) = & + phaseField_lastInc(i,1:grid(1),1:grid(2),1:grid(3)) + call utilities_scalarFFTforward() + call utilities_diffusion(0.5_pReal*maxval(geomSize/real(grid,pReal))* & + params%phaseFieldData(i)%diffusion,params%timeinc) + call utilities_scalarFFTbackward() + f_scal(9+i,1:grid(1),1:grid(2),1:grid(3)) = & + phaseField_real(1:grid(1),1:grid(2),1:grid(3)) + + phaseFieldRHS(i,1:grid(1),1:grid(2),1:grid(3)) = & + - params%phaseFieldData(i)%mobility* & + sum(residual_F* & + (F-reshape(spread(spread(spread(math_I3,3,grid(1)),4,grid(2)),5,grid(3)),[9,grid(1),grid(2),grid(3)])),dim=1) & + - params%phaseFieldData(i)%diffusion*(x_scal(9+i,1:grid(1),1:grid(2),1:grid(3)) - 1.0_pReal)/ & + 2.0_pReal/maxval(geomSize/real(grid,pReal)) + phaseField_real = 0.0_pReal + phaseField_real(1:grid(1),1:grid(2),1:grid(3)) = & + params%timeinc*params%phaseFieldData(i)%mobility* & + (phaseFieldRHS_lastInc(i,1:grid(1),1:grid(2),1:grid(3)) + & + phaseFieldRHS (i,1:grid(1),1:grid(2),1:grid(3)))/2.0_pReal + call utilities_scalarFFTforward() + call utilities_diffusion(2.0_pReal*maxval(geomSize/real(grid,pReal))* & + params%phaseFieldData(i)%diffusion,params%timeinc/2.0_pReal) + call utilities_scalarFFTbackward() + f_scal(9+i,1:grid(1),1:grid(2),1:grid(3)) = & + x_scal(9+i,1:grid(1),1:grid(2),1:grid(3)) - & + f_scal(9+i,1:grid(1),1:grid(2),1:grid(3)) - & + phaseField_real(1:grid(1),1:grid(2),1:grid(3)) + err_phaseField(i) = maxval(abs(f_scal(9+i,1:grid(1),1:grid(2),1:grid(3)))) + phaseField_Avg(i) = sum(x_scal(9+i,1:grid(1),1:grid(2),1:grid(3)))*wgt + + end select + enddo + !-------------------------------------------------------------------------------------------------- ! constructing residual residual_F = reshape(field_real(1:grid(1),1:grid(2),1:grid(3),1:3,1:3),& [9,grid(1),grid(2),grid(3)],order=[2,3,4,1]) - -!-------------------------------------------------------------------------------------------------- -! constructing phase field residual - diffusivity = 400.0; mobility = 4e6 !sample coefficients for copper... need to clear this up with better mode of input - phaseField_real = 0.0_pReal - phaseField_real(1:grid(1),1:grid(2),1:grid(3)) = phaseField_lastInc - call utilities_scalarFFTforward() - call utilities_diffusion(diffusivity/mobility,params%timeinc) - call utilities_scalarFFTbackward() - residual_phaseField = phaseField_real(1:grid(1),1:grid(2),1:grid(3)) - phaseField_real = 0.0_pReal - phaseField_real(1:grid(1),1:grid(2),1:grid(3)) = & - (heatSource_lastInc + reshape(materialpoint_heat(1,1_pInt:product(grid)),[grid(1),grid(2),grid(3)]))* & - params%timeinc/mobility/2.0_pReal - call utilities_scalarFFTforward() - call utilities_diffusion(diffusivity/mobility,params%timeinc/2.0_pReal) - call utilities_scalarFFTbackward() - residual_phaseField = phaseField - & - (residual_phaseField + phaseField_real(1:grid(1),1:grid(2),1:grid(3))) - - err_phaseField = maxval(abs(residual_phaseField)); phaseField_Avg = sum(phaseField)*wgt - + end subroutine BasicPETSc_formResidual @@ -601,9 +672,8 @@ subroutine BasicPETSc_converged(snes_local,PETScIter,xnorm,snorm,fnorm,reason,du err_divPrev = err_div; err_div = err_divDummy converged: if ((totalIter >= itmin .and. & - all([ err_div/divTol, & - err_stress/stressTol, & - err_phaseField/phaseField_Avg/1.0e-3 ] < 1.0_pReal)) & + all([ err_div/divTol, err_stress/stressTol] < 1.0_pReal) .and. & + maxval(err_phaseField/phaseField_Avg) < 1.0e-3_pReal) & .or. terminallyIll) then reason = 1 elseif (totalIter >= itmax) then converged @@ -620,7 +690,7 @@ subroutine BasicPETSc_converged(snes_local,PETScIter,xnorm,snorm,fnorm,reason,du write(6,'(a,f12.2,a,es8.2,a,es9.2,a)') ' error stress BC = ', & err_stress/stressTol, ' (',err_stress, ' Pa, tol =',stressTol,')' write(6,'(a,f10.2,a,es8.2,a,es9.2,a)') ' error phase field = ', & - err_phaseField/phaseField_Avg/1.0e-3, ' (',err_phaseField/phaseField_Avg, ' Pa, tol =',1.0e-3,')' + maxval(err_phaseField/phaseField_Avg)/1.0e-3, ' (',maxval(err_phaseField/phaseField_Avg), ' Pa, tol =',1.0e-3,')' write(6,'(/,a)') ' ===========================================================================' flush(6) diff --git a/code/DAMASK_spectral_utilities.f90 b/code/DAMASK_spectral_utilities.f90 index 31ecbc59a..5ce3e6ec4 100644 --- a/code/DAMASK_spectral_utilities.f90 +++ b/code/DAMASK_spectral_utilities.f90 @@ -36,6 +36,7 @@ module DAMASK_spectral_utilities #include #endif logical, public :: cutBack =.false. !< cut back of BVP solver in case convergence is not achieved or a material point is terminally ill + integer(pInt), public, parameter :: maxPhaseFields = 2_pInt !-------------------------------------------------------------------------------------------------- ! grid related information information integer(pInt), public, dimension(3) :: grid !< grid points as specified in geometry file @@ -99,6 +100,14 @@ module DAMASK_spectral_utilities logical, dimension(3,3) :: maskLogical = .false. character(len=64) :: myType = 'None' end type tBoundaryCondition + + type, public :: phaseFieldDataBin !< set of parameters defining a phase field + real(pReal) :: diffusion = 0.0_pReal, & !< thermal conductivity + mobility = 0.0_pReal, & !< thermal mobility + phaseField0 = 0.0_pReal !< homogeneous damage field starting condition + logical :: active = .false. + character(len=64) :: label = '' + end type phaseFieldDataBin public :: & utilities_init, & @@ -116,8 +125,7 @@ module DAMASK_spectral_utilities utilities_constitutiveResponse, & utilities_calculateRate, & utilities_forwardField, & - utilities_destroy, & - utilities_temperatureUpdate + utilities_destroy private :: & utilities_getFilter @@ -612,9 +620,6 @@ subroutine utilities_diffusion(coefficient,timeinc) real(pReal),intent(in) :: timeinc, coefficient integer(pInt) :: i, j, k integer(pInt), dimension(3) :: k_s - - write(6,'(/,a)') ' ... doing diffusion .......................................................' - flush(6) !-------------------------------------------------------------------------------------------------- ! do the actual spectral method calculation (mechanical equilibrium) @@ -1116,66 +1121,6 @@ real(pReal) function utilities_getFilter(k) end function utilities_getFilter -!-------------------------------------------------------------------------------------------------- -!> @brief calculates filter for fourier convolution depending on type given in numerics.config -!-------------------------------------------------------------------------------------------------- -subroutine utilities_temperatureUpdate(timeinc) - use crystallite, only: & - crystallite_temperature - use homogenization, only: & - materialpoint_heat - - implicit none - real(pReal),intent(in) :: timeinc - integer :: & - x,y,z,e - real(pReal) :: & - a - forall(e=1_pInt:product(grid)) & - crystallite_temperature(1,e) = crystallite_temperature(1,e) + materialpoint_heat(1,e)*timeinc - e = 0_pInt - z = 0_pInt - y = 0_pInt - x = 0_pInt - -!< 6 or less neighboring IPs as [element_num, IP_index, neighbor_index that points to me] - do z = 0_pInt,grid(3)-1_pInt - do y = 0_pInt,grid(2)-1_pInt - do x = 0_pInt,grid(1)-1_pInt - e = e + 1_pInt - a = 0.0_pReal - a = a+ crystallite_temperature(1, z * grid(1) * grid(2) & - + y * grid(1) & - + modulo(x+1_pInt,grid(1)) & - + 1_pInt) - a = a+ crystallite_temperature(1,z * grid(1) * grid(2) & - + y * grid(1) & - + modulo(x-1_pInt,grid(1)) & - + 1_pInt) - a = a+ crystallite_temperature(1,z * grid(1) * grid(2) & - + modulo(y+1_pInt,grid(2)) * grid(1) & - + x & - + 1_pInt) - a = a+ crystallite_temperature(1,z * grid(1) * grid(2) & - + modulo(y-1_pInt,grid(2)) * grid(1) & - + x & - + 1_pInt) - a = a+ crystallite_temperature(1, modulo(z+1_pInt,grid(3)) * grid(1) * grid(2) & - + y * grid(1) & - + x & - + 1_pInt) - a = a+ crystallite_temperature(1,modulo(z-1_pInt,grid(3)) * grid(1) * grid(2) & - + y * grid(1) & - + x & - + 1_pInt) - crystallite_temperature(1,e) = (crystallite_temperature(1,e)+a)/7.0_pReal - enddo - enddo - enddo - -end subroutine utilities_temperatureUpdate - - !-------------------------------------------------------------------------------------------------- !> @brief cleans up !--------------------------------------------------------------------------------------------------