diff --git a/PRIVATE b/PRIVATE index 29ef436ac..62df7f24f 160000 --- a/PRIVATE +++ b/PRIVATE @@ -1 +1 @@ -Subproject commit 29ef436acca5417aebc945b688642c34697af911 +Subproject commit 62df7f24f2a95fda255f7d20b130afcfeecb1b4a diff --git a/examples/mesh/tensionY_mono.load b/examples/mesh/tensionY_mono.load deleted file mode 100644 index e7ab69fc4..000000000 --- a/examples/mesh/tensionY_mono.load +++ /dev/null @@ -1,14 +0,0 @@ -# initial elastic step -$Loadcase 1 t 0.0005 N 1 f_out 1 - Face 3 Y -0.025 - Face 4 X 0.0 - Face 4 Y 0.0 - Face 4 Z 0.0 -$EndLoadcase -# plastic step -$Loadcase 2 t 1.0 N 10 f_out 2 - Face 3 Y -0.025 - Face 4 X 0.0 - Face 4 Y 0.0 - Face 4 Z 0.0 -$EndLoadcase diff --git a/examples/mesh/tensionY_mono.yaml b/examples/mesh/tensionY_mono.yaml new file mode 100644 index 000000000..fb7da1b86 --- /dev/null +++ b/examples/mesh/tensionY_mono.yaml @@ -0,0 +1,22 @@ +--- +loadstep: + - boundary_conditions: + mechanical: + - dot_u: ['x', -0.025, 'x'] + tag: 3 + - dot_u: [0.0, 0.0, 0.0] + tag: 4 + discretization: + t: 0.0005 + N: 1 + f_out: 1 + - boundary_conditions: + mechanical: + - dot_u: ['x', -0.025, 'x'] + tag: 3 + - dot_u: [0.0, 0.0, 0.0] + tag: 4 + discretization: + t: 1.0 + N: 10 + f_out: 2 diff --git a/examples/mesh/tensionZ_3g.load b/examples/mesh/tensionZ_3g.load deleted file mode 100644 index b888873ea..000000000 --- a/examples/mesh/tensionZ_3g.load +++ /dev/null @@ -1,14 +0,0 @@ -# initial elastic step -$Loadcase 1 t 0.0005 N 1 f_out 1 - Face 1 Z 0.01 - Face 2 X 0.0 - Face 2 Y 0.0 - Face 2 Z 0.0 -$EndLoadcase -# plastic step -$Loadcase 2 t 1.0 N 10 f_out 2 - Face 1 Z 0.01 - Face 2 X 0.0 - Face 2 Y 0.0 - Face 2 Z 0.0 -$EndLoadcase diff --git a/examples/mesh/tensionZ_3g.yaml b/examples/mesh/tensionZ_3g.yaml new file mode 100644 index 000000000..e10f80229 --- /dev/null +++ b/examples/mesh/tensionZ_3g.yaml @@ -0,0 +1,22 @@ +--- +loadstep: + - boundary_conditions: + mechanical: + - dot_u: ['x', 'x', 0.01] + tag: 1 + - dot_u: [0.0, 0.0, 0.0] + tag: 2 + discretization: + t: 0.0005 + N: 1 + f_out: 1 + - boundary_conditions: + mechanical: + - dot_u: ['x', 'x', 0.01] + tag: 1 + - dot_u: [0.0, 0.0, 0.0] + tag: 2 + discretization: + t: 1.0 + N: 10 + f_out: 2 diff --git a/examples/mesh/tensionZ_5g.load b/examples/mesh/tensionZ_5g.load deleted file mode 100644 index 475931523..000000000 --- a/examples/mesh/tensionZ_5g.load +++ /dev/null @@ -1,18 +0,0 @@ -# initial elastic step -$Loadcase 1 t 0.0005 N 1 f_out 1 - Face 1 X 0.0 - Face 1 Y 0.0 - Face 1 Z 0.0 - Face 2 X 0.0 - Face 2 Y 0.0 - Face 2 Z 0.0025 -$EndLoadcase -# plastic step -$Loadcase 2 t 1.0 N 10 f_out 2 - Face 1 X 0.0 - Face 1 Y 0.0 - Face 1 Z 0.0 - Face 2 X 0.0 - Face 2 Y 0.0 - Face 2 Z 0.0025 -$EndLoadcase diff --git a/examples/mesh/tensionZ_5g.yaml b/examples/mesh/tensionZ_5g.yaml new file mode 100644 index 000000000..f43e1e326 --- /dev/null +++ b/examples/mesh/tensionZ_5g.yaml @@ -0,0 +1,22 @@ +--- +loadstep: + - boundary_conditions: + mechanical: + - dot_u: [0.0, 0.0, 0.0] + tag: 1 + - dot_u: [0.0, 0.0, 0.0025] + tag: 2 + discretization: + t: 0.0005 + N: 1 + f_out: 1 + - boundary_conditions: + mechanical: + - dot_u: [0.0, 0.0, 0.0] + tag: 1 + - dot_u: [0.0, 0.0, 0.0025] + tag: 2 + discretization: + t: 1.0 + N: 10 + f_out: 2 diff --git a/src/mesh/DAMASK_mesh.f90 b/src/mesh/DAMASK_mesh.f90 index 370dc3535..b8eb64d7d 100644 --- a/src/mesh/DAMASK_mesh.f90 +++ b/src/mesh/DAMASK_mesh.f90 @@ -23,38 +23,32 @@ program DAMASK_mesh implicit none(type,external) type :: tLoadCase - real(pREAL) :: time = 0.0_pREAL !< length of increment - integer :: incs = 0, & !< number of increments - outputfrequency = 1 !< frequency of result writes - logical :: followFormerTrajectory = .true. !< follow trajectory of former loadcase - integer, allocatable, dimension(:) :: faceID - type(tFieldBC), allocatable, dimension(:) :: fieldBC + real(pREAL) :: t = 0.0_pREAL !< length of increment + integer :: N = 0, & !< number of increments + f_out = 1 !< frequency of result writes + logical :: estimate_rate = .true. !< follow trajectory of former loadcase + integer, allocatable, dimension(:) :: tag + type(tMechBC) :: mechBC end type tLoadCase -!-------------------------------------------------------------------------------------------------- -! variables related to information from load case and geom file - integer, allocatable, dimension(:) :: chunkPos ! this is longer than needed for geometry parsing - integer :: & - N_def = 0 !< # of rate of deformation specifiers found in load case file - character(len=:), allocatable :: & - line !-------------------------------------------------------------------------------------------------- ! loop variables, convergence etc. integer, parameter :: & subStepFactor = 2 !< for each substep, divide the last time increment by 2.0 real(pREAL) :: & - time = 0.0_pREAL, & !< elapsed time - time0 = 0.0_pREAL, & !< begin of interval - timeinc = 0.0_pREAL, & !< current time interval - timeIncOld = 0.0_pREAL, & !< previous time interval - remainingLoadCaseTime = 0.0_pREAL !< remaining time of current load case + t = 0.0_pREAL, & !< elapsed time + t_0 = 0.0_pREAL, & !< begin of interval + Delta_t = 0.0_pREAL, & !< current time interval + Delta_t_prev = 0.0_pREAL, & !< previous time interval + t_remaining = 0.0_pREAL !< remaining time of current load case logical :: & guess, & !< guess along former trajectory stagIterate integer :: & l, & i, & + m, & errorID, & cutBackLevel = 0, & !< cut back level \f$ t = \frac{t_{inc}}{2^l} \f$ stepFraction = 0, & !< fraction of current time interval @@ -67,8 +61,16 @@ program DAMASK_mesh component type(tDict), pointer :: & num_solver, & - num_mesh - character(len=pSTRLEN), dimension(:), allocatable :: fileContent + num_mesh, & + load, & + load_step, & + step_bc, & + mech_BC, & + step_discretization + type(tList), pointer :: & + load_steps, & + mech_u, & + step_mech character(len=pSTRLEN) :: & incInfo, & loadcase_string @@ -83,6 +85,9 @@ program DAMASK_mesh integer(kind(COMPONENT_UNDEFINED_ID)) :: ID external :: & quit + character(len=:), allocatable :: & + fileContent, fname + !-------------------------------------------------------------------------------------------------- ! init DAMASK (all modules) @@ -104,135 +109,101 @@ program DAMASK_mesh CHKERRA(err_PETSc) allocate(solres(1)) -!-------------------------------------------------------------------------------------------------- -! reading basic information from load case file and allocate data structure containing load cases - fileContent = IO_readlines(trim(CLI_loadFile)) - do l = 1, size(fileContent) - line = fileContent(l) - if (IO_isBlank(line)) cycle ! skip empty lines + if (worldrank == 0) then + fileContent = IO_read(CLI_loadFile) + fname = CLI_loadFile + if (scan(fname,'/') /= 0) fname = fname(scan(fname,'/',.true.)+1:) + call result_openJobFile(parallel=.false.) + call result_addSetupFile(fileContent,fname,'load case definition (mesh solver)') + call result_closeJobFile() + end if - chunkPos = IO_strPos(line) - do i = 1, chunkPos(1) ! reading compulsory parameters for loadcase - select case (IO_strValue(line,chunkPos,i)) - case('$Loadcase') - N_def = N_def + 1 - end select - end do ! count all identifiers to allocate memory and do sanity check - end do + call parallelization_bcast_str(fileContent) + load => YAML_parse_str_asDict(fileContent) + load_steps => load%get_list('loadstep') - if (N_def < 1) call IO_error(error_ID = 837) - allocate(loadCases(N_def)) + allocate(loadCases(load_steps%length)) - do i = 1, size(loadCases) - allocate(loadCases(i)%fieldBC(1)) - loadCases(i)%fieldBC(1)%ID = FIELD_MECH_ID - end do - - do i = 1, size(loadCases) - loadCases(i)%fieldBC(1)%nComponents = dimPlex !< X, Y (, Z) displacements - allocate(loadCases(i)%fieldBC(1)%componentBC(loadCases(i)%fieldBC(1)%nComponents)) - do component = 1, loadCases(i)%fieldBC(1)%nComponents + do l = 1, load_steps%length + load_step => load_steps%get_dict(l) + step_bc => load_step%get_dict('boundary_conditions') + step_mech => step_bc%get_list('mechanical') + loadCases(l)%mechBC%nComponents = dimPlex !< X, Y (, Z) displacements + allocate(loadCases(l)%mechBC%componentBC(dimPlex)) + do component = 1, dimPlex select case (component) case (1) - loadCases(i)%fieldBC(1)%componentBC(component)%ID = COMPONENT_MECH_X_ID + loadCases(l)%mechBC%componentBC(component)%ID = COMPONENT_MECH_X_ID case (2) - loadCases(i)%fieldBC(1)%componentBC(component)%ID = COMPONENT_MECH_Y_ID + loadCases(l)%mechBC%componentBC(component)%ID = COMPONENT_MECH_Y_ID case (3) - loadCases(i)%fieldBC(1)%componentBC(component)%ID = COMPONENT_MECH_Z_ID + loadCases(l)%mechBC%componentBC(component)%ID = COMPONENT_MECH_Z_ID end select end do - do component = 1, loadCases(i)%fieldBC(1)%nComponents - allocate(loadCases(i)%fieldBC(1)%componentBC(component)%Value(mesh_Nboundaries), source = 0.0_pREAL) - allocate(loadCases(i)%fieldBC(1)%componentBC(component)%Mask (mesh_Nboundaries), source = .false.) + do component = 1, dimPlex + allocate(loadCases(l)%mechBC%componentBC(component)%Value(mesh_Nboundaries), source = 0.0_pREAL) + allocate(loadCases(l)%mechBC%componentBC(component)%Mask (mesh_Nboundaries), source = .false.) end do - end do -!-------------------------------------------------------------------------------------------------- -! reading the load case and assign values to the allocated data structure - do l = 1, size(fileContent) - line = fileContent(l) - if (IO_isBlank(line)) cycle ! skip empty lines - - chunkPos = IO_strPos(line) - do i = 1, chunkPos(1) - select case (IO_strValue(line,chunkPos,i)) -!-------------------------------------------------------------------------------------------------- -! loadcase information - case('$Loadcase') - currentLoadCase = IO_intValue(line,chunkPos,i+1) - case('Face') - currentFace = IO_intValue(line,chunkPos,i+1) - currentFaceSet = -1 - do faceSet = 1, mesh_Nboundaries - if (mesh_boundaries(faceSet) == currentFace) currentFaceSet = faceSet - end do - if (currentFaceSet < 0) call IO_error(error_ID = 837, ext_msg = 'invalid BC') - case('t') - loadCases(currentLoadCase)%time = IO_realValue(line,chunkPos,i+1) - case('N') - loadCases(currentLoadCase)%incs = IO_intValue(line,chunkPos,i+1) - case('f_out') - loadCases(currentLoadCase)%outputfrequency = IO_intValue(line,chunkPos,i+1) - case('estimate_rate') - loadCases(currentLoadCase)%followFormerTrajectory = .false. ! do not continue to predict deformation along former trajectory - -!-------------------------------------------------------------------------------------------------- -! boundary condition information - case('X','Y','Z') - select case(IO_strValue(line,chunkPos,i)) - case('X') - ID = COMPONENT_MECH_X_ID - case('Y') - ID = COMPONENT_MECH_Y_ID - case('Z') - ID = COMPONENT_MECH_Z_ID - end select - - do component = 1, loadcases(currentLoadCase)%fieldBC(1)%nComponents - if (loadCases(currentLoadCase)%fieldBC(1)%componentBC(component)%ID == ID) then - loadCases(currentLoadCase)%fieldBC(1)%componentBC(component)%Mask (currentFaceSet) = & - .true. - loadCases(currentLoadCase)%fieldBC(1)%componentBC(component)%Value(currentFaceSet) = & - IO_realValue(line,chunkPos,i+1) - end if - end do - end select + do m = 1, step_mech%length + mech_BC => step_mech%get_dict(m) + currentFaceSet = -1 + do faceSet = 1, mesh_Nboundaries + if (mesh_boundaries(faceSet) == mech_BC%get_asInt('tag')) currentFaceSet = faceSet + end do + if (currentFaceSet < 0) call IO_error(error_ID = 837, ext_msg = 'invalid BC') + do component = 1, dimPlex + mech_u => mech_BC%get_list('dot_u') + if (mech_u%get_asStr(component) /= 'x') then + loadCases(l)%mechBC%componentBC(component)%Mask(currentFaceSet) = .true. + loadCases(l)%mechBC%componentBC(component)%Value(currentFaceSet) = mech_u%get_asReal(component) + end if + end do end do + step_discretization => load_step%get_dict('discretization') + loadCases(l)%t = step_discretization%get_asReal('t') + loadCases(l)%N = step_discretization%get_asInt ('N') + + if (load_step%get_asStr('f_out',defaultVal='n/a') == 'none') then + loadCases(l)%f_out = huge(0) + else + loadCases(l)%f_out = load_step%get_asInt('f_out', defaultVal=1) + end if + loadCases(l)%estimate_rate = (load_step%get_asBool('estimate_rate',defaultVal=.true.) .and. l>1) end do !-------------------------------------------------------------------------------------------------- ! consistency checks and output of load case - loadCases(1)%followFormerTrajectory = .false. ! cannot guess along trajectory for first inc of first currentLoadCase errorID = 0 - checkLoadcases: do currentLoadCase = 1, size(loadCases) - write (loadcase_string, '(i0)' ) currentLoadCase - print'(/,1x,a,1x,i0)', 'load case:', currentLoadCase - if (.not. loadCases(currentLoadCase)%followFormerTrajectory) & + checkLoadcases: do l = 1, load_steps%length + write (loadcase_string, '(i0)' ) l + print'(/,1x,a,1x,i0)', 'load case:', l + if (.not. loadCases(l)%estimate_rate) & print'(2x,a)', 'drop guessing along trajectory' print'(2x,a)', 'Field '//trim(FIELD_MECH_label) do faceSet = 1, mesh_Nboundaries - do component = 1, loadCases(currentLoadCase)%fieldBC(1)%nComponents - if (loadCases(currentLoadCase)%fieldBC(1)%componentBC(component)%Mask(faceSet)) & + do component = 1, loadCases(l)%mechBC%nComponents + if (loadCases(l)%mechBC%componentBC(component)%Mask(faceSet)) & print'(a,i2,a,i2,a,f12.7)', & ' Face ', mesh_boundaries(faceSet), & ' Component ', component, & - ' Value ', loadCases(currentLoadCase)%fieldBC(1)%componentBC(component)%Value(faceSet) + ' Value ', loadCases(l)%mechBC%componentBC(component)%Value(faceSet) end do end do - print'(2x,a,f12.6)', 'time: ', loadCases(currentLoadCase)%time - if (loadCases(currentLoadCase)%incs < 1) errorID = 835 ! non-positive incs count - print'(2x,a,i5)', 'increments: ', loadCases(currentLoadCase)%incs - if (loadCases(currentLoadCase)%outputfrequency < 1) errorID = 836 ! non-positive result frequency + print'(2x,a,f12.6)', 'time: ', loadCases(l)%t + if (loadCases(l)%N < 1) errorID = 835 ! non-positive incs count + print'(2x,a,i5)', 'increments: ', loadCases(l)%N + if (loadCases(l)%f_out < 1) errorID = 836 ! non-positive result frequency print'(2x,a,i5)', 'output frequency: ', & - loadCases(currentLoadCase)%outputfrequency - if (errorID > 0) call IO_error(error_ID = errorID, ext_msg = loadcase_string) ! exit with error message + loadCases(l)%f_out + if (errorID > 0) call IO_error(error_ID = errorID, ext_msg = loadcase_string) ! exit with error message end do checkLoadcases !-------------------------------------------------------------------------------------------------- ! doing initialization depending on active solvers call FEM_Utilities_init(num_mesh) - call FEM_mechanical_init(loadCases(1)%fieldBC(1),num_mesh) + call FEM_mechanical_init(loadCases(1)%mechBC,num_mesh) call config_numerics_deallocate() if (worldrank == 0) then @@ -244,46 +215,46 @@ program DAMASK_mesh flush(IO_STDOUT) call materialpoint_result(0,0.0_pREAL) - loadCaseLooping: do currentLoadCase = 1, size(loadCases) - time0 = time ! load case start time - guess = loadCases(currentLoadCase)%followFormerTrajectory ! change of load case? homogeneous guess for the first inc + loadCaseLooping: do l = 1, load_steps%length + t_0 = t ! load case start time + guess = loadCases(l)%estimate_rate ! change of load case? homogeneous guess for the first inc - incLooping: do inc = 1, loadCases(currentLoadCase)%incs + incLooping: do inc = 1, loadCases(l)%N totalIncsCounter = totalIncsCounter + 1 !-------------------------------------------------------------------------------------------------- ! forwarding time - timeIncOld = timeinc ! last timeinc that brought former inc to an end - timeinc = loadCases(currentLoadCase)%time/real(loadCases(currentLoadCase)%incs,pREAL) - timeinc = timeinc * real(subStepFactor,pREAL)**real(-cutBackLevel,pREAL) ! depending on cut back level, decrease time step + Delta_t_prev = Delta_t ! last timeinc that brought former inc to an end + Delta_t = loadCases(l)%t/real(loadCases(l)%N,pREAL) + Delta_t = Delta_t * real(subStepFactor,pREAL)**real(-cutBackLevel,pREAL) ! depending on cut back level, decrease time step stepFraction = 0 ! fraction scaled by stepFactor**cutLevel subStepLooping: do while (stepFraction < subStepFactor**cutBackLevel) - remainingLoadCaseTime = loadCases(currentLoadCase)%time+time0 - time - time = time + timeinc ! forward target time + t_remaining = loadCases(l)%t + t_0 - t + t = t + Delta_t ! forward target time stepFraction = stepFraction + 1 ! count step !-------------------------------------------------------------------------------------------------- ! report begin of new step print'(/,1x,a)', '###########################################################################' print'(1x,a,es12.5,6(a,i0))',& - 'Time', time, & - 's: Increment ', inc, '/', loadCases(currentLoadCase)%incs,& + 'Time', t, & + 's: Increment ', inc, '/', loadCases(l)%N,& '-', stepFraction, '/', subStepFactor**cutBackLevel,& - ' of load case ', currentLoadCase,'/',size(loadCases) + ' of load case ', l,'/',load_steps%length write(incInfo,'(4(a,i0))') & - 'Increment ',totalIncsCounter,'/',sum(loadCases%incs),& + 'Increment ',totalIncsCounter,'/',sum(loadCases%N),& '-',stepFraction, '/', subStepFactor**cutBackLevel flush(IO_STDOUT) - call FEM_mechanical_forward(guess,timeinc,timeIncOld,loadCases(currentLoadCase)%fieldBC(1)) + call FEM_mechanical_forward(guess,Delta_t,Delta_t_prev,loadCases(l)%mechBC) !-------------------------------------------------------------------------------------------------- ! solve fields stagIter = 0 stagIterate = .true. do while (stagIterate) - solres(1) = FEM_mechanical_solution(incInfo,timeinc,timeIncOld,loadCases(currentLoadCase)%fieldBC(1)) + solres(1) = FEM_mechanical_solution(incInfo,Delta_t,Delta_t_prev,loadCases(l)%mechBC) if (.not. solres(1)%converged) exit stagIter = stagIter + 1 @@ -294,13 +265,13 @@ program DAMASK_mesh ! check solution cutBack = .False. - if (.not. all(solres(:)%converged .and. solres(:)%stagConverged)) then ! no solution found + if (.not. all(solres(:)%converged .and. solres(:)%stagConverged)) then ! no solution found if (cutBackLevel < maxCutBack) then ! do cut back cutBack = .True. stepFraction = (stepFraction - 1) * subStepFactor ! adjust to new denominator cutBackLevel = cutBackLevel + 1 - time = time - timeinc ! rewind time - timeinc = timeinc/2.0_pREAL + t = t - Delta_t ! rewind time + Delta_t = Delta_t/2.0_pREAL print'(/,1x,a)', 'cutting back' else ! default behavior, exit if spectral solver does not converge if (worldrank == 0) close(statUnit) @@ -308,10 +279,10 @@ program DAMASK_mesh end if else guess = .true. ! start guessing after first converged (sub)inc - timeIncOld = timeinc + Delta_t_prev = Delta_t end if if (.not. cutBack .and. worldrank == 0) then - write(statUnit,*) totalIncsCounter, time, cutBackLevel, & + write(statUnit,*) totalIncsCounter, t, cutBackLevel, & solres%converged, solres%iterationsNeeded ! write statistics about accepted solution flush(statUnit) end if @@ -325,10 +296,10 @@ program DAMASK_mesh print'(/,1x,a,1x,i0,1x,a)', 'increment', totalIncsCounter, 'NOT converged' end if; flush(IO_STDOUT) - if (mod(inc,loadCases(currentLoadCase)%outputFrequency) == 0) then ! at output frequency + if (mod(inc,loadCases(l)%f_out) == 0) then ! at output frequency print'(/,1x,a)', '... saving results ........................................................' call FEM_mechanical_updateCoords() - call materialpoint_result(totalIncsCounter,time) + call materialpoint_result(totalIncsCounter,t) end if diff --git a/src/mesh/FEM_utilities.f90 b/src/mesh/FEM_utilities.f90 index b1c218172..15a2168a8 100644 --- a/src/mesh/FEM_utilities.f90 +++ b/src/mesh/FEM_utilities.f90 @@ -38,10 +38,6 @@ module FEM_utilities character(len=*), parameter, public :: & FIELD_MECH_label = 'mechanical' - enum, bind(c); enumerator :: & - FIELD_UNDEFINED_ID, & - FIELD_MECH_ID - end enum enum, bind(c); enumerator :: & COMPONENT_UNDEFINED_ID, & COMPONENT_MECH_X_ID, & @@ -64,11 +60,10 @@ module FEM_utilities logical, allocatable, dimension(:) :: Mask end type tComponentBC - type, public :: tFieldBC - integer(kind(FIELD_UNDEFINED_ID)) :: ID + type, public :: tMechBC integer :: nComponents = 0 type(tComponentBC), allocatable, dimension(:) :: componentBC - end type tFieldBC + end type tMechBC external :: & ! ToDo: write interfaces PetscSectionGetFieldComponents, & @@ -79,7 +74,6 @@ module FEM_utilities FEM_utilities_init, & utilities_constitutiveResponse, & utilities_projectBCValues, & - FIELD_MECH_ID, & COMPONENT_UNDEFINED_ID, & COMPONENT_MECH_X_ID, & COMPONENT_MECH_Y_ID, & @@ -142,9 +136,9 @@ end subroutine FEM_utilities_init !-------------------------------------------------------------------------------------------------- !> @brief calculates constitutive response !-------------------------------------------------------------------------------------------------- -subroutine utilities_constitutiveResponse(timeinc,P_av,forwardData) +subroutine utilities_constitutiveResponse(Delta_t,P_av,forwardData) - real(pREAL), intent(in) :: timeinc !< loading time + real(pREAL), intent(in) :: Delta_t !< loading time logical, intent(in) :: forwardData !< age results real(pREAL),intent(out), dimension(3,3) :: P_av !< average PK stress @@ -152,9 +146,9 @@ subroutine utilities_constitutiveResponse(timeinc,P_av,forwardData) print'(/,1x,a)', '... evaluating constitutive response ......................................' - call homogenization_mechanical_response(timeinc,1,mesh_maxNips*mesh_NcpElems) ! calculate P field + call homogenization_mechanical_response(Delta_t,1,mesh_maxNips*mesh_NcpElems) ! calculate P field if (.not. terminallyIll) & - call homogenization_mechanical_response2(timeinc,[1,mesh_maxNips],[1,mesh_NcpElems]) + call homogenization_mechanical_response2(Delta_t,[1,mesh_maxNips],[1,mesh_NcpElems]) cutBack = .false. P_av = sum(homogenization_P,dim=3) * wgt @@ -168,7 +162,7 @@ end subroutine utilities_constitutiveResponse !-------------------------------------------------------------------------------------------------- !> @brief Project BC values to local vector !-------------------------------------------------------------------------------------------------- -subroutine utilities_projectBCValues(localVec,section,field,comp,bcPointsIS,BCValue,BCDotValue,timeinc) +subroutine utilities_projectBCValues(localVec,section,field,comp,bcPointsIS,BCValue,BCDotValue,Delta_t) Vec :: localVec PetscInt :: field, comp, nBcPoints, point, dof, numDof, numComp, offset @@ -176,7 +170,7 @@ subroutine utilities_projectBCValues(localVec,section,field,comp,bcPointsIS,BCVa IS :: bcPointsIS PetscInt, pointer :: bcPoints(:) real(pREAL), pointer :: localArray(:) - real(pREAL) :: BCValue,BCDotValue,timeinc + real(pREAL) :: BCValue,BCDotValue,Delta_t PetscErrorCode :: err_PETSc @@ -193,7 +187,7 @@ subroutine utilities_projectBCValues(localVec,section,field,comp,bcPointsIS,BCVa call PetscSectionGetFieldOffset(section,bcPoints(point),field,offset,err_PETSc) CHKERRQ(err_PETSc) do dof = offset+comp+1, offset+numDof, numComp - localArray(dof) = localArray(dof) + BCValue + BCDotValue*timeinc + localArray(dof) = localArray(dof) + BCValue + BCDotValue*Delta_t end do end do call VecRestoreArrayF90(localVec,localArray,err_PETSc) diff --git a/src/mesh/mesh_mech_FEM.f90 b/src/mesh/mesh_mech_FEM.f90 index 2d5556e63..824ebecdd 100644 --- a/src/mesh/mesh_mech_FEM.f90 +++ b/src/mesh/mesh_mech_FEM.f90 @@ -36,8 +36,8 @@ module mesh_mechanical_FEM !-------------------------------------------------------------------------------------------------- ! derived types type tSolutionParams - type(tFieldBC) :: fieldBC - real(pREAL) :: timeinc + type(tMechBC) :: mechBC + real(pREAL) :: Delta_t end type tSolutionParams type(tSolutionParams) :: params @@ -97,9 +97,9 @@ contains !-------------------------------------------------------------------------------------------------- !> @brief allocates all neccessary fields and fills them with data !-------------------------------------------------------------------------------------------------- -subroutine FEM_mechanical_init(fieldBC,num_mesh) +subroutine FEM_mechanical_init(mechBC,num_mesh) - type(tFieldBC), intent(in) :: fieldBC + type(tMechBC), intent(in) :: mechBC type(tDict), pointer, intent(in) :: num_mesh DM :: mechanical_mesh @@ -209,14 +209,14 @@ subroutine FEM_mechanical_init(fieldBC,num_mesh) end do numBC = 0 do field = 1, dimPlex; do faceSet = 1, mesh_Nboundaries - if (fieldBC%componentBC(field)%Mask(faceSet)) numBC = numBC + 1 + if (mechBC%componentBC(field)%Mask(faceSet)) numBC = numBC + 1 end do; end do allocate(pbcField(numBC), source=0_pPETSCINT) allocate(pbcComps(numBC)) allocate(pbcPoints(numBC)) numBC = 0 do field = 1, dimPlex; do faceSet = 1, mesh_Nboundaries - if (fieldBC%componentBC(field)%Mask(faceSet)) then + if (mechBC%componentBC(field)%Mask(faceSet)) then numBC = numBC + 1 call ISCreateGeneral(PETSC_COMM_WORLD,1_pPETSCINT,[field-1],PETSC_COPY_VALUES,pbcComps(numBC),err_PETSc) CHKERRQ(err_PETSc) @@ -320,15 +320,15 @@ end subroutine FEM_mechanical_init !> @brief solution for the FEM load step !-------------------------------------------------------------------------------------------------- type(tSolutionState) function FEM_mechanical_solution( & - incInfoIn,timeinc,timeinc_old,fieldBC) + incInfoIn,Delta_t,Delta_t_prev,mechBC) !-------------------------------------------------------------------------------------------------- ! input data for solution real(pREAL), intent(in) :: & - timeinc, & !< increment in time for current solution - timeinc_old !< increment in time of last increment - type(tFieldBC), intent(in) :: & - fieldBC + Delta_t, & !< increment in time for current solution + Delta_t_prev !< increment in time of last increment + type(tMechBC), intent(in) :: & + mechBC character(len=*), intent(in) :: & incInfoIn @@ -339,8 +339,8 @@ type(tSolutionState) function FEM_mechanical_solution( & FEM_mechanical_solution%converged = .false. !-------------------------------------------------------------------------------------------------- ! set module wide availabe data - params%timeinc = timeinc - params%fieldBC = fieldBC + params%Delta_t = Delta_t + params%mechBC = mechBC call SNESSolve(mechanical_snes,PETSC_NULL_VEC,solution,err_PETSc) ! solve mechanical_snes based on solution guess (result in solution) CHKERRQ(err_PETSc) @@ -407,13 +407,13 @@ subroutine FEM_mechanical_formResidual(dm_local,xx_local,f_local,dummy,err_PETSc call VecWAXPY(x_local,1.0_pREAL,xx_local,solution_local,err_PETSc) CHKERRQ(err_PETSc) do field = 1_pPETSCINT, dimPlex; do face = 1_pPETSCINT, mesh_Nboundaries - if (params%fieldBC%componentBC(field)%Mask(face)) then + if (params%mechBC%componentBC(field)%Mask(face)) then call DMGetStratumSize(dm_local,'Face Sets',mesh_boundaries(face),bcSize,err_PETSc) if (bcSize > 0) then call DMGetStratumIS(dm_local,'Face Sets',mesh_boundaries(face),bcPoints,err_PETSc) CHKERRQ(err_PETSc) call utilities_projectBCValues(x_local,section,0_pPETSCINT,field-1,bcPoints, & - 0.0_pREAL,params%fieldBC%componentBC(field)%Value(face),params%timeinc) + 0.0_pREAL,params%mechBC%componentBC(field)%Value(face),params%Delta_t) call ISDestroy(bcPoints,err_PETSc) CHKERRQ(err_PETSc) end if @@ -459,7 +459,7 @@ subroutine FEM_mechanical_formResidual(dm_local,xx_local,f_local,dummy,err_PETSc !-------------------------------------------------------------------------------------------------- ! evaluate constitutive response - call utilities_constitutiveResponse(params%timeinc,P_av,ForwardData) + call utilities_constitutiveResponse(params%Delta_t,P_av,ForwardData) call MPI_Allreduce(MPI_IN_PLACE,terminallyIll,1_MPI_INTEGER_KIND,MPI_LOGICAL,MPI_LOR,MPI_COMM_WORLD,err_MPI) if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error' ForwardData = .false. @@ -557,13 +557,13 @@ subroutine FEM_mechanical_formJacobian(dm_local,xx_local,Jac_pre,Jac,dummy,err_P call VecWAXPY(x_local,1.0_pREAL,xx_local,solution_local,err_PETSc) CHKERRQ(err_PETSc) do field = 1, dimPlex; do face = 1, mesh_Nboundaries - if (params%fieldBC%componentBC(field)%Mask(face)) then + if (params%mechBC%componentBC(field)%Mask(face)) then call DMGetStratumSize(dm_local,'Face Sets',mesh_boundaries(face),bcSize,err_PETSc) if (bcSize > 0) then call DMGetStratumIS(dm_local,'Face Sets',mesh_boundaries(face),bcPoints,err_PETSc) CHKERRQ(err_PETSc) call utilities_projectBCValues(x_local,section,0_pPETSCINT,field-1,bcPoints, & - 0.0_pREAL,params%fieldBC%componentBC(field)%Value(face),params%timeinc) + 0.0_pREAL,params%mechBC%componentBC(field)%Value(face),params%Delta_t) call ISDestroy(bcPoints,err_PETSc) CHKERRQ(err_PETSc) end if @@ -665,13 +665,13 @@ end subroutine FEM_mechanical_formJacobian !-------------------------------------------------------------------------------------------------- !> @brief forwarding routine !-------------------------------------------------------------------------------------------------- -subroutine FEM_mechanical_forward(guess,timeinc,timeinc_old,fieldBC) +subroutine FEM_mechanical_forward(guess,Delta_t,Delta_t_prev,mechBC) - type(tFieldBC), intent(in) :: & - fieldBC + type(tMechBC), intent(in) :: & + mechBC real(pREAL), intent(in) :: & - timeinc_old, & - timeinc + Delta_t_prev, & + Delta_t logical, intent(in) :: & guess @@ -702,13 +702,13 @@ subroutine FEM_mechanical_forward(guess,timeinc,timeinc_old,fieldBC) call VecAXPY(solution_local,1.0_pREAL,x_local,err_PETSc) CHKERRQ(err_PETSc) do field = 1, dimPlex; do face = 1, mesh_Nboundaries - if (fieldBC%componentBC(field)%Mask(face)) then + if (mechBC%componentBC(field)%Mask(face)) then call DMGetStratumSize(dm_local,'Face Sets',mesh_boundaries(face),bcSize,err_PETSc) if (bcSize > 0) then call DMGetStratumIS(dm_local,'Face Sets',mesh_boundaries(face),bcPoints,err_PETSc) CHKERRQ(err_PETSc) call utilities_projectBCValues(solution_local,section,0_pPETSCINT,field-1,bcPoints, & - 0.0_pREAL,fieldBC%componentBC(field)%Value(face),timeinc_old) + 0.0_pREAL,mechBC%componentBC(field)%Value(face),Delta_t_prev) call ISDestroy(bcPoints,err_PETSc) CHKERRQ(err_PETSc) end if @@ -721,12 +721,12 @@ subroutine FEM_mechanical_forward(guess,timeinc,timeinc_old,fieldBC) ! update rate and forward last inc call VecCopy(solution,solution_rate,err_PETSc) CHKERRQ(err_PETSc) - call VecScale(solution_rate,timeinc_old**(-1),err_PETSc) + call VecScale(solution_rate,Delta_t_prev**(-1),err_PETSc) CHKERRQ(err_PETSc) end if call VecCopy(solution_rate,solution,err_PETSc) CHKERRQ(err_PETSc) - call VecScale(solution,timeinc,err_PETSc) + call VecScale(solution,Delta_t,err_PETSc) CHKERRQ(err_PETSc) end subroutine FEM_mechanical_forward