223031012
/
DAMASK_EICMD
mirror of https://github.com/achalhp/DAMASK_EICMD.git
1
0
Fork 0
Code Issues Projects Releases Wiki Activity

* checked and corrected parallelization of code. compiles now successfully, but simply aborts computation with first parallel directive without any comment :-((( 

* also put a call to constitutive_microstructure at the start of each crystallite_integration subroutine like it was before. need that for nonlocal model in case of crystallite cutback
This commit is contained in:
Christoph Kords 2010-11-03 17:22:48 +00:00
parent 0dd99cb965
commit 405d5529e7
5 changed files with 940 additions and 919 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

8
code/CPFEM.f90
View File

@ -317,9 +317,8 @@ subroutine CPFEM_general(mode, ffn, ffn1, Temperature, dt, element, IP, cauchySt
! CPFEM_odd_jacobian
cp_en = mesh_FEasCP('elem',element)
selectiveDebugger = (cp_en == debug_e .and. IP == debug_i)
if (selectiveDebugger) then
if (selectiveDebugger .and. cp_en == debug_e .and. IP == debug_i) then
!$OMP CRITICAL (write2out)
write(6,*)
write(6,'(a)') '#######################################################'
@ -350,7 +349,7 @@ subroutine CPFEM_general(mode, ffn, ffn1, Temperature, dt, element, IP, cauchySt
j = 1:mesh_maxNips, &
k = 1:mesh_NcpElems ) &
constitutive_state0(i,j,k)%p = constitutive_state(i,j,k)%p ! microstructure of crystallites
if (selectiveDebugger) then
if (selectiveDebugger .and. cp_en == debug_e .and. IP == debug_i) then
!$OMP CRITICAL (write2out)
write(6,'(a,x,i8,x,i2,/,4(3(e20.8,x),/))') '<< cpfem >> AGED state of grain 1, element ip',&
cp_en,IP, constitutive_state(1,IP,cp_en)%p
@ -538,8 +537,7 @@ subroutine CPFEM_general(mode, ffn, ffn1, Temperature, dt, element, IP, cauchySt
pstress(:,:) = materialpoint_P(:,:,IP,cp_en)
dPdF(:,:,:,:) = materialpoint_dPdF(:,:,:,:,IP,cp_en)
selectiveDebugger = (cp_en == debug_e .and. IP == debug_i)
if (selectiveDebugger .and. mode < 6) then
if (selectiveDebugger .and. cp_en == debug_e .and. IP == debug_i .and. mode < 6) then
!$OMP CRITICAL (write2out)
write(6,'(a,x,i2,x,a,x,i4,/,6(f10.3,x)/)') 'stress/MPa at ip', IP, 'el', cp_en, cauchyStress/1e6
write(6,'(a,x,i2,x,a,x,i4,/,6(6(f10.3,x)/))') 'jacobian/GPa at ip', IP, 'el', cp_en, jacobian/1e9

1
code/constitutive.f90
View File

@ -142,7 +142,6 @@ subroutine constitutive_init()
myNgrains = homogenization_Ngrains(mesh_element(3,e))
do i = 1,FE_Nips(mesh_element(2,e)) ! loop over IPs
do g = 1,myNgrains ! loop over grains
selectiveDebugger = (e == debug_e .and. i == debug_i .and. g == debug_g)
myInstance = phase_constitutionInstance(material_phase(g,i,e))
select case(phase_constitution(material_phase(g,i,e)))

80
code/constitutive_nonlocal.f90
View File

@ -762,8 +762,7 @@ use math, only: math_Plain3333to99, &
math_det3x3, &
pi
use debug, only: debugger, &
verboseDebugger, &
selectiveDebugger
verboseDebugger
use mesh, only: mesh_NcpElems, &
mesh_maxNips, &
mesh_maxNipNeighbors, &
@ -964,25 +963,6 @@ do s = 1,ns
Tdislocation_v = Tdislocation_v + math_Mandel33to6( math_mul33x33(transpose(lattice2slip), math_mul33x33(sigma, lattice2slip) ) )
! if (selectiveDebugger .and. s==1) then
! write(6,*)
! write(6,'(a20,i1,x,i2,x,i5)') '::: microstructure ',g,ip,el
! write(6,*)
! write(6,'(a,/,3(3(f10.3,x)/))') 'position difference lattice / mu:', &
! transpose(neighboring_position((/1,3,5/),:)-neighboring_position((/2,4,6/),:)) * 1e6
! write(6,'(a,/,3(3(f10.3,x)/))') 'position difference slip system/ mu:', &
! math_mul33x33(lattice2slip,transpose(neighboring_position((/1,3,5/),:)-neighboring_position((/2,4,6/),:))) * 1e6
! write(6,*)
! write(6,'(a,/,2(3(e10.3,x)/))') 'excess dislo difference:', rhoExcessDifference
! write(6,*)
! write(6,'(a,/,2(3(e10.3,x)/))') 'disloGradients:', disloGradients
! write(6,*)
! write(6,'(a,/,3(21x,3(f10.4,x)/))') 'sigma / MPa:', transpose(sigma) * 1e-6
! write(6,'(a,/,3(21x,3(f10.4,x)/))') '2ndPK / MPa:', &
! transpose( math_mul33x33(transpose(lattice2slip), math_mul33x33(sigma, lattice2slip) ) ) * 1e-6
! write(6,*)
! endif
enddo
!**********************************************************************
@ -1009,8 +989,10 @@ use prec, only: pReal, &
use math, only: math_mul6x6, &
math_Mandel6to33
use debug, only: debugger, &
selectiveDebugger, &
verboseDebugger
verboseDebugger, &
debug_g, &
debug_i, &
debug_e
use mesh, only: mesh_NcpElems, &
mesh_maxNips
use material, only: homogenization_maxNgrains, &
@ -1080,7 +1062,7 @@ if ( Temperature > 0.0_pReal ) then
enddo
endif
!if (verboseDebugger .and. selectiveDebugger) then
!if (verboseDebugger .and. s) then
! !$OMP CRITICAL (write2out)
! write(6,*) '::: kinetics',g,ip,el
! write(6,*)
@ -1089,7 +1071,7 @@ endif
! write(6,'(a,/,12(f12.5,x),/)') 'tau / MPa', tau/1e6_pReal
! write(6,'(a,/,12(e12.5,x),/)') 'rhoForest / 1/m**2', rhoForest
! write(6,'(a,/,4(12(f12.5,x),/))') 'v / 1e-3m/s', constitutive_nonlocal_v(:,:,g,ip,el)*1e3
! !$OMPEND CRITICAL (write2out)
! !$OMP END CRITICAL (write2out)
!endif
endsubroutine
@ -1108,8 +1090,10 @@ use math, only: math_Plain3333to99, &
math_mul6x6, &
math_Mandel6to33
use debug, only: debugger, &
selectiveDebugger, &
verboseDebugger
verboseDebugger, &
debug_g, &
debug_i, &
debug_e
use mesh, only: mesh_NcpElems, &
mesh_maxNips
use material, only: homogenization_maxNgrains, &
@ -1201,7 +1185,7 @@ enddo
dLp_dTstar99 = math_Plain3333to99(dLp_dTstar3333)
!if (verboseDebugger .and. selectiveDebugger) then
!if (verboseDebugger .and. (debug_g==g .and. debug_i==i .and. debug_e==e)) then
! !$OMP CRITICAL (write2out)
! write(6,*) '::: LpandItsTangent',g,ip,el
! write(6,*)
@ -1210,7 +1194,7 @@ dLp_dTstar99 = math_Plain3333to99(dLp_dTstar3333)
! write(6,'(a,/,12(f12.5,x),/)') 'gdot total / 1e-3',gdotTotal*1e3_pReal
! write(6,'(a,/,3(3(f12.7,x)/))') 'Lp',Lp
! ! call flush(6)
! !$OMPEND CRITICAL (write2out)
! !$OMP END CRITICAL (write2out)
!endif
endsubroutine
@ -1228,7 +1212,9 @@ use prec, only: pReal, &
p_vec
use IO, only: IO_error
use debug, only: debugger, &
selectiveDebugger, &
debug_g, &
debug_i, &
debug_e, &
verboseDebugger
use math, only: math_norm3, &
math_mul6x6, &
@ -1358,11 +1344,11 @@ real(pReal) area, & ! area of
correction
logical, dimension(3) :: periodicSurfaceFlux ! flag indicating periodic fluxes at surfaces when surface normal points mainly in x, y and z direction respectively (in reference configuration)
if (verboseDebugger .and. selectiveDebugger) then
if (verboseDebugger .and. (debug_g==g .and. debug_i==ip .and. debug_e==el)) then
!$OMP CRITICAL (write2out)
write(6,*) '::: constitutive_nonlocal_dotState at ',g,ip,el
write(6,*)
!$OMPEND CRITICAL (write2out)
!$OMP END CRITICAL (write2out)
endif
select case(mesh_element(2,el))
@ -1416,12 +1402,12 @@ forall (s = 1:ns, t = 1:4, rhoSgl(s,t+4) * constitutive_nonlocal_v(s,t,g,ip,el)
gdot(s,t) = gdot(s,t) + abs(rhoSgl(s,t+4)) * constitutive_nonlocal_burgersPerSlipSystem(s,myInstance) &
* constitutive_nonlocal_v(s,t,g,ip,el)
if (verboseDebugger .and. selectiveDebugger) then
if (verboseDebugger .and. (debug_g==g .and. debug_i==ip .and. debug_e==el)) then
!$OMP CRITICAL (write2out)
write(6,'(a,/,10(12(e12.5,x),/))') 'rho / 1/m^2', rhoSgl, rhoDip
write(6,'(a,/,4(12(e12.5,x),/))') 'v / m/s', constitutive_nonlocal_v(:,:,g,ip,el)
write(6,'(a,/,4(12(e12.5,x),/))') 'gdot / 1/s',gdot
!$OMPEND CRITICAL (write2out)
!$OMP END CRITICAL (write2out)
endif
@ -1498,7 +1484,7 @@ detFe = math_det3x3(Fe(:,:,g,ip,el))
fluxdensity = rhoSgl(:,1:4) * constitutive_nonlocal_v(:,:,g,ip,el)
!if (selectiveDebugger) write(6,*) '--> dislocation flux <---'
!if ((debug_g==g .and. debug_i==ip .and. debug_e==el)) write(6,*) '--> dislocation flux <---'
do n = 1,FE_NipNeighbors(mesh_element(2,el)) ! loop through my neighbors
neighboring_el = mesh_ipNeighborhood(1,n,ip,el)
@ -1521,11 +1507,11 @@ do n = 1,FE_NipNeighbors(mesh_element(2,el))
surfaceNormal = surfaceNormal / math_norm3(surfaceNormal) ! normalize the surface normal to unit length
neighboring_rhoDotFlux = 0.0_pReal
! if (selectiveDebugger) write(6,'(a,x,i2)') 'neighbor',n
! if ((debug_g==g .and. debug_i==ip .and. debug_e==el)) write(6,'(a,x,i2)') 'neighbor',n
do s = 1,ns
! if (selectiveDebugger) write(6,'(a,x,i2)') ' system',s
! if ((debug_g==g .and. debug_i==ip .and. debug_e==el)) write(6,'(a,x,i2)') ' system',s
do t = 1,4
! if (selectiveDebugger) write(6,'(a,x,i2)') ' type',t
! if ((debug_g==g .and. debug_i==ip .and. debug_e==el)) write(6,'(a,x,i2)') ' type',t
c = (t + 1) / 2
topp = t + mod(t,2) - mod(t+1,2)
@ -1537,7 +1523,7 @@ do n = 1,FE_NipNeighbors(mesh_element(2,el))
if ( (opposite_el > 0 .and. opposite_ip > 0) &
.or. .not. all(periodicSurfaceFlux(maxloc(abs(mesh_ipAreaNormal(:,opposite_n,ip,el))))) ) then
rhoDotFlux(s,t) = rhoDotFlux(s,t) - lineLength / mesh_ipVolume(ip,el) ! subtract dislocation flux from cuurent mobile type
! if (selectiveDebugger) write(6,'(a,x,e12.5)') ' outgoing flux:', lineLength / mesh_ipVolume(ip,el)
! if ((debug_g==g .and. debug_i==ip .and. debug_e==el)) write(6,'(a,x,e12.5)') ' outgoing flux:', lineLength / mesh_ipVolume(ip,el)
endif
rhoDotFlux(s,t+4) = rhoDotFlux(s,t+4) + lineLength / mesh_ipVolume(ip,el) &
* (1.0_pReal - sum(constitutive_nonlocal_compatibility(c,:,s,n,ip,el)**2.0_pReal)) &
@ -1552,7 +1538,7 @@ do n = 1,FE_NipNeighbors(mesh_element(2,el))
neighboring_rhoDotFlux(:,topp) = neighboring_rhoDotFlux(:,topp) & ! ....transferring to opposite signed dislocation type at neighbor
+ lineLength / mesh_ipVolume(neighboring_ip,neighboring_el) &
* constitutive_nonlocal_compatibility(c,:,s,n,ip,el) ** 2.0_pReal
! if (selectiveDebugger) write(6,'(a,x,e12.5)') ' entering flux at neighbor:', lineLength / mesh_ipVolume(ip,el) &
! if ((debug_g==g .and. debug_i==ip .and. debug_e==el)) write(6,'(a,x,e12.5)') ' entering flux at neighbor:', lineLength / mesh_ipVolume(ip,el) &
! * sum(constitutive_nonlocal_compatibility(c,:,s,n,ip,el) ** 2.0_pReal)
endif
@ -1567,7 +1553,7 @@ do n = 1,FE_NipNeighbors(mesh_element(2,el))
constitutive_nonlocal_rhoDotFlux(:,:,g,neighboring_ip,neighboring_el) + neighboring_rhoDotFlux
dotState(g,neighboring_ip,neighboring_el)%p(1:10*ns) = &
dotState(g,neighboring_ip,neighboring_el)%p(1:10*ns) + reshape(neighboring_rhoDotFlux,(/10*ns/))
!$OMPEND CRITICAL (fluxes)
!$OMP END CRITICAL (fluxes)
else
neighboring_rhoDotFlux = 0.0_pReal
endif
@ -1577,7 +1563,7 @@ enddo ! neighbor loop
if (any(abs(rhoDotFlux) > 0.0_pReal)) then
!$OMP CRITICAL (fluxes)
constitutive_nonlocal_rhoDotFlux(:,:,g,ip,el) = constitutive_nonlocal_rhoDotFlux(:,:,g,ip,el) + rhoDotFlux
!$OMPEND CRITICAL (fluxes)
!$OMP END CRITICAL (fluxes)
endif
@ -1679,7 +1665,7 @@ forall (t = 1:10) &
! + rhoDotDipole2SingleStressChange(:,t)
! + rhoDotSingle2DipoleStressChange(:,t)
if (verboseDebugger .and. selectiveDebugger) then
if (verboseDebugger .and. (debug_g==g .and. debug_i==ip .and. debug_e==el)) then
!$OMP CRITICAL (write2out)
write(6,'(a,/,8(12(e12.5,x),/))') 'dislocation remobilization', rhoDotRemobilization(:,1:8) * timestep
write(6,'(a,/,4(12(e12.5,x),/))') 'dislocation multiplication', rhoDotMultiplication(:,1:4) * timestep
@ -1693,12 +1679,12 @@ if (verboseDebugger .and. selectiveDebugger) then
write(6,'(a,/,10(12(f12.7,x),/))') 'relative density change', rhoDot(:,1:8) * timestep / (abs(rhoSgl)+1.0e-10), &
rhoDot(:,9:10) * timestep / (rhoDip+1.0e-10)
write(6,*)
!$OMPEND CRITICAL (write2out)
!$OMP END CRITICAL (write2out)
endif
!$OMP CRITICAL (copy2dotState)
dotState(g,ip,el)%p(1:10*ns) = dotState(g,ip,el)%p(1:10*ns) + reshape(rhoDot,(/10*ns/))
!$OMPEND CRITICAL (copy2dotState)
!$OMP END CRITICAL (copy2dotState)
endsubroutine
@ -1734,8 +1720,7 @@ use lattice, only: lattice_sn, &
lattice_st
use debug, only: debugger, &
debug_e, debug_i, debug_g, &
verboseDebugger, &
selectiveDebugger
verboseDebugger
implicit none
@ -1777,7 +1762,6 @@ logical, dimension(maxval(constitutive_nonlocal_totalNslip)) :: &
compatibilityMask
selectiveDebugger = (debug_i==i .and. debug_e==e)
myPhase = material_phase(1,i,e)
myInstance = phase_constitutionInstance(myPhase)
myStructure = constitutive_nonlocal_structure(myInstance)

1682
code/crystallite.f90
View File

File diff suppressed because it is too large Load Diff

48
code/homogenization.f90
View File

@ -261,7 +261,6 @@ subroutine materialpoint_stressAndItsTangent(&
crystallite_stressAndItsTangent, &
crystallite_orientations
use debug, only: debugger, &
selectiveDebugger, &
verboseDebugger, &
debug_e, &
debug_i, &
@ -290,7 +289,7 @@ subroutine materialpoint_stressAndItsTangent(&
endif
!$OMP PARALLEL DO
!$OMP PARALLEL DO PRIVATE(myNgrains)
do e = FEsolving_execElem(1),FEsolving_execElem(2) ! iterate over elements to be processed
myNgrains = homogenization_Ngrains(mesh_element(3,e))
do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) ! iterate over IPs of this element to be processed
@ -324,20 +323,18 @@ subroutine materialpoint_stressAndItsTangent(&
do while (.not. terminallyIll .and. &
any(materialpoint_subStep(:,FEsolving_execELem(1):FEsolving_execElem(2)) > subStepMinHomog)) ! cutback loop for material points
!$OMP PARALLEL DO
!$OMP PARALLEL DO PRIVATE(myNgrains)
do e = FEsolving_execElem(1),FEsolving_execElem(2) ! iterate over elements to be processed
myNgrains = homogenization_Ngrains(mesh_element(3,e))
do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) ! iterate over IPs of this element to be processed
selectiveDebugger = (e == debug_e .and. i == debug_i)
if ( materialpoint_converged(i,e) ) then
if (verboseDebugger .and. selectiveDebugger) then
if (verboseDebugger .and. (e == debug_e .and. i == debug_i)) then
!$OMP CRITICAL (write2out)
write(6,'(a,x,f10.8,x,a,x,f10.8,x,a,/)') '°°° winding forward from', &
materialpoint_subFrac(i,e), 'to current materialpoint_subFrac', &
materialpoint_subFrac(i,e)+materialpoint_subStep(i,e),'in materialpoint_stressAndItsTangent'
!$OMPEND CRITICAL (write2out)
!$OMP END CRITICAL (write2out)
endif
! calculate new subStep and new subFrac
@ -363,7 +360,7 @@ subroutine materialpoint_stressAndItsTangent(&
!$OMP CRITICAL (distributionHomog)
debug_MaterialpointLoopDistribution(min(nHomog+1,NiterationHomog)) = &
debug_MaterialpointLoopDistribution(min(nHomog+1,NiterationHomog)) + 1
!$OMPEND CRITICAL (distributionHomog)
!$OMP END CRITICAL (distributionHomog)
endif
! materialpoint didn't converge, so we need a cutback here
@ -371,16 +368,18 @@ subroutine materialpoint_stressAndItsTangent(&
if ( (myNgrains == 1_pInt .and. materialpoint_subStep(i,e) <= 1.0 ) .or. & ! single grain already tried internal subStepping in crystallite
subStepSizeHomog * materialpoint_subStep(i,e) <= subStepMinHomog ) then ! would require too small subStep
! cutback makes no sense and...
terminallyIll = .true. ! ...one kills all
!$OMP CRITICAL (setTerminallyIll)
terminallyIll = .true. ! ...one kills all
!$OMP END CRITICAL (setTerminallyIll)
else ! cutback makes sense
materialpoint_subStep(i,e) = subStepSizeHomog * materialpoint_subStep(i,e) ! crystallite had severe trouble, so do a significant cutback
! <<modified to add more flexibility in cutback>>
if (verboseDebugger .and. selectiveDebugger) then
if (verboseDebugger .and. (e == debug_e .and. i == debug_i)) then
!$OMP CRITICAL (write2out)
write(6,'(a,x,f10.8,/)') '°°° cutback step in materialpoint_stressAndItsTangent with new materialpoint_subStep:',&
materialpoint_subStep(i,e)
!$OMPEND CRITICAL (write2out)
!$OMP END CRITICAL (write2out)
endif
! restore...
@ -405,7 +404,7 @@ subroutine materialpoint_stressAndItsTangent(&
endif
enddo ! loop IPs
enddo ! loop elements
!$OMP END PARALLEL DO
!$OMP END PARALLEL DO
!* Checks for cutback/substepping loops: added <<<updated 31.07.2009>>>
! write (6,'(a,/,8(L,x))') 'MP exceeds substep min',materialpoint_subStep(:,FEsolving_execELem(1):FEsolving_execElem(2)) > subStepMinHomog
@ -434,7 +433,7 @@ subroutine materialpoint_stressAndItsTangent(&
! homogenization_state
! results in crystallite_partionedF
!$OMP PARALLEL DO
!$OMP PARALLEL DO PRIVATE(myNgrains)
do e = FEsolving_execElem(1),FEsolving_execElem(2) ! iterate over elements to be processed
myNgrains = homogenization_Ngrains(mesh_element(3,e))
do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) ! iterate over IPs of this element to be processed
@ -448,7 +447,7 @@ subroutine materialpoint_stressAndItsTangent(&
endif
enddo
enddo
!$OMP END PARALLEL DO
!$OMP END PARALLEL DO
! write(6,'(a,/,125(8(8(l,x),2x),/))') 'crystallite request with updated partitioning', crystallite_requested
@ -462,7 +461,7 @@ subroutine materialpoint_stressAndItsTangent(&
! --+>> state update <<+--
!$OMP PARALLEL DO
!$OMP PARALLEL DO
do e = FEsolving_execElem(1),FEsolving_execElem(2) ! iterate over elements to be processed
do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) ! iterate over IPs of this element to be processed
if ( materialpoint_requested(i,e) .and. &
@ -473,13 +472,16 @@ subroutine materialpoint_stressAndItsTangent(&
materialpoint_doneAndHappy(:,i,e) = homogenization_updateState(i,e)
endif
materialpoint_converged(i,e) = all(materialpoint_doneAndHappy(:,i,e)) ! converged if done and happy
if (materialpoint_converged(i,e)) & ! added <<<updated 31.07.2009>>>
debug_MaterialpointStateLoopdistribution(NiterationMPstate) = &
debug_MaterialpointStateLoopdistribution(NiterationMPstate) + 1
if (materialpoint_converged(i,e)) then
!$OMP CRITICAL (distributionMPState)
debug_MaterialpointStateLoopdistribution(NiterationMPstate) = &
debug_MaterialpointStateLoopdistribution(NiterationMPstate) + 1
!$OMP END CRITICAL (distributionMPState)
endif
endif
enddo
enddo
!$OMP END PARALLEL DO
!$OMP END PARALLEL DO
! write(6,'(a,/,125(8(l,x),/))') 'material point done', materialpoint_doneAndHappy(1,:,:)
! write(6,'(a,/,125(8(l,x),/))') 'material point converged', materialpoint_converged
@ -493,13 +495,13 @@ subroutine materialpoint_stressAndItsTangent(&
if (.not. terminallyIll ) then
call crystallite_orientations() ! calculate crystal orientations
!$OMP PARALLEL DO
!$OMP PARALLEL DO
do e = FEsolving_execElem(1),FEsolving_execElem(2) ! iterate over elements to be processed
do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) ! iterate over IPs of this element to be processed
call homogenization_averageStressAndItsTangent(i,e)
call homogenization_averageTemperature(i,e)
enddo; enddo
!$OMP END PARALLEL DO
!$OMP END PARALLEL DO
if (debugger) then
write (6,*)
@ -533,7 +535,7 @@ subroutine materialpoint_postResults(dt)
real(pReal), intent(in) :: dt
integer(pInt) g,i,e,c,d,myNgrains,myCrystallite
!$OMP PARALLEL DO
!$OMP PARALLEL DO PRIVATE(myNgrains,myCrystallite,c,d)
do e = FEsolving_execElem(1),FEsolving_execElem(2) ! iterate over elements to be processed
myNgrains = homogenization_Ngrains(mesh_element(3,e))
myCrystallite = microstructure_crystallite(mesh_element(4,e))
@ -553,7 +555,7 @@ subroutine materialpoint_postResults(dt)
enddo
enddo
enddo
!$OMP END PARALLEL DO
!$OMP END PARALLEL DO
endsubroutine