From 94843becdafeab9d743d5496556a5b749c15cc3c Mon Sep 17 00:00:00 2001 From: Abisheik Panneerselvam Date: Wed, 14 Jul 2021 19:33:54 +0200 Subject: [PATCH 01/63] writing ip displacements (algorithm contain bugs) --- PRIVATE | 2 +- src/mesh/mesh_mech_FEM.f90 | 43 ++++++++++++++++++++++++++++++++++---- 2 files changed, 40 insertions(+), 5 deletions(-) diff --git a/PRIVATE b/PRIVATE index 1ca9a0e9f..918eed1d9 160000 --- a/PRIVATE +++ b/PRIVATE @@ -1 +1 @@ -Subproject commit 1ca9a0e9f2b333d3244b1ab44480380b3b22bebe +Subproject commit 918eed1d9f67eed75b4a4c66945c8c3053fa10fb diff --git a/src/mesh/mesh_mech_FEM.f90 b/src/mesh/mesh_mech_FEM.f90 index 6fa2f668b..43b997b70 100644 --- a/src/mesh/mesh_mech_FEM.f90 +++ b/src/mesh/mesh_mech_FEM.f90 @@ -670,7 +670,7 @@ end subroutine FEM_mechanical_converged !-------------------------------------------------------------------------------------------------- -!> @brief Calculate current coordinates (FEM nodal coordinates only at the moment) +!> @brief Calculate current coordinates (both nodal and ip coordinates) !-------------------------------------------------------------------------------------------------- subroutine FEM_mechanical_updateCoords() @@ -678,28 +678,63 @@ subroutine FEM_mechanical_updateCoords() nodeCoords_linear !< nodal coordinates (dimPlex*Nnodes) real(pReal), pointer, dimension(:,:) :: & nodeCoords !< nodal coordinates (3,Nnodes) + real(pReal), pointer, dimension(:,:,:) :: & + ipCoords !< ip coordinates (3,nQuadrature,mesh_NcpElems) DM :: dm_local Vec :: x_local PetscErrorCode :: ierr - PetscInt :: dimPlex, pStart, pEnd, p, s, e + PetscInt :: dimPlex, pStart, pEnd, p, s, e,& + cellStart, cellEnd, c, qPt, comp, nc, nqpoints,qOffset,nOffset,n,Nnodes PetscSection :: section + PetscFE :: mechFE + PetscDS :: mechDS + PetscQuadrature :: mechQuad + PetscReal, dimension(:), pointer :: qPoints, qWeights, basisField, basisFieldDer + PetscScalar, dimension(:), pointer :: x_scal call SNESGetDM(mechanical_snes,dm_local,ierr); CHKERRQ(ierr) + call DMGetDS(dm_local,mechQuad,ierr); CHKERRQ(ierr) call DMGetLocalSection(dm_local,section,ierr); CHKERRQ(ierr) call DMGetLocalVector(dm_local,x_local,ierr); CHKERRQ(ierr) call DMGetDimension(dm_local,dimPlex,ierr); CHKERRQ(ierr) + + ! write cell vertex displacements call DMPlexGetDepthStratum(dm_local,0,pStart,pEnd,ierr); CHKERRQ(ierr) allocate(nodeCoords(3,pStart:pEnd-1),source=0.0_pReal) call VecGetArrayF90(x_local,nodeCoords_linear,ierr); CHKERRQ(ierr) - do p=pStart, pEnd-1 call DMPlexGetPointLocal(dm_local, p, s, e, ierr); CHKERRQ(ierr) nodeCoords(1:dimPlex,p)=nodeCoords_linear(s+1:e) end do - call discretization_setNodeCoords(nodeCoords) call VecRestoreArrayF90(x_local,nodeCoords_linear,ierr); CHKERRQ(ierr) + + ! write ip displacements + call DMPlexGetHeightStratum(dm_local,0,cellStart,cellEnd,ierr); CHKERRQ(ierr) + call PetscDSGetTabulation(mechQuad,0,basisField,basisFieldDer,ierr); CHKERRQ(ierr) + allocate(ipCoords(0:2,0:nQuadrature-1,0:mesh_NcpElems-1),source=0.0_pReal) + ipCoords=0.0_pReal + do c=cellStart,cellEnd-1 + qOffset=0 + call DMPlexVecGetClosure(dm_local,section,x_local,c,x_scal,ierr); CHKERRQ(ierr) !< get nodal coordinates of each element + Nnodes=size(x_scal)/dimplex !< no. of nodes for each element + do qPt=0,nQuadrature-1 + qOffset=qPt*(cellDof*dimPlex) + do comp=0,dimPlex-1 !< loop over components + nOffset=0 + do n=0,Nnodes-1 + ipCoords(comp,qPt,c)=ipCoords(comp,qPt,c)+sum(basisField(qOffset*dimplex+1:qOffset*dimplex+2)*& + x_scal(nOffset+1:nOffset+2)) + qOffset=qOffset+dimplex + nOffset=nOffset+dimPlex + enddo + qOffset=qPt*(cellDof*dimPlex)+1 + enddo + enddo + call DMPlexVecRestoreClosure(dm_local,section,x_local,c,x_scal,ierr); CHKERRQ(ierr) + end do + call discretization_setIPcoords(reshape(ipCoords,[3,mesh_NcpElems*nQuadrature])) call DMRestoreLocalVector(dm_local,x_local,ierr); CHKERRQ(ierr) end subroutine FEM_mechanical_updateCoords From 4d6a9a51ed8be27ef3b2d81411369a26f5dd5d89 Mon Sep 17 00:00:00 2001 From: Sharan Roongta Date: Thu, 15 Jul 2021 17:30:00 +0200 Subject: [PATCH 02/63] correct logic --- PRIVATE | 2 +- src/mesh/mesh_mech_FEM.f90 | 28 +++++++++++++++++----------- 2 files changed, 18 insertions(+), 12 deletions(-) diff --git a/PRIVATE b/PRIVATE index 918eed1d9..8fec909d1 160000 --- a/PRIVATE +++ b/PRIVATE @@ -1 +1 @@ -Subproject commit 918eed1d9f67eed75b4a4c66945c8c3053fa10fb +Subproject commit 8fec909d1931b092b223b0560dd30c3339c6e5a7 diff --git a/src/mesh/mesh_mech_FEM.f90 b/src/mesh/mesh_mech_FEM.f90 index 43b997b70..c63f92a2c 100644 --- a/src/mesh/mesh_mech_FEM.f90 +++ b/src/mesh/mesh_mech_FEM.f90 @@ -681,10 +681,17 @@ subroutine FEM_mechanical_updateCoords() real(pReal), pointer, dimension(:,:,:) :: & ipCoords !< ip coordinates (3,nQuadrature,mesh_NcpElems) + integer :: & + qPt, & + comp, & + nc, & + qOffset, & + nOffset + DM :: dm_local Vec :: x_local PetscErrorCode :: ierr - PetscInt :: dimPlex, pStart, pEnd, p, s, e,& + PetscInt :: dimPlex, pStart, pEnd, p, s, e, q, & cellStart, cellEnd, c, qPt, comp, nc, nqpoints,qOffset,nOffset,n,Nnodes PetscSection :: section PetscFE :: mechFE @@ -713,23 +720,22 @@ subroutine FEM_mechanical_updateCoords() ! write ip displacements call DMPlexGetHeightStratum(dm_local,0,cellStart,cellEnd,ierr); CHKERRQ(ierr) call PetscDSGetTabulation(mechQuad,0,basisField,basisFieldDer,ierr); CHKERRQ(ierr) - allocate(ipCoords(0:2,0:nQuadrature-1,0:mesh_NcpElems-1),source=0.0_pReal) - ipCoords=0.0_pReal + allocate(ipCoords(3,nQuadrature,mesh_NcpElems),source=0.0_pReal) do c=cellStart,cellEnd-1 qOffset=0 call DMPlexVecGetClosure(dm_local,section,x_local,c,x_scal,ierr); CHKERRQ(ierr) !< get nodal coordinates of each element - Nnodes=size(x_scal)/dimplex !< no. of nodes for each element do qPt=0,nQuadrature-1 - qOffset=qPt*(cellDof*dimPlex) + qOffset= qPt * (size(basisField)/nQuadrature) do comp=0,dimPlex-1 !< loop over components nOffset=0 - do n=0,Nnodes-1 - ipCoords(comp,qPt,c)=ipCoords(comp,qPt,c)+sum(basisField(qOffset*dimplex+1:qOffset*dimplex+2)*& - x_scal(nOffset+1:nOffset+2)) - qOffset=qOffset+dimplex - nOffset=nOffset+dimPlex + q = comp + do n=0,nBasis-1 + ipCoords(comp+1,qPt+1,c+1)=ipCoords(comp+1,qPt+1,c+1)+& + sum(basisField(qOffset+(q*dimPlex)+1:qOffset+(q*dimPlex)+dimPlex)*& + x_scal(nOffset+1:nOffset+dimPlex)) + q = q+dimPlex + nOffset = nOffset+dimPlex enddo - qOffset=qPt*(cellDof*dimPlex)+1 enddo enddo call DMPlexVecRestoreClosure(dm_local,section,x_local,c,x_scal,ierr); CHKERRQ(ierr) From f05ac877eb9eb557dddb8254049c5f1f5f644c03 Mon Sep 17 00:00:00 2001 From: Sharan Roongta Date: Thu, 15 Jul 2021 18:53:51 +0200 Subject: [PATCH 03/63] ip displacement added to reference file --- PRIVATE | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/PRIVATE b/PRIVATE index 8fec909d1..feed4bff5 160000 --- a/PRIVATE +++ b/PRIVATE @@ -1 +1 @@ -Subproject commit 8fec909d1931b092b223b0560dd30c3339c6e5a7 +Subproject commit feed4bff505c7d556393d35ac9938379eda6d2e4 From 97203ff551e468277f9980a1b7f0e4217edabc66 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sun, 18 Jul 2021 23:30:30 +0200 Subject: [PATCH 04/63] document and keep definitions together --- src/mesh/FEM_quadrature.f90 | 7 +++++++ src/prec.f90 | 5 ----- 2 files changed, 7 insertions(+), 5 deletions(-) diff --git a/src/mesh/FEM_quadrature.f90 b/src/mesh/FEM_quadrature.f90 index ed0bbef9a..3af535c2f 100644 --- a/src/mesh/FEM_quadrature.f90 +++ b/src/mesh/FEM_quadrature.f90 @@ -20,6 +20,10 @@ module FEM_quadrature -1.0_pReal, 1.0_pReal, -1.0_pReal, & -1.0_pReal, -1.0_pReal, 1.0_pReal], shape=[3,4]) + type :: group_float !< variable length datatype used for storage of state + real(pReal), dimension(:), pointer :: p + end type group_float + integer, dimension(2:3,maxOrder), public, protected :: & FEM_nQuadrature !< number of quadrature points for a given spatial dimension(2-3) and interpolation order(1-maxOrder) type(group_float), dimension(2:3,maxOrder), public, protected :: & @@ -39,6 +43,9 @@ subroutine FEM_quadrature_init print'(/,a)', ' <<<+- FEM_quadrature init -+>>>'; flush(6) + print*, 'L. Zhang et al., Journal of Computational Mathematics 27(1):89-96, 2009' + print*, 'https://www.jstor.org/stable/43693493' + !-------------------------------------------------------------------------------------------------- ! 2D linear FEM_nQuadrature(2,1) = 1 diff --git a/src/prec.f90 b/src/prec.f90 index 7613cfe46..733286a8a 100644 --- a/src/prec.f90 +++ b/src/prec.f90 @@ -27,11 +27,6 @@ module prec real(pReal), parameter :: tol_math_check = 1.0e-8_pReal !< tolerance for internal math self-checks (rotation) - - type :: group_float !< variable length datatype used for storage of state - real(pReal), dimension(:), pointer :: p - end type group_float - type :: tState integer :: & sizeState = 0, & !< size of state From 28179963a292877022aab8161569038e0e33313b Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Mon, 19 Jul 2021 13:04:16 +0200 Subject: [PATCH 05/63] avoid failing tests --- PRIVATE | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/PRIVATE b/PRIVATE index feed4bff5..806d65296 160000 --- a/PRIVATE +++ b/PRIVATE @@ -1 +1 @@ -Subproject commit feed4bff505c7d556393d35ac9938379eda6d2e4 +Subproject commit 806d6529674b21931bb42957492b3f063032febc From a654cd4fb1caa30bff7178986eeb43acf376acf4 Mon Sep 17 00:00:00 2001 From: Abisheik Panneerselvam Date: Mon, 19 Jul 2021 17:07:43 +0200 Subject: [PATCH 06/63] To avoid fortran double variable truncation for quadrature points --- PRIVATE | 2 +- src/mesh/FEM_quadrature.f90 | 25 ++++++++++++++----------- src/mesh/mesh_mech_FEM.f90 | 8 +++----- 3 files changed, 18 insertions(+), 17 deletions(-) diff --git a/PRIVATE b/PRIVATE index 806d65296..2bafc2dbc 160000 --- a/PRIVATE +++ b/PRIVATE @@ -1 +1 @@ -Subproject commit 806d6529674b21931bb42957492b3f063032febc +Subproject commit 2bafc2dbc9b3b6682ef7c77c193dcdfe55d6b9c2 diff --git a/src/mesh/FEM_quadrature.f90 b/src/mesh/FEM_quadrature.f90 index 3af535c2f..5aa7ff947 100644 --- a/src/mesh/FEM_quadrature.f90 +++ b/src/mesh/FEM_quadrature.f90 @@ -71,12 +71,12 @@ subroutine FEM_quadrature_init FEM_nQuadrature(2,3) = 6 allocate(FEM_quadrature_weights(2,3)%p(6)) - FEM_quadrature_weights(2,3)%p(1:3) = 0.22338158967801146570_pReal - FEM_quadrature_weights(2,3)%p(4:6) = 0.10995174365532186764_pReal + FEM_quadrature_weights(2,3)%p(1:3) = 0.2233815896780115_pReal + FEM_quadrature_weights(2,3)%p(4:6) = 0.1099517436553218_pReal allocate(FEM_quadrature_points (2,3)%p(12)) - FEM_quadrature_points (2,3)%p(1:6) = permutationStar21([0.44594849091596488632_pReal]) - FEM_quadrature_points (2,3)%p(7:12)= permutationStar21([0.091576213509770743460_pReal]) + FEM_quadrature_points (2,3)%p(1:6) = permutationStar21([0.4459484909159649_pReal]) + FEM_quadrature_points (2,3)%p(7:12)= permutationStar21([0.0915762135097707_pReal]) !-------------------------------------------------------------------------------------------------- ! 2D quartic @@ -97,6 +97,7 @@ subroutine FEM_quadrature_init FEM_nQuadrature(2,5) = 16 allocate(FEM_quadrature_weights(2,5)%p(16)) + FEM_quadrature_weights(2,5)%p(1 ) = 0.14431560767779_pReal FEM_quadrature_weights(2,5)%p(2:4) = 0.09509163426728_pReal FEM_quadrature_weights(2,5)%p(5:7) = 0.10321737053472_pReal @@ -104,6 +105,7 @@ subroutine FEM_quadrature_init FEM_quadrature_weights(2,5)%p(11:16)= 0.02723031417443_pReal allocate(FEM_quadrature_points (2,5)%p(32)) + FEM_quadrature_points (2,5)%p(1:2) = permutationStar3([0.33333333333333_pReal]) FEM_quadrature_points (2,5)%p(3:8) = permutationStar21([0.45929258829272_pReal]) FEM_quadrature_points (2,5)%p(9:14) = permutationStar21([0.17056930775176_pReal]) @@ -128,21 +130,22 @@ subroutine FEM_quadrature_init FEM_quadrature_weights(3,2)%p(1:4) = 0.25_pReal allocate(FEM_quadrature_points (3,2)%p(12)) - FEM_quadrature_points (3,2)%p(1:12)= permutationStar31([0.13819660112501051518_pReal]) + + FEM_quadrature_points (3,2)%p(1:12)= permutationStar31([0.1381966011250105_pReal]) !-------------------------------------------------------------------------------------------------- ! 3D cubic FEM_nQuadrature(3,3) = 14 allocate(FEM_quadrature_weights(3,3)%p(14)) - FEM_quadrature_weights(3,3)%p(5:8) = 0.11268792571801585080_pReal - FEM_quadrature_weights(3,3)%p(1:4) = 0.073493043116361949544_pReal - FEM_quadrature_weights(3,3)%p(9:14) = 0.042546020777081466438_pReal + FEM_quadrature_weights(3,3)%p(5:8) = 0.11268792571801574_pReal + FEM_quadrature_weights(3,3)%p(1:4) = 0.0734930431163620_pReal + FEM_quadrature_weights(3,3)%p(9:14) = 0.0425460207770815_pReal allocate(FEM_quadrature_points (3,3)%p(42)) - FEM_quadrature_points (3,3)%p(1:12) = permutationStar31([0.092735250310891226402_pReal]) - FEM_quadrature_points (3,3)%p(13:24)= permutationStar31([0.31088591926330060980_pReal]) - FEM_quadrature_points (3,3)%p(25:42)= permutationStar22([0.045503704125649649492_pReal]) + FEM_quadrature_points (3,3)%p(1:12) = permutationStar31([0.0927352503108912_pReal]) + FEM_quadrature_points (3,3)%p(13:24)= permutationStar31([0.3108859192633006_pReal]) + FEM_quadrature_points (3,3)%p(25:42)= permutationStar22([0.0455037041256497_pReal]) !-------------------------------------------------------------------------------------------------- ! 3D quartic diff --git a/src/mesh/mesh_mech_FEM.f90 b/src/mesh/mesh_mech_FEM.f90 index 1bbb3d178..f63ec347f 100644 --- a/src/mesh/mesh_mech_FEM.f90 +++ b/src/mesh/mesh_mech_FEM.f90 @@ -690,13 +690,11 @@ subroutine FEM_mechanical_updateCoords() DM :: dm_local Vec :: x_local PetscErrorCode :: ierr - PetscInt :: dimPlex, pStart, pEnd, p, s, e, q, & - cellStart, cellEnd, c, nqpoints,n + PetscInt :: pStart, pEnd, p, s, e, q, & + cellStart, cellEnd, c, n PetscSection :: section - PetscFE :: mechFE - PetscDS :: mechDS PetscQuadrature :: mechQuad - PetscReal, dimension(:), pointer :: qPoints, qWeights, basisField, basisFieldDer + PetscReal, dimension(:), pointer :: basisField, basisFieldDer PetscScalar, dimension(:), pointer :: x_scal call SNESGetDM(mechanical_snes,dm_local,ierr); CHKERRQ(ierr) From 86d918367fe0d4b0ea95b8b3a6b8d3487a4d2bf8 Mon Sep 17 00:00:00 2001 From: Abisheik Panneerselvam Date: Tue, 20 Jul 2021 15:34:13 +0200 Subject: [PATCH 07/63] Round quadrature points --- src/mesh/FEM_quadrature.f90 | 49 ++++++++++++++++++++++--------------- 1 file changed, 29 insertions(+), 20 deletions(-) diff --git a/src/mesh/FEM_quadrature.f90 b/src/mesh/FEM_quadrature.f90 index 5aa7ff947..32516f4e1 100644 --- a/src/mesh/FEM_quadrature.f90 +++ b/src/mesh/FEM_quadrature.f90 @@ -71,10 +71,12 @@ subroutine FEM_quadrature_init FEM_nQuadrature(2,3) = 6 allocate(FEM_quadrature_weights(2,3)%p(6)) + FEM_quadrature_weights(2,3)%p(1:3) = 0.2233815896780115_pReal - FEM_quadrature_weights(2,3)%p(4:6) = 0.1099517436553218_pReal + FEM_quadrature_weights(2,3)%p(4:6) = 0.1099517436553219_pReal allocate(FEM_quadrature_points (2,3)%p(12)) + FEM_quadrature_points (2,3)%p(1:6) = permutationStar21([0.4459484909159649_pReal]) FEM_quadrature_points (2,3)%p(7:12)= permutationStar21([0.0915762135097707_pReal]) @@ -83,34 +85,35 @@ subroutine FEM_quadrature_init FEM_nQuadrature(2,4) = 12 allocate(FEM_quadrature_weights(2,4)%p(12)) - FEM_quadrature_weights(2,4)%p(1:3) = 0.11678627572638_pReal - FEM_quadrature_weights(2,4)%p(4:6) = 0.05084490637021_pReal - FEM_quadrature_weights(2,4)%p(7:12) = 0.08285107561837_pReal + + FEM_quadrature_weights(2,4)%p(1:3) = 0.1167862757263800_pReal + FEM_quadrature_weights(2,4)%p(4:6) = 0.0508449063702100_pReal + FEM_quadrature_weights(2,4)%p(7:12) = 0.0828510756183700_pReal allocate(FEM_quadrature_points (2,4)%p(24)) - FEM_quadrature_points (2,4)%p(1:6) = permutationStar21([0.24928674517091_pReal]) - FEM_quadrature_points (2,4)%p(7:12) = permutationStar21([0.06308901449150_pReal]) - FEM_quadrature_points (2,4)%p(13:24)= permutationStar111([0.31035245103378_pReal, 0.63650249912140_pReal]) + + FEM_quadrature_points (2,4)%p(1:6) = permutationStar21([0.2492867451709100_pReal]) + FEM_quadrature_points (2,4)%p(7:12) = permutationStar21([0.0630890144915000_pReal]) + FEM_quadrature_points (2,4)%p(13:24)= permutationStar111([0.3103524510337800_pReal, 0.6365024991214000_pReal]) !-------------------------------------------------------------------------------------------------- ! 2D quintic FEM_nQuadrature(2,5) = 16 allocate(FEM_quadrature_weights(2,5)%p(16)) - - FEM_quadrature_weights(2,5)%p(1 ) = 0.14431560767779_pReal - FEM_quadrature_weights(2,5)%p(2:4) = 0.09509163426728_pReal - FEM_quadrature_weights(2,5)%p(5:7) = 0.10321737053472_pReal - FEM_quadrature_weights(2,5)%p(8:10) = 0.03245849762320_pReal - FEM_quadrature_weights(2,5)%p(11:16)= 0.02723031417443_pReal + FEM_quadrature_weights(2,5)%p(1 ) = 0.1443156076777900_pReal + FEM_quadrature_weights(2,5)%p(2:4) = 0.0950916342672800_pReal + FEM_quadrature_weights(2,5)%p(5:7) = 0.1032173705347200_pReal + FEM_quadrature_weights(2,5)%p(8:10) = 0.0324584976232000_pReal + FEM_quadrature_weights(2,5)%p(11:16)= 0.0272303141744300_pReal allocate(FEM_quadrature_points (2,5)%p(32)) - FEM_quadrature_points (2,5)%p(1:2) = permutationStar3([0.33333333333333_pReal]) - FEM_quadrature_points (2,5)%p(3:8) = permutationStar21([0.45929258829272_pReal]) - FEM_quadrature_points (2,5)%p(9:14) = permutationStar21([0.17056930775176_pReal]) - FEM_quadrature_points (2,5)%p(15:20)= permutationStar21([0.05054722831703_pReal]) - FEM_quadrature_points (2,5)%p(21:32)= permutationStar111([0.26311282963464_pReal, 0.72849239295540_pReal]) + FEM_quadrature_points (2,5)%p(1:2) = permutationStar3([0.3333333333333300_pReal]) + FEM_quadrature_points (2,5)%p(3:8) = permutationStar21([0.4592925882927200_pReal]) + FEM_quadrature_points (2,5)%p(9:14) = permutationStar21([0.1705693077517600_pReal]) + FEM_quadrature_points (2,5)%p(15:20)= permutationStar21([0.0505472283170300_pReal]) + FEM_quadrature_points (2,5)%p(21:32)= permutationStar111([0.2631128296346400_pReal, 0.7284923929554000_pReal]) !-------------------------------------------------------------------------------------------------- ! 3D linear @@ -138,11 +141,13 @@ subroutine FEM_quadrature_init FEM_nQuadrature(3,3) = 14 allocate(FEM_quadrature_weights(3,3)%p(14)) - FEM_quadrature_weights(3,3)%p(5:8) = 0.11268792571801574_pReal + + FEM_quadrature_weights(3,3)%p(5:8) = 0.1126879257180159_pReal FEM_quadrature_weights(3,3)%p(1:4) = 0.0734930431163620_pReal FEM_quadrature_weights(3,3)%p(9:14) = 0.0425460207770815_pReal allocate(FEM_quadrature_points (3,3)%p(42)) + FEM_quadrature_points (3,3)%p(1:12) = permutationStar31([0.0927352503108912_pReal]) FEM_quadrature_points (3,3)%p(13:24)= permutationStar31([0.3108859192633006_pReal]) FEM_quadrature_points (3,3)%p(25:42)= permutationStar22([0.0455037041256497_pReal]) @@ -152,6 +157,7 @@ subroutine FEM_quadrature_init FEM_nQuadrature(3,4) = 35 allocate(FEM_quadrature_weights(3,4)%p(35)) + FEM_quadrature_weights(3,4)%p(1:4) = 0.0021900463965388_pReal FEM_quadrature_weights(3,4)%p(5:16) = 0.0143395670177665_pReal FEM_quadrature_weights(3,4)%p(17:22) = 0.0250305395686746_pReal @@ -159,6 +165,7 @@ subroutine FEM_quadrature_init FEM_quadrature_weights(3,4)%p(35) = 0.0931745731195340_pReal allocate(FEM_quadrature_points (3,4)%p(105)) + FEM_quadrature_points (3,4)%p(1:12) = permutationStar31([0.0267367755543735_pReal]) FEM_quadrature_points (3,4)%p(13:48) = permutationStar211([0.0391022406356488_pReal, 0.7477598884818090_pReal]) FEM_quadrature_points (3,4)%p(49:66) = permutationStar22([0.4547545999844830_pReal]) @@ -170,6 +177,7 @@ subroutine FEM_quadrature_init FEM_nQuadrature(3,5) = 56 allocate(FEM_quadrature_weights(3,5)%p(56)) + FEM_quadrature_weights(3,5)%p(1:4) = 0.0010373112336140_pReal FEM_quadrature_weights(3,5)%p(5:16) = 0.0096016645399480_pReal FEM_quadrature_weights(3,5)%p(17:28) = 0.0164493976798232_pReal @@ -178,6 +186,7 @@ subroutine FEM_quadrature_init FEM_quadrature_weights(3,5)%p(53:56) = 0.0366291366405108_pReal allocate(FEM_quadrature_points (3,5)%p(168)) + FEM_quadrature_points (3,5)%p(1:12) = permutationStar31([0.0149520651530592_pReal]) FEM_quadrature_points (3,5)%p(13:48) = permutationStar211([0.0340960211962615_pReal, 0.1518319491659370_pReal]) FEM_quadrature_points (3,5)%p(49:84) = permutationStar211([0.0462051504150017_pReal, 0.3549340560639790_pReal]) @@ -372,4 +381,4 @@ pure function permutationStar1111(point) result(qPt) end function permutationStar1111 -end module FEM_quadrature +end module FEM_quadrature \ No newline at end of file From 9427885b43b8b6ac66385de5da6c1e09d24a8005 Mon Sep 17 00:00:00 2001 From: Sharan Roongta Date: Wed, 21 Jul 2021 13:29:59 +0200 Subject: [PATCH 08/63] updated private --- PRIVATE | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/PRIVATE b/PRIVATE index 2bafc2dbc..4d696af2b 160000 --- a/PRIVATE +++ b/PRIVATE @@ -1 +1 @@ -Subproject commit 2bafc2dbc9b3b6682ef7c77c193dcdfe55d6b9c2 +Subproject commit 4d696af2bf78bb6d707fd51e958843b42b1279ff From f25dee4c3a7f1ddc73a577072fe7f2a87b78a50f Mon Sep 17 00:00:00 2001 From: Abisheik Panneerselvam Date: Wed, 21 Jul 2021 14:34:18 +0200 Subject: [PATCH 09/63] quad points updated and rounded based on literature values --- PRIVATE | 2 +- src/mesh/FEM_quadrature.f90 | 48 +++++++++++++++++-------------------- 2 files changed, 23 insertions(+), 27 deletions(-) diff --git a/PRIVATE b/PRIVATE index 4d696af2b..567067daf 160000 --- a/PRIVATE +++ b/PRIVATE @@ -1 +1 @@ -Subproject commit 4d696af2bf78bb6d707fd51e958843b42b1279ff +Subproject commit 567067daf05866c8e022fc0af6f441840f143b81 diff --git a/src/mesh/FEM_quadrature.f90 b/src/mesh/FEM_quadrature.f90 index 32516f4e1..e310cae06 100644 --- a/src/mesh/FEM_quadrature.f90 +++ b/src/mesh/FEM_quadrature.f90 @@ -70,50 +70,46 @@ subroutine FEM_quadrature_init ! 2D cubic FEM_nQuadrature(2,3) = 6 - allocate(FEM_quadrature_weights(2,3)%p(6)) - + allocate(FEM_quadrature_weights(2,3)%p(6)) FEM_quadrature_weights(2,3)%p(1:3) = 0.2233815896780115_pReal FEM_quadrature_weights(2,3)%p(4:6) = 0.1099517436553219_pReal allocate(FEM_quadrature_points (2,3)%p(12)) - FEM_quadrature_points (2,3)%p(1:6) = permutationStar21([0.4459484909159649_pReal]) - FEM_quadrature_points (2,3)%p(7:12)= permutationStar21([0.0915762135097707_pReal]) + FEM_quadrature_points (2,3)%p(7:12)= permutationStar21([0.09157621350977074_pReal]) !-------------------------------------------------------------------------------------------------- ! 2D quartic FEM_nQuadrature(2,4) = 12 - allocate(FEM_quadrature_weights(2,4)%p(12)) - - FEM_quadrature_weights(2,4)%p(1:3) = 0.1167862757263800_pReal - FEM_quadrature_weights(2,4)%p(4:6) = 0.0508449063702100_pReal - FEM_quadrature_weights(2,4)%p(7:12) = 0.0828510756183700_pReal + allocate(FEM_quadrature_weights(2,4)%p(12)) + FEM_quadrature_weights(2,4)%p(1:3) = 0.1167862757263794_pReal + FEM_quadrature_weights(2,4)%p(4:6) = 0.0508449063702068_pReal + FEM_quadrature_weights(2,4)%p(7:12) = 0.08285107561837358_pReal - allocate(FEM_quadrature_points (2,4)%p(24)) - - FEM_quadrature_points (2,4)%p(1:6) = permutationStar21([0.2492867451709100_pReal]) - FEM_quadrature_points (2,4)%p(7:12) = permutationStar21([0.0630890144915000_pReal]) - FEM_quadrature_points (2,4)%p(13:24)= permutationStar111([0.3103524510337800_pReal, 0.6365024991214000_pReal]) + allocate(FEM_quadrature_points (2,4)%p(24)) + FEM_quadrature_points (2,4)%p(1:6) = permutationStar21([0.2492867451709104_pReal]) + FEM_quadrature_points (2,4)%p(7:12) = permutationStar21([0.06308901449150223_pReal]) + FEM_quadrature_points (2,4)%p(13:24)= permutationStar111([0.3103524510337844_pReal, 0.05314504984481695_pReal]) !-------------------------------------------------------------------------------------------------- ! 2D quintic FEM_nQuadrature(2,5) = 16 allocate(FEM_quadrature_weights(2,5)%p(16)) - FEM_quadrature_weights(2,5)%p(1 ) = 0.1443156076777900_pReal - FEM_quadrature_weights(2,5)%p(2:4) = 0.0950916342672800_pReal - FEM_quadrature_weights(2,5)%p(5:7) = 0.1032173705347200_pReal - FEM_quadrature_weights(2,5)%p(8:10) = 0.0324584976232000_pReal - FEM_quadrature_weights(2,5)%p(11:16)= 0.0272303141744300_pReal + FEM_quadrature_weights(2,5)%p(1 ) = 0.1443156076777871_pReal + FEM_quadrature_weights(2,5)%p(2:4) = 0.09509163426728462_pReal + FEM_quadrature_weights(2,5)%p(5:7) = 0.1032173705347183_pReal + FEM_quadrature_weights(2,5)%p(8:10) = 0.03245849762319808_pReal + FEM_quadrature_weights(2,5)%p(11:16)= 0.02723031417443499_pReal allocate(FEM_quadrature_points (2,5)%p(32)) - FEM_quadrature_points (2,5)%p(1:2) = permutationStar3([0.3333333333333300_pReal]) - FEM_quadrature_points (2,5)%p(3:8) = permutationStar21([0.4592925882927200_pReal]) - FEM_quadrature_points (2,5)%p(9:14) = permutationStar21([0.1705693077517600_pReal]) - FEM_quadrature_points (2,5)%p(15:20)= permutationStar21([0.0505472283170300_pReal]) - FEM_quadrature_points (2,5)%p(21:32)= permutationStar111([0.2631128296346400_pReal, 0.7284923929554000_pReal]) + FEM_quadrature_points (2,5)%p(1:2) = permutationStar3([0.3333333333333333_pReal]) + FEM_quadrature_points (2,5)%p(3:8) = permutationStar21([0.4592925882927231_pReal]) + FEM_quadrature_points (2,5)%p(9:14) = permutationStar21([0.1705693077517602_pReal]) + FEM_quadrature_points (2,5)%p(15:20)= permutationStar21([0.0505472283170310_pReal]) + FEM_quadrature_points (2,5)%p(21:32)= permutationStar111([0.2631128296346381_pReal, 0.008394777409957605_pReal]) !-------------------------------------------------------------------------------------------------- ! 3D linear @@ -148,9 +144,9 @@ subroutine FEM_quadrature_init allocate(FEM_quadrature_points (3,3)%p(42)) - FEM_quadrature_points (3,3)%p(1:12) = permutationStar31([0.0927352503108912_pReal]) + FEM_quadrature_points (3,3)%p(1:12) = permutationStar31([0.09273525031089123_pReal]) FEM_quadrature_points (3,3)%p(13:24)= permutationStar31([0.3108859192633006_pReal]) - FEM_quadrature_points (3,3)%p(25:42)= permutationStar22([0.0455037041256497_pReal]) + FEM_quadrature_points (3,3)%p(25:42)= permutationStar22([0.04550370412564965_pReal]) !-------------------------------------------------------------------------------------------------- ! 3D quartic From b1ba64e6afbc7ba9e2b013c2750ccb7ef07d3d87 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 21 Jul 2021 17:51:58 +0200 Subject: [PATCH 10/63] understanding floating point precision double can use 16 or 17 significant figures - if the leading number is high (above 5 or so), only 16 figures are significant. - a trailing 1 is never significant --- src/mesh/FEM_quadrature.f90 | 70 ++++++++++++++++--------------------- 1 file changed, 31 insertions(+), 39 deletions(-) diff --git a/src/mesh/FEM_quadrature.f90 b/src/mesh/FEM_quadrature.f90 index e310cae06..f99bc58b3 100644 --- a/src/mesh/FEM_quadrature.f90 +++ b/src/mesh/FEM_quadrature.f90 @@ -51,65 +51,64 @@ subroutine FEM_quadrature_init FEM_nQuadrature(2,1) = 1 allocate(FEM_quadrature_weights(2,1)%p(1)) - FEM_quadrature_weights(2,1)%p(1) = 1.0_pReal + FEM_quadrature_weights(2,1)%p(1) = 1._pReal allocate(FEM_quadrature_points (2,1)%p(2)) - FEM_quadrature_points (2,1)%p(1:2) = permutationStar3([1.0_pReal/3.0_pReal]) + FEM_quadrature_points (2,1)%p(1:2) = permutationStar3([1._pReal/3._pReal]) !-------------------------------------------------------------------------------------------------- ! 2D quadratic FEM_nQuadrature(2,2) = 3 allocate(FEM_quadrature_weights(2,2)%p(3)) - FEM_quadrature_weights(2,2)%p(1:3) = 1.0_pReal/3.0_pReal + FEM_quadrature_weights(2,2)%p(1:3) = 1._pReal/3._pReal allocate(FEM_quadrature_points (2,2)%p(6)) - FEM_quadrature_points (2,2)%p(1:6) = permutationStar21([1.0_pReal/6.0_pReal]) + FEM_quadrature_points (2,2)%p(1:6) = permutationStar21([1._pReal/6._pReal]) !-------------------------------------------------------------------------------------------------- ! 2D cubic FEM_nQuadrature(2,3) = 6 allocate(FEM_quadrature_weights(2,3)%p(6)) - FEM_quadrature_weights(2,3)%p(1:3) = 0.2233815896780115_pReal - FEM_quadrature_weights(2,3)%p(4:6) = 0.1099517436553219_pReal + FEM_quadrature_weights(2,3)%p(1:3) = 2.2338158967801147e-1_pReal + FEM_quadrature_weights(2,3)%p(4:6) = 1.0995174365532187e-1_pReal allocate(FEM_quadrature_points (2,3)%p(12)) - FEM_quadrature_points (2,3)%p(1:6) = permutationStar21([0.4459484909159649_pReal]) - FEM_quadrature_points (2,3)%p(7:12)= permutationStar21([0.09157621350977074_pReal]) + FEM_quadrature_points (2,3)%p(1:6) = permutationStar21([4.4594849091596489e-1_pReal]) + FEM_quadrature_points (2,3)%p(7:12)= permutationStar21([9.157621350977074e-2_pReal]) !-------------------------------------------------------------------------------------------------- ! 2D quartic FEM_nQuadrature(2,4) = 12 allocate(FEM_quadrature_weights(2,4)%p(12)) - FEM_quadrature_weights(2,4)%p(1:3) = 0.1167862757263794_pReal - FEM_quadrature_weights(2,4)%p(4:6) = 0.0508449063702068_pReal - FEM_quadrature_weights(2,4)%p(7:12) = 0.08285107561837358_pReal + FEM_quadrature_weights(2,4)%p(1:3) = 1.1678627572637937e-1_pReal + FEM_quadrature_weights(2,4)%p(4:6) = 5.0844906370206817e-2_pReal + FEM_quadrature_weights(2,4)%p(7:12) = 8.285107561837358e-2_pReal allocate(FEM_quadrature_points (2,4)%p(24)) - FEM_quadrature_points (2,4)%p(1:6) = permutationStar21([0.2492867451709104_pReal]) - FEM_quadrature_points (2,4)%p(7:12) = permutationStar21([0.06308901449150223_pReal]) - FEM_quadrature_points (2,4)%p(13:24)= permutationStar111([0.3103524510337844_pReal, 0.05314504984481695_pReal]) + FEM_quadrature_points (2,4)%p(1:6) = permutationStar21([2.4928674517091042e-1_pReal]) + FEM_quadrature_points (2,4)%p(7:12) = permutationStar21([6.308901449150223e-2_pReal]) + FEM_quadrature_points (2,4)%p(13:24)= permutationStar111([3.1035245103378440e-1_pReal, 5.3145049844816947e-2_pReal]) !-------------------------------------------------------------------------------------------------- ! 2D quintic FEM_nQuadrature(2,5) = 16 allocate(FEM_quadrature_weights(2,5)%p(16)) - FEM_quadrature_weights(2,5)%p(1 ) = 0.1443156076777871_pReal - FEM_quadrature_weights(2,5)%p(2:4) = 0.09509163426728462_pReal - FEM_quadrature_weights(2,5)%p(5:7) = 0.1032173705347183_pReal - FEM_quadrature_weights(2,5)%p(8:10) = 0.03245849762319808_pReal - FEM_quadrature_weights(2,5)%p(11:16)= 0.02723031417443499_pReal + FEM_quadrature_weights(2,5)%p(1 ) = 1.4431560767778717e-1_pReal + FEM_quadrature_weights(2,5)%p(2:4) = 9.509163426728463e-2_pReal + FEM_quadrature_weights(2,5)%p(5:7) = 1.0321737053471825e-1_pReal + FEM_quadrature_weights(2,5)%p(8:10) = 3.2458497623198080e-2_pReal + FEM_quadrature_weights(2,5)%p(11:16)= 2.7230314174434994e-2_pReal allocate(FEM_quadrature_points (2,5)%p(32)) - - FEM_quadrature_points (2,5)%p(1:2) = permutationStar3([0.3333333333333333_pReal]) - FEM_quadrature_points (2,5)%p(3:8) = permutationStar21([0.4592925882927231_pReal]) - FEM_quadrature_points (2,5)%p(9:14) = permutationStar21([0.1705693077517602_pReal]) - FEM_quadrature_points (2,5)%p(15:20)= permutationStar21([0.0505472283170310_pReal]) - FEM_quadrature_points (2,5)%p(21:32)= permutationStar111([0.2631128296346381_pReal, 0.008394777409957605_pReal]) + FEM_quadrature_points (2,5)%p(1:2) = permutationStar3([1._pReal/3._pReal]) + FEM_quadrature_points (2,5)%p(3:8) = permutationStar21([4.5929258829272316e-1_pReal]) + FEM_quadrature_points (2,5)%p(9:14) = permutationStar21([1.705693077517602e-1_pReal]) + FEM_quadrature_points (2,5)%p(15:20)= permutationStar21([5.0547228317030975e-2_pReal]) + FEM_quadrature_points (2,5)%p(21:32)= permutationStar111([2.631128296346381e-1_pReal, 8.3947774099576053e-2_pReal]) !-------------------------------------------------------------------------------------------------- ! 3D linear @@ -129,7 +128,6 @@ subroutine FEM_quadrature_init FEM_quadrature_weights(3,2)%p(1:4) = 0.25_pReal allocate(FEM_quadrature_points (3,2)%p(12)) - FEM_quadrature_points (3,2)%p(1:12)= permutationStar31([0.1381966011250105_pReal]) !-------------------------------------------------------------------------------------------------- @@ -137,23 +135,20 @@ subroutine FEM_quadrature_init FEM_nQuadrature(3,3) = 14 allocate(FEM_quadrature_weights(3,3)%p(14)) - - FEM_quadrature_weights(3,3)%p(5:8) = 0.1126879257180159_pReal - FEM_quadrature_weights(3,3)%p(1:4) = 0.0734930431163620_pReal - FEM_quadrature_weights(3,3)%p(9:14) = 0.0425460207770815_pReal + FEM_quadrature_weights(3,3)%p(5:8) = 1.1268792571801585e-1_pReal + FEM_quadrature_weights(3,3)%p(1:4) = 7.3493043116361949e-2_pReal + FEM_quadrature_weights(3,3)%p(9:14) = 4.2546020777081467e-2_pReal allocate(FEM_quadrature_points (3,3)%p(42)) - - FEM_quadrature_points (3,3)%p(1:12) = permutationStar31([0.09273525031089123_pReal]) - FEM_quadrature_points (3,3)%p(13:24)= permutationStar31([0.3108859192633006_pReal]) - FEM_quadrature_points (3,3)%p(25:42)= permutationStar22([0.04550370412564965_pReal]) + FEM_quadrature_points (3,3)%p(1:12) = permutationStar31([9.273525031089123e-2_pReal]) + FEM_quadrature_points (3,3)%p(13:24)= permutationStar31([3.108859192633006e-1_pReal]) + FEM_quadrature_points (3,3)%p(25:42)= permutationStar22([4.5503704125649649e-2_pReal]) !-------------------------------------------------------------------------------------------------- ! 3D quartic FEM_nQuadrature(3,4) = 35 allocate(FEM_quadrature_weights(3,4)%p(35)) - FEM_quadrature_weights(3,4)%p(1:4) = 0.0021900463965388_pReal FEM_quadrature_weights(3,4)%p(5:16) = 0.0143395670177665_pReal FEM_quadrature_weights(3,4)%p(17:22) = 0.0250305395686746_pReal @@ -161,7 +156,6 @@ subroutine FEM_quadrature_init FEM_quadrature_weights(3,4)%p(35) = 0.0931745731195340_pReal allocate(FEM_quadrature_points (3,4)%p(105)) - FEM_quadrature_points (3,4)%p(1:12) = permutationStar31([0.0267367755543735_pReal]) FEM_quadrature_points (3,4)%p(13:48) = permutationStar211([0.0391022406356488_pReal, 0.7477598884818090_pReal]) FEM_quadrature_points (3,4)%p(49:66) = permutationStar22([0.4547545999844830_pReal]) @@ -173,7 +167,6 @@ subroutine FEM_quadrature_init FEM_nQuadrature(3,5) = 56 allocate(FEM_quadrature_weights(3,5)%p(56)) - FEM_quadrature_weights(3,5)%p(1:4) = 0.0010373112336140_pReal FEM_quadrature_weights(3,5)%p(5:16) = 0.0096016645399480_pReal FEM_quadrature_weights(3,5)%p(17:28) = 0.0164493976798232_pReal @@ -182,7 +175,6 @@ subroutine FEM_quadrature_init FEM_quadrature_weights(3,5)%p(53:56) = 0.0366291366405108_pReal allocate(FEM_quadrature_points (3,5)%p(168)) - FEM_quadrature_points (3,5)%p(1:12) = permutationStar31([0.0149520651530592_pReal]) FEM_quadrature_points (3,5)%p(13:48) = permutationStar211([0.0340960211962615_pReal, 0.1518319491659370_pReal]) FEM_quadrature_points (3,5)%p(49:84) = permutationStar211([0.0462051504150017_pReal, 0.3549340560639790_pReal]) @@ -377,4 +369,4 @@ pure function permutationStar1111(point) result(qPt) end function permutationStar1111 -end module FEM_quadrature \ No newline at end of file +end module FEM_quadrature From 000007ba596390ac52427e812805567897865cb5 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 21 Jul 2021 18:06:37 +0200 Subject: [PATCH 11/63] better test automatically --- src/mesh/FEM_quadrature.f90 | 21 ++++++++++++++++++++- 1 file changed, 20 insertions(+), 1 deletion(-) diff --git a/src/mesh/FEM_quadrature.f90 b/src/mesh/FEM_quadrature.f90 index f99bc58b3..31d4acec0 100644 --- a/src/mesh/FEM_quadrature.f90 +++ b/src/mesh/FEM_quadrature.f90 @@ -39,7 +39,7 @@ contains !-------------------------------------------------------------------------------------------------- !> @brief initializes FEM interpolation data !-------------------------------------------------------------------------------------------------- -subroutine FEM_quadrature_init +subroutine FEM_quadrature_init() print'(/,a)', ' <<<+- FEM_quadrature init -+>>>'; flush(6) @@ -182,6 +182,8 @@ subroutine FEM_quadrature_init FEM_quadrature_points (3,5)%p(121:156)= permutationStar211([0.3523052600879940_pReal, 0.0992057202494530_pReal]) FEM_quadrature_points (3,5)%p(157:168)= permutationStar31([0.1344783347929940_pReal]) + call selfTest + end subroutine FEM_quadrature_init @@ -369,4 +371,21 @@ pure function permutationStar1111(point) result(qPt) end function permutationStar1111 + +!-------------------------------------------------------------------------------------------------- +!> @brief Check correctness of quadrature weights. +!-------------------------------------------------------------------------------------------------- +subroutine selfTest + + integer :: o, d + + + do d = lbound(FEM_quadrature_weights,1), ubound(FEM_quadrature_weights,1) + do o = lbound(FEM_quadrature_weights(d,:),1), ubound(FEM_quadrature_weights(d,:),1) + if (dNeq(sum(FEM_quadrature_weights(d,o)%p),1.0_pReal,5e-15_pReal)) error stop 'quadrature weights' + enddo + enddo + +end subroutine selfTest + end module FEM_quadrature From c56979e2ad7d1b08accc75081ee98fc3e4840668 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 21 Jul 2021 18:09:44 +0200 Subject: [PATCH 12/63] single source of truth --- src/mesh/FEM_quadrature.f90 | 20 ++++++++++---------- 1 file changed, 10 insertions(+), 10 deletions(-) diff --git a/src/mesh/FEM_quadrature.f90 b/src/mesh/FEM_quadrature.f90 index 31d4acec0..2498e3947 100644 --- a/src/mesh/FEM_quadrature.f90 +++ b/src/mesh/FEM_quadrature.f90 @@ -50,7 +50,7 @@ subroutine FEM_quadrature_init() ! 2D linear FEM_nQuadrature(2,1) = 1 - allocate(FEM_quadrature_weights(2,1)%p(1)) + allocate(FEM_quadrature_weights(2,1)%p(FEM_nQuadrature(2,1))) FEM_quadrature_weights(2,1)%p(1) = 1._pReal allocate(FEM_quadrature_points (2,1)%p(2)) @@ -60,7 +60,7 @@ subroutine FEM_quadrature_init() ! 2D quadratic FEM_nQuadrature(2,2) = 3 - allocate(FEM_quadrature_weights(2,2)%p(3)) + allocate(FEM_quadrature_weights(2,2)%p(FEM_nQuadrature(2,2))) FEM_quadrature_weights(2,2)%p(1:3) = 1._pReal/3._pReal allocate(FEM_quadrature_points (2,2)%p(6)) @@ -70,7 +70,7 @@ subroutine FEM_quadrature_init() ! 2D cubic FEM_nQuadrature(2,3) = 6 - allocate(FEM_quadrature_weights(2,3)%p(6)) + allocate(FEM_quadrature_weights(2,3)%p(FEM_nQuadrature(2,3))) FEM_quadrature_weights(2,3)%p(1:3) = 2.2338158967801147e-1_pReal FEM_quadrature_weights(2,3)%p(4:6) = 1.0995174365532187e-1_pReal @@ -82,7 +82,7 @@ subroutine FEM_quadrature_init() ! 2D quartic FEM_nQuadrature(2,4) = 12 - allocate(FEM_quadrature_weights(2,4)%p(12)) + allocate(FEM_quadrature_weights(2,4)%p(FEM_nQuadrature(2,4))) FEM_quadrature_weights(2,4)%p(1:3) = 1.1678627572637937e-1_pReal FEM_quadrature_weights(2,4)%p(4:6) = 5.0844906370206817e-2_pReal FEM_quadrature_weights(2,4)%p(7:12) = 8.285107561837358e-2_pReal @@ -96,7 +96,7 @@ subroutine FEM_quadrature_init() ! 2D quintic FEM_nQuadrature(2,5) = 16 - allocate(FEM_quadrature_weights(2,5)%p(16)) + allocate(FEM_quadrature_weights(2,5)%p(FEM_nQuadrature(2,5))) FEM_quadrature_weights(2,5)%p(1 ) = 1.4431560767778717e-1_pReal FEM_quadrature_weights(2,5)%p(2:4) = 9.509163426728463e-2_pReal FEM_quadrature_weights(2,5)%p(5:7) = 1.0321737053471825e-1_pReal @@ -114,7 +114,7 @@ subroutine FEM_quadrature_init() ! 3D linear FEM_nQuadrature(3,1) = 1 - allocate(FEM_quadrature_weights(3,1)%p(1)) + allocate(FEM_quadrature_weights(3,1)%p(FEM_nQuadrature(3,1))) FEM_quadrature_weights(3,1)%p(1) = 1.0_pReal allocate(FEM_quadrature_points (3,1)%p(3)) @@ -124,7 +124,7 @@ subroutine FEM_quadrature_init() ! 3D quadratic FEM_nQuadrature(3,2) = 4 - allocate(FEM_quadrature_weights(3,2)%p(4)) + allocate(FEM_quadrature_weights(3,2)%p(FEM_nQuadrature(3,2))) FEM_quadrature_weights(3,2)%p(1:4) = 0.25_pReal allocate(FEM_quadrature_points (3,2)%p(12)) @@ -134,7 +134,7 @@ subroutine FEM_quadrature_init() ! 3D cubic FEM_nQuadrature(3,3) = 14 - allocate(FEM_quadrature_weights(3,3)%p(14)) + allocate(FEM_quadrature_weights(3,3)%p(FEM_nQuadrature(3,3))) FEM_quadrature_weights(3,3)%p(5:8) = 1.1268792571801585e-1_pReal FEM_quadrature_weights(3,3)%p(1:4) = 7.3493043116361949e-2_pReal FEM_quadrature_weights(3,3)%p(9:14) = 4.2546020777081467e-2_pReal @@ -148,7 +148,7 @@ subroutine FEM_quadrature_init() ! 3D quartic FEM_nQuadrature(3,4) = 35 - allocate(FEM_quadrature_weights(3,4)%p(35)) + allocate(FEM_quadrature_weights(3,4)%p(FEM_nQuadrature(3,4))) FEM_quadrature_weights(3,4)%p(1:4) = 0.0021900463965388_pReal FEM_quadrature_weights(3,4)%p(5:16) = 0.0143395670177665_pReal FEM_quadrature_weights(3,4)%p(17:22) = 0.0250305395686746_pReal @@ -166,7 +166,7 @@ subroutine FEM_quadrature_init() ! 3D quintic FEM_nQuadrature(3,5) = 56 - allocate(FEM_quadrature_weights(3,5)%p(56)) + allocate(FEM_quadrature_weights(3,5)%p(FEM_nQuadrature(3,5))) FEM_quadrature_weights(3,5)%p(1:4) = 0.0010373112336140_pReal FEM_quadrature_weights(3,5)%p(5:16) = 0.0096016645399480_pReal FEM_quadrature_weights(3,5)%p(17:28) = 0.0164493976798232_pReal From 5a04a1d661a384c8469f9b22bc3f03e9b5a40f15 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 21 Jul 2021 18:59:34 +0200 Subject: [PATCH 13/63] let the computer do the counting --- src/mesh/FEM_quadrature.f90 | 104 +++++++++++++++++------------------- 1 file changed, 50 insertions(+), 54 deletions(-) diff --git a/src/mesh/FEM_quadrature.f90 b/src/mesh/FEM_quadrature.f90 index 2498e3947..e76e9849f 100644 --- a/src/mesh/FEM_quadrature.f90 +++ b/src/mesh/FEM_quadrature.f90 @@ -9,7 +9,7 @@ module FEM_quadrature private integer, parameter :: & - maxOrder = 5 !< current max interpolation set at cubic (intended to be arbitrary) + maxOrder = 5 !< maximum integration order real(pReal), dimension(2,3), parameter :: & triangle = reshape([-1.0_pReal, -1.0_pReal, & 1.0_pReal, -1.0_pReal, & @@ -21,7 +21,7 @@ module FEM_quadrature -1.0_pReal, -1.0_pReal, 1.0_pReal], shape=[3,4]) type :: group_float !< variable length datatype used for storage of state - real(pReal), dimension(:), pointer :: p + real(pReal), dimension(:), allocatable :: p end type group_float integer, dimension(2:3,maxOrder), public, protected :: & @@ -53,8 +53,7 @@ subroutine FEM_quadrature_init() allocate(FEM_quadrature_weights(2,1)%p(FEM_nQuadrature(2,1))) FEM_quadrature_weights(2,1)%p(1) = 1._pReal - allocate(FEM_quadrature_points (2,1)%p(2)) - FEM_quadrature_points (2,1)%p(1:2) = permutationStar3([1._pReal/3._pReal]) + FEM_quadrature_points (2,1)%p = permutationStar3([1._pReal/3._pReal]) !-------------------------------------------------------------------------------------------------- ! 2D quadratic @@ -63,8 +62,7 @@ subroutine FEM_quadrature_init() allocate(FEM_quadrature_weights(2,2)%p(FEM_nQuadrature(2,2))) FEM_quadrature_weights(2,2)%p(1:3) = 1._pReal/3._pReal - allocate(FEM_quadrature_points (2,2)%p(6)) - FEM_quadrature_points (2,2)%p(1:6) = permutationStar21([1._pReal/6._pReal]) + FEM_quadrature_points (2,2)%p = permutationStar21([1._pReal/6._pReal]) !-------------------------------------------------------------------------------------------------- ! 2D cubic @@ -74,9 +72,9 @@ subroutine FEM_quadrature_init() FEM_quadrature_weights(2,3)%p(1:3) = 2.2338158967801147e-1_pReal FEM_quadrature_weights(2,3)%p(4:6) = 1.0995174365532187e-1_pReal - allocate(FEM_quadrature_points (2,3)%p(12)) - FEM_quadrature_points (2,3)%p(1:6) = permutationStar21([4.4594849091596489e-1_pReal]) - FEM_quadrature_points (2,3)%p(7:12)= permutationStar21([9.157621350977074e-2_pReal]) + FEM_quadrature_points (2,3)%p = [ & + permutationStar21([4.4594849091596489e-1_pReal]), & + permutationStar21([9.157621350977074e-2_pReal]) ] !-------------------------------------------------------------------------------------------------- ! 2D quartic @@ -87,10 +85,10 @@ subroutine FEM_quadrature_init() FEM_quadrature_weights(2,4)%p(4:6) = 5.0844906370206817e-2_pReal FEM_quadrature_weights(2,4)%p(7:12) = 8.285107561837358e-2_pReal - allocate(FEM_quadrature_points (2,4)%p(24)) - FEM_quadrature_points (2,4)%p(1:6) = permutationStar21([2.4928674517091042e-1_pReal]) - FEM_quadrature_points (2,4)%p(7:12) = permutationStar21([6.308901449150223e-2_pReal]) - FEM_quadrature_points (2,4)%p(13:24)= permutationStar111([3.1035245103378440e-1_pReal, 5.3145049844816947e-2_pReal]) + FEM_quadrature_points (2,4)%p = [ & + permutationStar21([2.4928674517091042e-1_pReal]), & + permutationStar21([6.308901449150223e-2_pReal]), & + permutationStar111([3.1035245103378440e-1_pReal, 5.3145049844816947e-2_pReal]) ] !-------------------------------------------------------------------------------------------------- ! 2D quintic @@ -103,22 +101,21 @@ subroutine FEM_quadrature_init() FEM_quadrature_weights(2,5)%p(8:10) = 3.2458497623198080e-2_pReal FEM_quadrature_weights(2,5)%p(11:16)= 2.7230314174434994e-2_pReal - allocate(FEM_quadrature_points (2,5)%p(32)) - FEM_quadrature_points (2,5)%p(1:2) = permutationStar3([1._pReal/3._pReal]) - FEM_quadrature_points (2,5)%p(3:8) = permutationStar21([4.5929258829272316e-1_pReal]) - FEM_quadrature_points (2,5)%p(9:14) = permutationStar21([1.705693077517602e-1_pReal]) - FEM_quadrature_points (2,5)%p(15:20)= permutationStar21([5.0547228317030975e-2_pReal]) - FEM_quadrature_points (2,5)%p(21:32)= permutationStar111([2.631128296346381e-1_pReal, 8.3947774099576053e-2_pReal]) + FEM_quadrature_points (2,5)%p = [ & + permutationStar3([1._pReal/3._pReal]), & + permutationStar21([4.5929258829272316e-1_pReal]), & + permutationStar21([1.705693077517602e-1_pReal]), & + permutationStar21([5.0547228317030975e-2_pReal]), & + permutationStar111([2.631128296346381e-1_pReal, 8.3947774099576053e-2_pReal]) ] !-------------------------------------------------------------------------------------------------- ! 3D linear FEM_nQuadrature(3,1) = 1 allocate(FEM_quadrature_weights(3,1)%p(FEM_nQuadrature(3,1))) - FEM_quadrature_weights(3,1)%p(1) = 1.0_pReal + FEM_quadrature_weights(3,1)%p(1) = 1.0_pReal - allocate(FEM_quadrature_points (3,1)%p(3)) - FEM_quadrature_points (3,1)%p(1:3)= permutationStar4([0.25_pReal]) + FEM_quadrature_points (3,1)%p = permutationStar4([0.25_pReal]) !-------------------------------------------------------------------------------------------------- ! 3D quadratic @@ -127,60 +124,59 @@ subroutine FEM_quadrature_init() allocate(FEM_quadrature_weights(3,2)%p(FEM_nQuadrature(3,2))) FEM_quadrature_weights(3,2)%p(1:4) = 0.25_pReal - allocate(FEM_quadrature_points (3,2)%p(12)) - FEM_quadrature_points (3,2)%p(1:12)= permutationStar31([0.1381966011250105_pReal]) + FEM_quadrature_points (3,2)%p = permutationStar31([0.1381966011250105_pReal]) !-------------------------------------------------------------------------------------------------- ! 3D cubic FEM_nQuadrature(3,3) = 14 allocate(FEM_quadrature_weights(3,3)%p(FEM_nQuadrature(3,3))) - FEM_quadrature_weights(3,3)%p(5:8) = 1.1268792571801585e-1_pReal FEM_quadrature_weights(3,3)%p(1:4) = 7.3493043116361949e-2_pReal + FEM_quadrature_weights(3,3)%p(5:8) = 1.1268792571801585e-1_pReal FEM_quadrature_weights(3,3)%p(9:14) = 4.2546020777081467e-2_pReal - allocate(FEM_quadrature_points (3,3)%p(42)) - FEM_quadrature_points (3,3)%p(1:12) = permutationStar31([9.273525031089123e-2_pReal]) - FEM_quadrature_points (3,3)%p(13:24)= permutationStar31([3.108859192633006e-1_pReal]) - FEM_quadrature_points (3,3)%p(25:42)= permutationStar22([4.5503704125649649e-2_pReal]) + FEM_quadrature_points (3,3)%p = [ & + permutationStar31([9.273525031089123e-2_pReal]), & + permutationStar31([3.108859192633006e-1_pReal]), & + permutationStar22([4.5503704125649649e-2_pReal]) ] !-------------------------------------------------------------------------------------------------- ! 3D quartic FEM_nQuadrature(3,4) = 35 allocate(FEM_quadrature_weights(3,4)%p(FEM_nQuadrature(3,4))) - FEM_quadrature_weights(3,4)%p(1:4) = 0.0021900463965388_pReal - FEM_quadrature_weights(3,4)%p(5:16) = 0.0143395670177665_pReal - FEM_quadrature_weights(3,4)%p(17:22) = 0.0250305395686746_pReal - FEM_quadrature_weights(3,4)%p(23:34) = 0.0479839333057554_pReal - FEM_quadrature_weights(3,4)%p(35) = 0.0931745731195340_pReal + FEM_quadrature_weights(3,4)%p(1:4) = 0.0021900463965388_pReal + FEM_quadrature_weights(3,4)%p(5:16) = 0.0143395670177665_pReal + FEM_quadrature_weights(3,4)%p(17:22) = 0.0250305395686746_pReal + FEM_quadrature_weights(3,4)%p(23:34) = 0.0479839333057554_pReal + FEM_quadrature_weights(3,4)%p(35) = 0.0931745731195340_pReal - allocate(FEM_quadrature_points (3,4)%p(105)) - FEM_quadrature_points (3,4)%p(1:12) = permutationStar31([0.0267367755543735_pReal]) - FEM_quadrature_points (3,4)%p(13:48) = permutationStar211([0.0391022406356488_pReal, 0.7477598884818090_pReal]) - FEM_quadrature_points (3,4)%p(49:66) = permutationStar22([0.4547545999844830_pReal]) - FEM_quadrature_points (3,4)%p(67:102) = permutationStar211([0.2232010379623150_pReal, 0.0504792790607720_pReal]) - FEM_quadrature_points (3,4)%p(103:105)= permutationStar4([0.25_pReal]) + FEM_quadrature_points (3,4)%p = [ & + permutationStar31([0.0267367755543735_pReal]), & + permutationStar211([0.0391022406356488_pReal, 0.7477598884818090_pReal]), & + permutationStar22([0.4547545999844830_pReal]), & + permutationStar211([0.2232010379623150_pReal, 0.0504792790607720_pReal]), & + permutationStar4([0.25_pReal]) ] !-------------------------------------------------------------------------------------------------- ! 3D quintic FEM_nQuadrature(3,5) = 56 allocate(FEM_quadrature_weights(3,5)%p(FEM_nQuadrature(3,5))) - FEM_quadrature_weights(3,5)%p(1:4) = 0.0010373112336140_pReal - FEM_quadrature_weights(3,5)%p(5:16) = 0.0096016645399480_pReal - FEM_quadrature_weights(3,5)%p(17:28) = 0.0164493976798232_pReal - FEM_quadrature_weights(3,5)%p(29:40) = 0.0153747766513310_pReal - FEM_quadrature_weights(3,5)%p(41:52) = 0.0293520118375230_pReal - FEM_quadrature_weights(3,5)%p(53:56) = 0.0366291366405108_pReal + FEM_quadrature_weights(3,5)%p(1:4) = 0.0010373112336140_pReal + FEM_quadrature_weights(3,5)%p(5:16) = 0.0096016645399480_pReal + FEM_quadrature_weights(3,5)%p(17:28) = 0.0164493976798232_pReal + FEM_quadrature_weights(3,5)%p(29:40) = 0.0153747766513310_pReal + FEM_quadrature_weights(3,5)%p(41:52) = 0.0293520118375230_pReal + FEM_quadrature_weights(3,5)%p(53:56) = 0.0366291366405108_pReal - allocate(FEM_quadrature_points (3,5)%p(168)) - FEM_quadrature_points (3,5)%p(1:12) = permutationStar31([0.0149520651530592_pReal]) - FEM_quadrature_points (3,5)%p(13:48) = permutationStar211([0.0340960211962615_pReal, 0.1518319491659370_pReal]) - FEM_quadrature_points (3,5)%p(49:84) = permutationStar211([0.0462051504150017_pReal, 0.3549340560639790_pReal]) - FEM_quadrature_points (3,5)%p(85:120) = permutationStar211([0.2281904610687610_pReal, 0.0055147549744775_pReal]) - FEM_quadrature_points (3,5)%p(121:156)= permutationStar211([0.3523052600879940_pReal, 0.0992057202494530_pReal]) - FEM_quadrature_points (3,5)%p(157:168)= permutationStar31([0.1344783347929940_pReal]) + FEM_quadrature_points (3,5)%p = [ & + permutationStar31([0.0149520651530592_pReal]), & + permutationStar211([0.0340960211962615_pReal, 0.1518319491659370_pReal]), & + permutationStar211([0.0462051504150017_pReal, 0.3549340560639790_pReal]), & + permutationStar211([0.2281904610687610_pReal, 0.0055147549744775_pReal]), & + permutationStar211([0.3523052600879940_pReal, 0.0992057202494530_pReal]), & + permutationStar31([0.1344783347929940_pReal]) ] call selfTest @@ -218,7 +214,7 @@ pure function permutationStar21(point) result(qPt) temp(:,2) = [point(1), 1.0_pReal - 2.0_pReal*point(1), point(1)] temp(:,3) = [1.0_pReal - 2.0_pReal*point(1), point(1), point(1)] - qPt = reshape(matmul(triangle, temp),[6]) + qPt = reshape(matmul(triangle, temp),[6]) end function permutationStar21 From 23077fd58c13f8f584a8e8ff30b6b5fa212ab1ce Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 21 Jul 2021 19:42:53 +0200 Subject: [PATCH 14/63] testing coordinates I don't understand why this pattern exists (with w(2) = 3, w(3) = 2), but it exists --- src/mesh/FEM_quadrature.f90 | 15 ++++++++++++--- 1 file changed, 12 insertions(+), 3 deletions(-) diff --git a/src/mesh/FEM_quadrature.f90 b/src/mesh/FEM_quadrature.f90 index e76e9849f..53aadbe6f 100644 --- a/src/mesh/FEM_quadrature.f90 +++ b/src/mesh/FEM_quadrature.f90 @@ -373,12 +373,21 @@ end function permutationStar1111 !-------------------------------------------------------------------------------------------------- subroutine selfTest - integer :: o, d - + integer :: o, d, n + real(pReal), dimension(2:3), parameter :: w = [3.0_pReal,2.0_pReal] do d = lbound(FEM_quadrature_weights,1), ubound(FEM_quadrature_weights,1) do o = lbound(FEM_quadrature_weights(d,:),1), ubound(FEM_quadrature_weights(d,:),1) - if (dNeq(sum(FEM_quadrature_weights(d,o)%p),1.0_pReal,5e-15_pReal)) error stop 'quadrature weights' + if (dNeq(sum(FEM_quadrature_weights(d,o)%p),1.0_pReal,5e-15_pReal)) & + error stop 'quadrature weights' + enddo + enddo + + do d = lbound(FEM_quadrature_points,1), ubound(FEM_quadrature_points,1) + do o = lbound(FEM_quadrature_points(d,:),1), ubound(FEM_quadrature_points(d,:),1) + n = size(FEM_quadrature_points(d,o)%p,1)/d + if (any(dNeq(sum(reshape(FEM_quadrature_points(d,o)%p,[d,n]),2),-real(n,pReal)/w(d),1.e-14_pReal))) & + error stop 'quadrature points' enddo enddo From b817c620a371e7dcb32b9172b756a5bd59e83879 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 21 Jul 2021 20:35:42 +0200 Subject: [PATCH 15/63] let the computer count --- src/mesh/FEM_quadrature.f90 | 147 ++++++++++++++++++++---------------- 1 file changed, 82 insertions(+), 65 deletions(-) diff --git a/src/mesh/FEM_quadrature.f90 b/src/mesh/FEM_quadrature.f90 index 53aadbe6f..278bfb5db 100644 --- a/src/mesh/FEM_quadrature.f90 +++ b/src/mesh/FEM_quadrature.f90 @@ -51,7 +51,7 @@ subroutine FEM_quadrature_init() FEM_nQuadrature(2,1) = 1 allocate(FEM_quadrature_weights(2,1)%p(FEM_nQuadrature(2,1))) - FEM_quadrature_weights(2,1)%p(1) = 1._pReal + FEM_quadrature_weights(2,1)%p(1) = 1._pReal FEM_quadrature_points (2,1)%p = permutationStar3([1._pReal/3._pReal]) @@ -95,11 +95,11 @@ subroutine FEM_quadrature_init() FEM_nQuadrature(2,5) = 16 allocate(FEM_quadrature_weights(2,5)%p(FEM_nQuadrature(2,5))) - FEM_quadrature_weights(2,5)%p(1 ) = 1.4431560767778717e-1_pReal - FEM_quadrature_weights(2,5)%p(2:4) = 9.509163426728463e-2_pReal - FEM_quadrature_weights(2,5)%p(5:7) = 1.0321737053471825e-1_pReal - FEM_quadrature_weights(2,5)%p(8:10) = 3.2458497623198080e-2_pReal - FEM_quadrature_weights(2,5)%p(11:16)= 2.7230314174434994e-2_pReal + FEM_quadrature_weights(2,5)%p(1:1) = 1.4431560767778717e-1_pReal + FEM_quadrature_weights(2,5)%p(2:4) = 9.509163426728463e-2_pReal + FEM_quadrature_weights(2,5)%p(5:7) = 1.0321737053471825e-1_pReal + FEM_quadrature_weights(2,5)%p(8:10) = 3.2458497623198080e-2_pReal + FEM_quadrature_weights(2,5)%p(11:16) = 2.7230314174434994e-2_pReal FEM_quadrature_points (2,5)%p = [ & permutationStar3([1._pReal/3._pReal]), & @@ -193,7 +193,9 @@ pure function permutationStar3(point) result(qPt) real(pReal), dimension(3,1) :: temp - temp(:,1) = [point(1), point(1), point(1)] + + temp = reshape([ & + point(1), point(1), point(1)],shape(temp)) qPt = reshape(matmul(triangle, temp),[2]) @@ -210,9 +212,11 @@ pure function permutationStar21(point) result(qPt) real(pReal), dimension(3,3) :: temp - temp(:,1) = [point(1), point(1), 1.0_pReal - 2.0_pReal*point(1)] - temp(:,2) = [point(1), 1.0_pReal - 2.0_pReal*point(1), point(1)] - temp(:,3) = [1.0_pReal - 2.0_pReal*point(1), point(1), point(1)] + + temp = reshape([ & + point(1), point(1), 1.0_pReal - 2.0_pReal*point(1), & + point(1), 1.0_pReal - 2.0_pReal*point(1), point(1), & + 1.0_pReal - 2.0_pReal*point(1), point(1), point(1)],shape(temp)) qPt = reshape(matmul(triangle, temp),[6]) @@ -229,12 +233,14 @@ pure function permutationStar111(point) result(qPt) real(pReal), dimension(3,6) :: temp - temp(:,1) = [point(1), point(2), 1.0_pReal - point(1) - point(2)] - temp(:,2) = [point(1), 1.0_pReal - point(1) - point(2), point(2)] - temp(:,3) = [point(2), point(1), 1.0_pReal - point(1) - point(2)] - temp(:,4) = [point(2), 1.0_pReal - point(1) - point(2), point(1)] - temp(:,5) = [1.0_pReal - point(1) - point(2), point(2), point(1)] - temp(:,6) = [1.0_pReal - point(1) - point(2), point(1), point(2)] + + temp = reshape([ & + point(1), point(2), 1.0_pReal - point(1) - point(2), & + point(1), 1.0_pReal - point(1) - point(2), point(2), & + point(2), point(1), 1.0_pReal - point(1) - point(2), & + point(2), 1.0_pReal - point(1) - point(2), point(1), & + 1.0_pReal - point(1) - point(2), point(2), point(1), & + 1.0_pReal - point(1) - point(2), point(1), point(2)],shape(temp)) qPt = reshape(matmul(triangle, temp),[12]) @@ -251,7 +257,9 @@ pure function permutationStar4(point) result(qPt) real(pReal), dimension(4,1) :: temp - temp(:,1) = [point(1), point(1), point(1), point(1)] + + temp = reshape([ & + point(1), point(1), point(1), point(1)],shape(temp)) qPt = reshape(matmul(tetrahedron, temp),[3]) @@ -268,10 +276,12 @@ pure function permutationStar31(point) result(qPt) real(pReal), dimension(4,4) :: temp - temp(:,1) = [point(1), point(1), point(1), 1.0_pReal - 3.0_pReal*point(1)] - temp(:,2) = [point(1), point(1), 1.0_pReal - 3.0_pReal*point(1), point(1)] - temp(:,3) = [point(1), 1.0_pReal - 3.0_pReal*point(1), point(1), point(1)] - temp(:,4) = [1.0_pReal - 3.0_pReal*point(1), point(1), point(1), point(1)] + + temp = reshape([ & + point(1), point(1), point(1), 1.0_pReal - 3.0_pReal*point(1), & + point(1), point(1), 1.0_pReal - 3.0_pReal*point(1), point(1), & + point(1), 1.0_pReal - 3.0_pReal*point(1), point(1), point(1), & + 1.0_pReal - 3.0_pReal*point(1), point(1), point(1), point(1)],shape(temp)) qPt = reshape(matmul(tetrahedron, temp),[12]) @@ -281,19 +291,21 @@ end function permutationStar31 !-------------------------------------------------------------------------------------------------- !> @brief star 22 permutation of input !-------------------------------------------------------------------------------------------------- -pure function permutationStar22(point) result(qPt) +function permutationStar22(point) result(qPt) real(pReal), dimension(18) :: qPt real(pReal), dimension(1), intent(in) :: point real(pReal), dimension(4,6) :: temp - temp(:,1) = [point(1), point(1), 0.5_pReal - point(1), 0.5_pReal - point(1)] - temp(:,2) = [point(1), 0.5_pReal - point(1), point(1), 0.5_pReal - point(1)] - temp(:,3) = [0.5_pReal - point(1), point(1), point(1), 0.5_pReal - point(1)] - temp(:,4) = [0.5_pReal - point(1), point(1), 0.5_pReal - point(1), point(1)] - temp(:,5) = [0.5_pReal - point(1), 0.5_pReal - point(1), point(1), point(1)] - temp(:,6) = [point(1), 0.5_pReal - point(1), 0.5_pReal - point(1), point(1)] + + temp = reshape([ & + point(1), point(1), 0.5_pReal - point(1), 0.5_pReal - point(1), & + point(1), 0.5_pReal - point(1), point(1), 0.5_pReal - point(1), & + 0.5_pReal - point(1), point(1), point(1), 0.5_pReal - point(1), & + 0.5_pReal - point(1), point(1), 0.5_pReal - point(1), point(1), & + 0.5_pReal - point(1), 0.5_pReal - point(1), point(1), point(1), & + point(1), 0.5_pReal - point(1), 0.5_pReal - point(1), point(1)],shape(temp)) qPt = reshape(matmul(tetrahedron, temp),[18]) @@ -310,18 +322,20 @@ pure function permutationStar211(point) result(qPt) real(pReal), dimension(4,12) :: temp - temp(:,1 ) = [point(1), point(1), point(2), 1.0_pReal - 2.0_pReal*point(1) - point(2)] - temp(:,2 ) = [point(1), point(1), 1.0_pReal - 2.0_pReal*point(1) - point(2), point(2)] - temp(:,3 ) = [point(1), point(2), point(1), 1.0_pReal - 2.0_pReal*point(1) - point(2)] - temp(:,4 ) = [point(1), point(2), 1.0_pReal - 2.0_pReal*point(1) - point(2), point(1)] - temp(:,5 ) = [point(1), 1.0_pReal - 2.0_pReal*point(1) - point(2), point(1), point(2)] - temp(:,6 ) = [point(1), 1.0_pReal - 2.0_pReal*point(1) - point(2), point(2), point(1)] - temp(:,7 ) = [point(2), point(1), point(1), 1.0_pReal - 2.0_pReal*point(1) - point(2)] - temp(:,8 ) = [point(2), point(1), 1.0_pReal - 2.0_pReal*point(1) - point(2), point(1)] - temp(:,9 ) = [point(2), 1.0_pReal - 2.0_pReal*point(1) - point(2), point(1), point(1)] - temp(:,10) = [1.0_pReal - 2.0_pReal*point(1) - point(2), point(1), point(1), point(2)] - temp(:,11) = [1.0_pReal - 2.0_pReal*point(1) - point(2), point(1), point(2), point(1)] - temp(:,12) = [1.0_pReal - 2.0_pReal*point(1) - point(2), point(2), point(1), point(1)] + + temp = reshape([ & + point(1), point(1), point(2), 1.0_pReal - 2.0_pReal*point(1) - point(2), & + point(1), point(1), 1.0_pReal - 2.0_pReal*point(1) - point(2), point(2), & + point(1), point(2), point(1), 1.0_pReal - 2.0_pReal*point(1) - point(2), & + point(1), point(2), 1.0_pReal - 2.0_pReal*point(1) - point(2), point(1), & + point(1), 1.0_pReal - 2.0_pReal*point(1) - point(2), point(1), point(2), & + point(1), 1.0_pReal - 2.0_pReal*point(1) - point(2), point(2), point(1), & + point(2), point(1), point(1), 1.0_pReal - 2.0_pReal*point(1) - point(2), & + point(2), point(1), 1.0_pReal - 2.0_pReal*point(1) - point(2), point(1), & + point(2), 1.0_pReal - 2.0_pReal*point(1) - point(2), point(1), point(1), & + 1.0_pReal - 2.0_pReal*point(1) - point(2), point(1), point(1), point(2), & + 1.0_pReal - 2.0_pReal*point(1) - point(2), point(1), point(2), point(1), & + 1.0_pReal - 2.0_pReal*point(1) - point(2), point(2), point(1), point(1)],shape(temp)) qPt = reshape(matmul(tetrahedron, temp),[36]) @@ -338,30 +352,32 @@ pure function permutationStar1111(point) result(qPt) real(pReal), dimension(4,24) :: temp - temp(:,1 ) = [point(1), point(2), point(3), 1.0_pReal - point(1) - point(2)- point(3)] - temp(:,2 ) = [point(1), point(2), 1.0_pReal - point(1) - point(2)- point(3), point(3)] - temp(:,3 ) = [point(1), point(3), point(2), 1.0_pReal - point(1) - point(2)- point(3)] - temp(:,4 ) = [point(1), point(3), 1.0_pReal - point(1) - point(2)- point(3), point(2)] - temp(:,5 ) = [point(1), 1.0_pReal - point(1) - point(2)- point(3), point(2), point(3)] - temp(:,6 ) = [point(1), 1.0_pReal - point(1) - point(2)- point(3), point(3), point(2)] - temp(:,7 ) = [point(2), point(1), point(3), 1.0_pReal - point(1) - point(2)- point(3)] - temp(:,8 ) = [point(2), point(1), 1.0_pReal - point(1) - point(2)- point(3), point(3)] - temp(:,9 ) = [point(2), point(3), point(1), 1.0_pReal - point(1) - point(2)- point(3)] - temp(:,10) = [point(2), point(3), 1.0_pReal - point(1) - point(2)- point(3), point(1)] - temp(:,11) = [point(2), 1.0_pReal - point(1) - point(2)- point(3), point(1), point(3)] - temp(:,12) = [point(2), 1.0_pReal - point(1) - point(2)- point(3), point(3), point(1)] - temp(:,13) = [point(3), point(1), point(2), 1.0_pReal - point(1) - point(2)- point(3)] - temp(:,14) = [point(3), point(1), 1.0_pReal - point(1) - point(2)- point(3), point(2)] - temp(:,15) = [point(3), point(2), point(1), 1.0_pReal - point(1) - point(2)- point(3)] - temp(:,16) = [point(3), point(2), 1.0_pReal - point(1) - point(2)- point(3), point(1)] - temp(:,17) = [point(3), 1.0_pReal - point(1) - point(2)- point(3), point(1), point(2)] - temp(:,18) = [point(3), 1.0_pReal - point(1) - point(2)- point(3), point(2), point(1)] - temp(:,19) = [1.0_pReal - point(1) - point(2)- point(3), point(1), point(2), point(3)] - temp(:,20) = [1.0_pReal - point(1) - point(2)- point(3), point(1), point(3), point(2)] - temp(:,21) = [1.0_pReal - point(1) - point(2)- point(3), point(2), point(1), point(3)] - temp(:,22) = [1.0_pReal - point(1) - point(2)- point(3), point(2), point(3), point(1)] - temp(:,23) = [1.0_pReal - point(1) - point(2)- point(3), point(3), point(1), point(2)] - temp(:,24) = [1.0_pReal - point(1) - point(2)- point(3), point(3), point(2), point(1)] + + temp = reshape([ & + point(1), point(2), point(3), 1.0_pReal - point(1) - point(2)- point(3), & + point(1), point(2), 1.0_pReal - point(1) - point(2)- point(3), point(3), & + point(1), point(3), point(2), 1.0_pReal - point(1) - point(2)- point(3), & + point(1), point(3), 1.0_pReal - point(1) - point(2)- point(3), point(2), & + point(1), 1.0_pReal - point(1) - point(2)- point(3), point(2), point(3), & + point(1), 1.0_pReal - point(1) - point(2)- point(3), point(3), point(2), & + point(2), point(1), point(3), 1.0_pReal - point(1) - point(2)- point(3), & + point(2), point(1), 1.0_pReal - point(1) - point(2)- point(3), point(3), & + point(2), point(3), point(1), 1.0_pReal - point(1) - point(2)- point(3), & + point(2), point(3), 1.0_pReal - point(1) - point(2)- point(3), point(1), & + point(2), 1.0_pReal - point(1) - point(2)- point(3), point(1), point(3), & + point(2), 1.0_pReal - point(1) - point(2)- point(3), point(3), point(1), & + point(3), point(1), point(2), 1.0_pReal - point(1) - point(2)- point(3), & + point(3), point(1), 1.0_pReal - point(1) - point(2)- point(3), point(2), & + point(3), point(2), point(1), 1.0_pReal - point(1) - point(2)- point(3), & + point(3), point(2), 1.0_pReal - point(1) - point(2)- point(3), point(1), & + point(3), 1.0_pReal - point(1) - point(2)- point(3), point(1), point(2), & + point(3), 1.0_pReal - point(1) - point(2)- point(3), point(2), point(1), & + 1.0_pReal - point(1) - point(2)- point(3), point(1), point(2), point(3), & + 1.0_pReal - point(1) - point(2)- point(3), point(1), point(3), point(2), & + 1.0_pReal - point(1) - point(2)- point(3), point(2), point(1), point(3), & + 1.0_pReal - point(1) - point(2)- point(3), point(2), point(3), point(1), & + 1.0_pReal - point(1) - point(2)- point(3), point(3), point(1), point(2), & + 1.0_pReal - point(1) - point(2)- point(3), point(3), point(2), point(1)],shape(temp)) qPt = reshape(matmul(tetrahedron, temp),[72]) @@ -369,13 +385,14 @@ end function permutationStar1111 !-------------------------------------------------------------------------------------------------- -!> @brief Check correctness of quadrature weights. +!> @brief Check correctness of quadrature weights and points. !-------------------------------------------------------------------------------------------------- subroutine selfTest integer :: o, d, n real(pReal), dimension(2:3), parameter :: w = [3.0_pReal,2.0_pReal] + do d = lbound(FEM_quadrature_weights,1), ubound(FEM_quadrature_weights,1) do o = lbound(FEM_quadrature_weights(d,:),1), ubound(FEM_quadrature_weights(d,:),1) if (dNeq(sum(FEM_quadrature_weights(d,o)%p),1.0_pReal,5e-15_pReal)) & From 58a79219666525ae8be6d86051e190da983f3573 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 21 Jul 2021 22:45:04 +0200 Subject: [PATCH 16/63] no need for temp, pack instead of reshape for 1D array --- src/mesh/FEM_quadrature.f90 | 68 ++++++++++--------------------------- 1 file changed, 18 insertions(+), 50 deletions(-) diff --git a/src/mesh/FEM_quadrature.f90 b/src/mesh/FEM_quadrature.f90 index 278bfb5db..bf1f1567f 100644 --- a/src/mesh/FEM_quadrature.f90 +++ b/src/mesh/FEM_quadrature.f90 @@ -20,12 +20,12 @@ module FEM_quadrature -1.0_pReal, 1.0_pReal, -1.0_pReal, & -1.0_pReal, -1.0_pReal, 1.0_pReal], shape=[3,4]) - type :: group_float !< variable length datatype used for storage of state + type :: group_float !< variable length datatype real(pReal), dimension(:), allocatable :: p end type group_float integer, dimension(2:3,maxOrder), public, protected :: & - FEM_nQuadrature !< number of quadrature points for a given spatial dimension(2-3) and interpolation order(1-maxOrder) + FEM_nQuadrature !< number of quadrature points for spatial dimension(2-3) and interpolation order (1-maxOrder) type(group_float), dimension(2:3,maxOrder), public, protected :: & FEM_quadrature_weights, & !< quadrature weights for each quadrature rule FEM_quadrature_points !< quadrature point coordinates (in simplical system) for each quadrature rule @@ -191,13 +191,9 @@ pure function permutationStar3(point) result(qPt) real(pReal), dimension(2) :: qPt real(pReal), dimension(1), intent(in) :: point - real(pReal), dimension(3,1) :: temp - - temp = reshape([ & - point(1), point(1), point(1)],shape(temp)) - - qPt = reshape(matmul(triangle, temp),[2]) + qPt = pack(matmul(triangle,reshape([ & + point(1), point(1), point(1)],[3,1])),.true.) end function permutationStar3 @@ -210,15 +206,11 @@ pure function permutationStar21(point) result(qPt) real(pReal), dimension(6) :: qPt real(pReal), dimension(1), intent(in) :: point - real(pReal), dimension(3,3) :: temp - - temp = reshape([ & + qPt = pack(matmul(triangle,reshape([ & point(1), point(1), 1.0_pReal - 2.0_pReal*point(1), & point(1), 1.0_pReal - 2.0_pReal*point(1), point(1), & - 1.0_pReal - 2.0_pReal*point(1), point(1), point(1)],shape(temp)) - - qPt = reshape(matmul(triangle, temp),[6]) + 1.0_pReal - 2.0_pReal*point(1), point(1), point(1)],[3,3])),.true.) end function permutationStar21 @@ -231,18 +223,14 @@ pure function permutationStar111(point) result(qPt) real(pReal), dimension(12) :: qPt real(pReal), dimension(2), intent(in) :: point - real(pReal), dimension(3,6) :: temp - - temp = reshape([ & + qPt = pack(matmul(triangle,reshape([ & point(1), point(2), 1.0_pReal - point(1) - point(2), & point(1), 1.0_pReal - point(1) - point(2), point(2), & point(2), point(1), 1.0_pReal - point(1) - point(2), & point(2), 1.0_pReal - point(1) - point(2), point(1), & 1.0_pReal - point(1) - point(2), point(2), point(1), & - 1.0_pReal - point(1) - point(2), point(1), point(2)],shape(temp)) - - qPt = reshape(matmul(triangle, temp),[12]) + 1.0_pReal - point(1) - point(2), point(1), point(2)],[3,6])),.true.) end function permutationStar111 @@ -255,13 +243,9 @@ pure function permutationStar4(point) result(qPt) real(pReal), dimension(3) :: qPt real(pReal), dimension(1), intent(in) :: point - real(pReal), dimension(4,1) :: temp - - temp = reshape([ & - point(1), point(1), point(1), point(1)],shape(temp)) - - qPt = reshape(matmul(tetrahedron, temp),[3]) + qPt = pack(matmul(tetrahedron,reshape([ & + point(1), point(1), point(1), point(1)],[4,1])),.true.) end function permutationStar4 @@ -274,16 +258,12 @@ pure function permutationStar31(point) result(qPt) real(pReal), dimension(12) :: qPt real(pReal), dimension(1), intent(in) :: point - real(pReal), dimension(4,4) :: temp - - temp = reshape([ & + qPt = pack(matmul(tetrahedron,reshape([ & point(1), point(1), point(1), 1.0_pReal - 3.0_pReal*point(1), & point(1), point(1), 1.0_pReal - 3.0_pReal*point(1), point(1), & point(1), 1.0_pReal - 3.0_pReal*point(1), point(1), point(1), & - 1.0_pReal - 3.0_pReal*point(1), point(1), point(1), point(1)],shape(temp)) - - qPt = reshape(matmul(tetrahedron, temp),[12]) + 1.0_pReal - 3.0_pReal*point(1), point(1), point(1), point(1)],[4,4])),.true.) end function permutationStar31 @@ -296,18 +276,14 @@ function permutationStar22(point) result(qPt) real(pReal), dimension(18) :: qPt real(pReal), dimension(1), intent(in) :: point - real(pReal), dimension(4,6) :: temp - - temp = reshape([ & + qPt = pack(matmul(tetrahedron,reshape([ & point(1), point(1), 0.5_pReal - point(1), 0.5_pReal - point(1), & point(1), 0.5_pReal - point(1), point(1), 0.5_pReal - point(1), & 0.5_pReal - point(1), point(1), point(1), 0.5_pReal - point(1), & 0.5_pReal - point(1), point(1), 0.5_pReal - point(1), point(1), & 0.5_pReal - point(1), 0.5_pReal - point(1), point(1), point(1), & - point(1), 0.5_pReal - point(1), 0.5_pReal - point(1), point(1)],shape(temp)) - - qPt = reshape(matmul(tetrahedron, temp),[18]) + point(1), 0.5_pReal - point(1), 0.5_pReal - point(1), point(1)],[4,6])),.true.) end function permutationStar22 @@ -320,10 +296,8 @@ pure function permutationStar211(point) result(qPt) real(pReal), dimension(36) :: qPt real(pReal), dimension(2), intent(in) :: point - real(pReal), dimension(4,12) :: temp - - temp = reshape([ & + qPt = pack(matmul(tetrahedron,reshape([ & point(1), point(1), point(2), 1.0_pReal - 2.0_pReal*point(1) - point(2), & point(1), point(1), 1.0_pReal - 2.0_pReal*point(1) - point(2), point(2), & point(1), point(2), point(1), 1.0_pReal - 2.0_pReal*point(1) - point(2), & @@ -335,9 +309,7 @@ pure function permutationStar211(point) result(qPt) point(2), 1.0_pReal - 2.0_pReal*point(1) - point(2), point(1), point(1), & 1.0_pReal - 2.0_pReal*point(1) - point(2), point(1), point(1), point(2), & 1.0_pReal - 2.0_pReal*point(1) - point(2), point(1), point(2), point(1), & - 1.0_pReal - 2.0_pReal*point(1) - point(2), point(2), point(1), point(1)],shape(temp)) - - qPt = reshape(matmul(tetrahedron, temp),[36]) + 1.0_pReal - 2.0_pReal*point(1) - point(2), point(2), point(1), point(1)],[4,12])),.true.) end function permutationStar211 @@ -350,10 +322,8 @@ pure function permutationStar1111(point) result(qPt) real(pReal), dimension(72) :: qPt real(pReal), dimension(3), intent(in) :: point - real(pReal), dimension(4,24) :: temp - - temp = reshape([ & + qPt = pack(matmul(tetrahedron,reshape([ & point(1), point(2), point(3), 1.0_pReal - point(1) - point(2)- point(3), & point(1), point(2), 1.0_pReal - point(1) - point(2)- point(3), point(3), & point(1), point(3), point(2), 1.0_pReal - point(1) - point(2)- point(3), & @@ -377,9 +347,7 @@ pure function permutationStar1111(point) result(qPt) 1.0_pReal - point(1) - point(2)- point(3), point(2), point(1), point(3), & 1.0_pReal - point(1) - point(2)- point(3), point(2), point(3), point(1), & 1.0_pReal - point(1) - point(2)- point(3), point(3), point(1), point(2), & - 1.0_pReal - point(1) - point(2)- point(3), point(3), point(2), point(1)],shape(temp)) - - qPt = reshape(matmul(tetrahedron, temp),[72]) + 1.0_pReal - point(1) - point(2)- point(3), point(3), point(2), point(1)],[4,24])),.true.) end function permutationStar1111 From de7ef430958183843867ea2b42430f16919af16f Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 21 Jul 2021 23:43:55 +0200 Subject: [PATCH 17/63] documenting/silencing gfortran --- src/mesh/FEM_quadrature.f90 | 6 +++--- src/mesh/discretization_mesh.f90 | 6 ++---- src/quit.f90 | 5 ++++- 3 files changed, 9 insertions(+), 8 deletions(-) diff --git a/src/mesh/FEM_quadrature.f90 b/src/mesh/FEM_quadrature.f90 index bf1f1567f..0c470ca6b 100644 --- a/src/mesh/FEM_quadrature.f90 +++ b/src/mesh/FEM_quadrature.f90 @@ -124,7 +124,7 @@ subroutine FEM_quadrature_init() allocate(FEM_quadrature_weights(3,2)%p(FEM_nQuadrature(3,2))) FEM_quadrature_weights(3,2)%p(1:4) = 0.25_pReal - FEM_quadrature_points (3,2)%p = permutationStar31([0.1381966011250105_pReal]) + FEM_quadrature_points (3,2)%p = permutationStar31([1.3819660112501052e-1_pReal]) !-------------------------------------------------------------------------------------------------- ! 3D cubic @@ -141,7 +141,7 @@ subroutine FEM_quadrature_init() permutationStar22([4.5503704125649649e-2_pReal]) ] !-------------------------------------------------------------------------------------------------- -! 3D quartic +! 3D quartic (lower precision/unknown source) FEM_nQuadrature(3,4) = 35 allocate(FEM_quadrature_weights(3,4)%p(FEM_nQuadrature(3,4))) @@ -159,7 +159,7 @@ subroutine FEM_quadrature_init() permutationStar4([0.25_pReal]) ] !-------------------------------------------------------------------------------------------------- -! 3D quintic +! 3D quintic (lower precision/unknown source) FEM_nQuadrature(3,5) = 56 allocate(FEM_quadrature_weights(3,5)%p(FEM_nQuadrature(3,5))) diff --git a/src/mesh/discretization_mesh.f90 b/src/mesh/discretization_mesh.f90 index 49eaf0e18..422f4ca85 100644 --- a/src/mesh/discretization_mesh.f90 +++ b/src/mesh/discretization_mesh.f90 @@ -71,16 +71,14 @@ subroutine discretization_mesh_init(restart) logical, intent(in) :: restart - integer, allocatable, dimension(:) :: chunkPos integer :: dimPlex, & mesh_Nnodes, & !< total number of nodes in mesh - j, l, & + j, & debug_element, debug_ip PetscSF :: sf DM :: globalMesh PetscInt :: nFaceSets PetscInt, pointer, dimension(:) :: pFaceSets - character(len=pStringLen), dimension(:), allocatable :: fileContent IS :: faceSetIS PetscErrorCode :: ierr integer, dimension(:), allocatable :: & @@ -88,7 +86,7 @@ subroutine discretization_mesh_init(restart) class(tNode), pointer :: & num_mesh integer :: integrationOrder !< order of quadrature rule required - type(tvec) :: coords_node0 + type(tvec) :: coords_node0 print'(/,a)', ' <<<+- discretization_mesh init -+>>>' diff --git a/src/quit.f90 b/src/quit.f90 index 786241d85..6361ff783 100644 --- a/src/quit.f90 +++ b/src/quit.f90 @@ -6,7 +6,10 @@ !-------------------------------------------------------------------------------------------------- subroutine quit(stop_id) #include - use PetscSys + use PETScSys +#if (PETSC_VERSION_MAJOR==3 && PETSC_VERSION_MINOR>14) && !defined(PETSC_HAVE_MPI_F90MODULE_VISIBILITY) + use MPI_f08 +#endif use HDF5 implicit none From 9767866c1e8d944483fb40857458dda05cb55fab Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Thu, 22 Jul 2021 07:24:36 +0200 Subject: [PATCH 18/63] cleaned+updated test --- PRIVATE | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/PRIVATE b/PRIVATE index 567067daf..96cf7af03 160000 --- a/PRIVATE +++ b/PRIVATE @@ -1 +1 @@ -Subproject commit 567067daf05866c8e022fc0af6f441840f143b81 +Subproject commit 96cf7af039b3153f91c037bcce62eb456af51f4b From 173a5f8e55e365a6f875096686953b6d8e1db717 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Thu, 22 Jul 2021 07:48:01 +0200 Subject: [PATCH 19/63] less public/unneeded variables --- src/mesh/DAMASK_mesh.f90 | 1 + src/mesh/FEM_utilities.f90 | 17 ++++------------- src/mesh/discretization_mesh.f90 | 10 +++++----- 3 files changed, 10 insertions(+), 18 deletions(-) diff --git a/src/mesh/DAMASK_mesh.f90 b/src/mesh/DAMASK_mesh.f90 index 342dec5c0..2eb6a7d91 100644 --- a/src/mesh/DAMASK_mesh.f90 +++ b/src/mesh/DAMASK_mesh.f90 @@ -21,6 +21,7 @@ program DAMASK_mesh use mesh_mechanical_FEM implicit none + integer :: nActiveFields = 0 !-------------------------------------------------------------------------------------------------- ! variables related to information from load case and geom file diff --git a/src/mesh/FEM_utilities.f90 b/src/mesh/FEM_utilities.f90 index 5bb6d5f44..1c073dd69 100644 --- a/src/mesh/FEM_utilities.f90 +++ b/src/mesh/FEM_utilities.f90 @@ -23,14 +23,8 @@ module FEM_utilities implicit none private -!-------------------------------------------------------------------------------------------------- - logical, public :: cutBack = .false. !< cut back of BVP solver in case convergence is not achieved or a material point is terminally ill - integer, public, parameter :: maxFields = 6 - integer, public :: nActiveFields = 0 - -!-------------------------------------------------------------------------------------------------- -! grid related information information - real(pReal), public :: wgt !< weighting factor 1/Nelems + logical, public :: cutBack = .false. !< cut back of BVP solver in case convergence is not achieved or a material point is terminally ill + real(pReal), public, protected :: wgt !< weighting factor 1/Nelems !-------------------------------------------------------------------------------------------------- @@ -49,10 +43,6 @@ module FEM_utilities COMPONENT_MECH_Z_ID end enum -!-------------------------------------------------------------------------------------------------- -! variables controlling debugging - logical :: & - debugPETSc !< use some in debug defined options for more verbose PETSc solution !-------------------------------------------------------------------------------------------------- ! derived types @@ -109,8 +99,9 @@ subroutine FEM_utilities_init integer :: structOrder !< order of displacement shape functions character(len=*), parameter :: & PETSCDEBUG = ' -snes_view -snes_monitor ' - PetscErrorCode :: ierr + logical :: debugPETSc !< use some in debug defined options for more verbose PETSc solution + print'(/,a)', ' <<<+- FEM_utilities init -+>>>' diff --git a/src/mesh/discretization_mesh.f90 b/src/mesh/discretization_mesh.f90 index 422f4ca85..e5f989484 100644 --- a/src/mesh/discretization_mesh.f90 +++ b/src/mesh/discretization_mesh.f90 @@ -40,6 +40,11 @@ module discretization_mesh mesh_maxNips !< max number of IPs in any CP element !!!! BEGIN DEPRECATED !!!!! + DM, public :: geomMesh + + PetscInt, dimension(:), allocatable, public, protected :: & + mesh_boundaries + real(pReal), dimension(:,:), allocatable :: & mesh_ipVolume, & !< volume associated with IP (initially!) mesh_node0 !< node x,y,z coordinates (initially!) @@ -50,11 +55,6 @@ module discretization_mesh real(pReal), dimension(:,:,:), allocatable :: & mesh_ipCoordinates !< IP x,y,z coordinates (after deformation!) - DM, public :: geomMesh - - PetscInt, dimension(:), allocatable, public, protected :: & - mesh_boundaries - public :: & discretization_mesh_init, & mesh_FEM_build_ipVolumes, & From eb834b635dee5e301dcea4d0d8fef29e01a38b79 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Thu, 22 Jul 2021 07:58:34 +0200 Subject: [PATCH 20/63] unify notation with grid consistent style: Symbols, not descriptions. also removed untested/unused loginc functionality. Once load cases are written in YAML, we can introduce the scaling as in DAMASK_grid --- PRIVATE | 2 +- examples/mesh/tensionZ.load | 4 +-- src/mesh/DAMASK_mesh.f90 | 70 +++++++++++++++---------------------- src/mesh/FEM_utilities.f90 | 10 ------ 4 files changed, 32 insertions(+), 54 deletions(-) diff --git a/PRIVATE b/PRIVATE index 96cf7af03..d92b030e5 160000 --- a/PRIVATE +++ b/PRIVATE @@ -1 +1 @@ -Subproject commit 96cf7af039b3153f91c037bcce62eb456af51f4b +Subproject commit d92b030e5777e718c77edc2e1e93abfa0981b024 diff --git a/examples/mesh/tensionZ.load b/examples/mesh/tensionZ.load index eb9a7c426..b3588daef 100644 --- a/examples/mesh/tensionZ.load +++ b/examples/mesh/tensionZ.load @@ -1,11 +1,11 @@ #initial elastic step -$Loadcase 1 time 0.0005 incs 1 frequency 5 +$Loadcase 1 t 0.0005 N 1 f_out 5 Face 1 X 0.01 Face 2 X 0.0 Face 2 Y 0.0 Face 2 Z 0.0 $EndLoadcase -$Loadcase 2 time 10.0 incs 200 frequency 5 +$Loadcase 2 t 10.0 N 200 f_out 5 Face 1 X 0.01 Face 2 X 0.0 Face 2 Y 0.0 diff --git a/src/mesh/DAMASK_mesh.f90 b/src/mesh/DAMASK_mesh.f90 index 2eb6a7d91..d1568e87f 100644 --- a/src/mesh/DAMASK_mesh.f90 +++ b/src/mesh/DAMASK_mesh.f90 @@ -23,6 +23,15 @@ program DAMASK_mesh implicit none integer :: nActiveFields = 0 + 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 + 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 @@ -104,8 +113,8 @@ program DAMASK_mesh chunkPos = IO_stringPos(line) do i = 1, chunkPos(1) ! reading compulsory parameters for loadcase - select case (IO_lc(IO_stringValue(line,chunkPos,i))) - case('$loadcase') + select case (IO_stringValue(line,chunkPos,i)) + case('$Loadcase') N_def = N_def + 1 end select enddo ! count all identifiers to allocate memory and do sanity check @@ -152,39 +161,36 @@ program DAMASK_mesh chunkPos = IO_stringPos(line) do i = 1, chunkPos(1) - select case (IO_lc(IO_stringValue(line,chunkPos,i))) + select case (IO_stringValue(line,chunkPos,i)) !-------------------------------------------------------------------------------------------------- ! loadcase information - case('$loadcase') + case('$Loadcase') currentLoadCase = IO_intValue(line,chunkPos,i+1) - case('face') + case('Face') currentFace = IO_intValue(line,chunkPos,i+1) currentFaceSet = -1 do faceSet = 1, mesh_Nboundaries if (mesh_boundaries(faceSet) == currentFace) currentFaceSet = faceSet enddo if (currentFaceSet < 0) call IO_error(error_ID = 837, ext_msg = 'invalid BC') - case('t','time','delta') ! increment time + case('t') loadCases(currentLoadCase)%time = IO_floatValue(line,chunkPos,i+1) - case('n','incs','increments','steps') ! number of increments + case('N') loadCases(currentLoadCase)%incs = IO_intValue(line,chunkPos,i+1) - case('logincs','logincrements','logsteps') ! number of increments (switch to log time scaling) - loadCases(currentLoadCase)%incs = IO_intValue(line,chunkPos,i+1) - loadCases(currentLoadCase)%logscale = 1 - case('freq','frequency','outputfreq') ! frequency of result writings + case('f_out') loadCases(currentLoadCase)%outputfrequency = IO_intValue(line,chunkPos,i+1) - case('guessreset','dropguessing') + 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_lc(IO_stringValue(line,chunkPos,i))) - case('x') + case('X','Y','Z') + select case(IO_stringValue(line,chunkPos,i)) + case('X') ID = COMPONENT_MECH_X_ID - case('y') + case('Y') ID = COMPONENT_MECH_Y_ID - case('z') + case('Z') ID = COMPONENT_MECH_Z_ID end select @@ -223,8 +229,8 @@ program DAMASK_mesh do component = 1, loadCases(currentLoadCase)%fieldBC(field)%nComponents if (loadCases(currentLoadCase)%fieldBC(field)%componentBC(component)%Mask(faceSet)) & print'(a,i2,a,i2,a,f12.7)', ' Face ', mesh_boundaries(faceSet), & - ' Component ', component, & - ' Value ', loadCases(currentLoadCase)%fieldBC(field)% & + ' Component ', component, & + ' Value ', loadCases(currentLoadCase)%fieldBC(field)% & componentBC(component)%Value(faceSet) enddo enddo @@ -267,32 +273,14 @@ program DAMASK_mesh !-------------------------------------------------------------------------------------------------- ! forwarding time timeIncOld = timeinc ! last timeinc that brought former inc to an end - if (loadCases(currentLoadCase)%logscale == 0) then ! linear scale - timeinc = loadCases(currentLoadCase)%time/real(loadCases(currentLoadCase)%incs,pReal) - else - if (currentLoadCase == 1) then ! 1st load case of logarithmic scale - if (inc == 1) then ! 1st inc of 1st load case of logarithmic scale - timeinc = loadCases(1)%time*(2.0_pReal**real( 1-loadCases(1)%incs ,pReal)) ! assume 1st inc is equal to 2nd - else ! not-1st inc of 1st load case of logarithmic scale - timeinc = loadCases(1)%time*(2.0_pReal**real(inc-1-loadCases(1)%incs ,pReal)) - endif - else ! not-1st load case of logarithmic scale - timeinc = time0 * & - ( (1.0_pReal + loadCases(currentLoadCase)%time/time0 )**(real( inc,pReal)/& - real(loadCases(currentLoadCase)%incs ,pReal))& - -(1.0_pReal + loadCases(currentLoadCase)%time/time0 )**(real( inc-1 ,pReal)/& - real(loadCases(currentLoadCase)%incs ,pReal))) - endif - endif + 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 - - - stepFraction = 0 ! fraction scaled by stepFactor**cutLevel + 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 - stepFraction = stepFraction + 1 ! count step + time = time + timeinc ! forward target time + stepFraction = stepFraction + 1 ! count step !-------------------------------------------------------------------------------------------------- ! report begin of new step diff --git a/src/mesh/FEM_utilities.f90 b/src/mesh/FEM_utilities.f90 index 1c073dd69..319997ad6 100644 --- a/src/mesh/FEM_utilities.f90 +++ b/src/mesh/FEM_utilities.f90 @@ -64,16 +64,6 @@ module FEM_utilities type(tComponentBC), allocatable :: componentBC(:) end type tFieldBC - type, public :: tLoadCase - real(pReal) :: time = 0.0_pReal !< length of increment - integer :: incs = 0, & !< number of increments - outputfrequency = 1, & !< frequency of result writes - logscale = 0 !< linear/logarithmic time inc flag - logical :: followFormerTrajectory = .true. !< follow trajectory of former loadcase - integer, allocatable, dimension(:) :: faceID - type(tFieldBC), allocatable, dimension(:) :: fieldBC - end type tLoadCase - public :: & FEM_utilities_init, & utilities_constitutiveResponse, & From 5b66db8a394ffe34b303039fe8cd17b370e60a30 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Thu, 22 Jul 2021 10:12:43 +0200 Subject: [PATCH 21/63] only mechanics at the moment will be extended, but most likely differently --- src/mesh/DAMASK_mesh.f90 | 113 ++++++++++++------------------------- src/mesh/FEM_utilities.f90 | 8 +-- 2 files changed, 41 insertions(+), 80 deletions(-) diff --git a/src/mesh/DAMASK_mesh.f90 b/src/mesh/DAMASK_mesh.f90 index d1568e87f..a8446d9fa 100644 --- a/src/mesh/DAMASK_mesh.f90 +++ b/src/mesh/DAMASK_mesh.f90 @@ -21,7 +21,6 @@ program DAMASK_mesh use mesh_mechanical_FEM implicit none - integer :: nActiveFields = 0 type :: tLoadCase real(pReal) :: time = 0.0_pReal !< length of increment @@ -78,7 +77,7 @@ program DAMASK_mesh type(tLoadCase), allocatable, dimension(:) :: loadCases !< array of all load cases type(tSolutionState), allocatable, dimension(:) :: solres - PetscInt :: faceSet, currentFaceSet, field, dimPlex + PetscInt :: faceSet, currentFaceSet, dimPlex PetscErrorCode :: ierr integer(kind(COMPONENT_UNDEFINED_ID)) :: ID external :: & @@ -101,8 +100,7 @@ program DAMASK_mesh ! reading basic information from load case file and allocate data structure containing load cases call DMGetDimension(geomMesh,dimPlex,ierr) !< dimension of mesh (2D or 3D) CHKERRA(ierr) - nActiveFields = 1 - allocate(solres(nActiveFields)) + allocate(solres(1)) !-------------------------------------------------------------------------------------------------- ! reading basic information from load case file and allocate data structure containing load cases @@ -124,32 +122,26 @@ program DAMASK_mesh allocate(loadCases(N_def)) do i = 1, size(loadCases) - allocate(loadCases(i)%fieldBC(nActiveFields)) - field = 1 - loadCases(i)%fieldBC(field)%ID = FIELD_MECH_ID + allocate(loadCases(i)%fieldBC(1)) + loadCases(i)%fieldBC(1)%ID = FIELD_MECH_ID enddo do i = 1, size(loadCases) - do field = 1, nActiveFields - select case (loadCases(i)%fieldBC(field)%ID) - case(FIELD_MECH_ID) - loadCases(i)%fieldBC(field)%nComponents = dimPlex !< X, Y (, Z) displacements - allocate(loadCases(i)%fieldBC(field)%componentBC(loadCases(i)%fieldBC(field)%nComponents)) - do component = 1, loadCases(i)%fieldBC(field)%nComponents - select case (component) - case (1) - loadCases(i)%fieldBC(field)%componentBC(component)%ID = COMPONENT_MECH_X_ID - case (2) - loadCases(i)%fieldBC(field)%componentBC(component)%ID = COMPONENT_MECH_Y_ID - case (3) - loadCases(i)%fieldBC(field)%componentBC(component)%ID = COMPONENT_MECH_Z_ID - end select - enddo + 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 + select case (component) + case (1) + loadCases(i)%fieldBC(1)%componentBC(component)%ID = COMPONENT_MECH_X_ID + case (2) + loadCases(i)%fieldBC(1)%componentBC(component)%ID = COMPONENT_MECH_Y_ID + case (3) + loadCases(i)%fieldBC(1)%componentBC(component)%ID = COMPONENT_MECH_Z_ID end select - do component = 1, loadCases(i)%fieldBC(field)%nComponents - allocate(loadCases(i)%fieldBC(field)%componentBC(component)%Value(mesh_Nboundaries), source = 0.0_pReal) - allocate(loadCases(i)%fieldBC(field)%componentBC(component)%Mask (mesh_Nboundaries), source = .false.) - enddo + enddo + 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.) enddo enddo @@ -194,16 +186,12 @@ program DAMASK_mesh ID = COMPONENT_MECH_Z_ID end select - do field = 1, nActiveFields - if (loadCases(currentLoadCase)%fieldBC(field)%ID == FIELD_MECH_ID) then - do component = 1, loadcases(currentLoadCase)%fieldBC(field)%nComponents - if (loadCases(currentLoadCase)%fieldBC(field)%componentBC(component)%ID == ID) then - loadCases(currentLoadCase)%fieldBC(field)%componentBC(component)%Mask (currentFaceSet) = & - .true. - loadCases(currentLoadCase)%fieldBC(field)%componentBC(component)%Value(currentFaceSet) = & - IO_floatValue(line,chunkPos,i+1) - endif - enddo + 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_floatValue(line,chunkPos,i+1) endif enddo end select @@ -219,21 +207,16 @@ program DAMASK_mesh print'(a,i0)', ' load case: ', currentLoadCase if (.not. loadCases(currentLoadCase)%followFormerTrajectory) & print'(a)', ' drop guessing along trajectory' - do field = 1, nActiveFields - select case (loadCases(currentLoadCase)%fieldBC(field)%ID) - case(FIELD_MECH_ID) - print'(a)', ' Field '//trim(FIELD_MECH_label) + print'(a)', ' Field '//trim(FIELD_MECH_label) - end select - do faceSet = 1, mesh_Nboundaries - do component = 1, loadCases(currentLoadCase)%fieldBC(field)%nComponents - if (loadCases(currentLoadCase)%fieldBC(field)%componentBC(component)%Mask(faceSet)) & - print'(a,i2,a,i2,a,f12.7)', ' Face ', mesh_boundaries(faceSet), & - ' Component ', component, & - ' Value ', loadCases(currentLoadCase)%fieldBC(field)% & - componentBC(component)%Value(faceSet) - enddo - enddo + do faceSet = 1, mesh_Nboundaries + do component = 1, loadCases(currentLoadCase)%fieldBC(1)%nComponents + if (loadCases(currentLoadCase)%fieldBC(1)%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) + enddo enddo print'(a,f12.6)', ' time: ', loadCases(currentLoadCase)%time if (loadCases(currentLoadCase)%incs < 1) errorID = 835 ! non-positive incs count @@ -247,12 +230,7 @@ program DAMASK_mesh !-------------------------------------------------------------------------------------------------- ! doing initialization depending on active solvers call FEM_Utilities_init - do field = 1, nActiveFields - select case (loadCases(1)%fieldBC(field)%ID) - case(FIELD_MECH_ID) - call FEM_mechanical_init(loadCases(1)%fieldBC(field)) - end select - enddo + call FEM_mechanical_init(loadCases(1)%fieldBC(1)) if (worldrank == 0) then open(newunit=statUnit,file=trim(getSolverJobName())//'.sta',form='FORMATTED',status='REPLACE') @@ -295,33 +273,16 @@ program DAMASK_mesh '-',stepFraction, '/', subStepFactor**cutBackLevel flush(IO_STDOUT) -!-------------------------------------------------------------------------------------------------- -! forward fields - do field = 1, nActiveFields - select case (loadCases(currentLoadCase)%fieldBC(field)%ID) - case(FIELD_MECH_ID) - call FEM_mechanical_forward (& - guess,timeinc,timeIncOld,loadCases(currentLoadCase)%fieldBC(field)) - - end select - enddo + call FEM_mechanical_forward(guess,timeinc,timeIncOld,loadCases(currentLoadCase)%fieldBC(1)) !-------------------------------------------------------------------------------------------------- ! solve fields stagIter = 0 stagIterate = .true. do while (stagIterate) - do field = 1, nActiveFields - select case (loadCases(currentLoadCase)%fieldBC(field)%ID) - case(FIELD_MECH_ID) - solres(field) = FEM_mechanical_solution (& - incInfo,timeinc,timeIncOld,loadCases(currentLoadCase)%fieldBC(field)) + solres(1) = FEM_mechanical_solution(incInfo,timeinc,timeIncOld,loadCases(currentLoadCase)%fieldBC(1)) + if(.not. solres(1)%converged) exit - end select - - if(.not. solres(field)%converged) exit ! no solution found - - enddo stagIter = stagIter + 1 stagIterate = stagIter < stagItMax & .and. all(solres(:)%converged) & diff --git a/src/mesh/FEM_utilities.f90 b/src/mesh/FEM_utilities.f90 index 319997ad6..64bcc3896 100644 --- a/src/mesh/FEM_utilities.f90 +++ b/src/mesh/FEM_utilities.f90 @@ -53,15 +53,15 @@ module FEM_utilities end type tSolutionState type, public :: tComponentBC - integer(kind(COMPONENT_UNDEFINED_ID)) :: ID - real(pReal), allocatable, dimension(:) :: Value - logical, allocatable, dimension(:) :: Mask + integer(kind(COMPONENT_UNDEFINED_ID)) :: ID + real(pReal), allocatable, dimension(:) :: Value + logical, allocatable, dimension(:) :: Mask end type tComponentBC type, public :: tFieldBC integer(kind(FIELD_UNDEFINED_ID)) :: ID integer :: nComponents = 0 - type(tComponentBC), allocatable :: componentBC(:) + type(tComponentBC), allocatable, dimension(:) :: componentBC end type tFieldBC public :: & From aa03e5930632147cf984710d0fef6f41c68eb37d Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Thu, 22 Jul 2021 11:22:25 +0200 Subject: [PATCH 22/63] same statement, different syntax but silences gfortran runtime warning with automated LHS, compiling in DEBUG mode (tested for mesh solver) gives: home/m/DAMASK/src/phase_mechanical.f90:1010:43: runtime error: signed integer overflow: -9223372036854775808 - 1 cannot be represented in type 'integer(kind=8)' --- src/phase_mechanical.f90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/phase_mechanical.f90 b/src/phase_mechanical.f90 index 8b04a0009..3c14accb5 100644 --- a/src/phase_mechanical.f90 +++ b/src/phase_mechanical.f90 @@ -1007,7 +1007,7 @@ module function phase_mechanical_constitutive(Delta_t,co,ip,el) result(converged subLi0 = phase_mechanical_Li0(ph)%data(1:3,1:3,en) subLp0 = phase_mechanical_Lp0(ph)%data(1:3,1:3,en) - subState0 = plasticState(ph)%State0(:,en) + allocate(subState0,source=plasticState(ph)%State0(:,en)) subFp0 = phase_mechanical_Fp0(ph)%data(1:3,1:3,en) subFi0 = phase_mechanical_Fi0(ph)%data(1:3,1:3,en) subF0 = phase_mechanical_F0(ph)%data(1:3,1:3,en) From 69cc0b528bf8a02695a59ed315cc316c852a07fb Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Thu, 22 Jul 2021 11:38:03 +0200 Subject: [PATCH 23/63] do calculation for initialization --- src/mesh/mesh_mech_FEM.f90 | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/src/mesh/mesh_mech_FEM.f90 b/src/mesh/mesh_mech_FEM.f90 index 1468e3fcc..6bcf44c3a 100644 --- a/src/mesh/mesh_mech_FEM.f90 +++ b/src/mesh/mesh_mech_FEM.f90 @@ -109,7 +109,7 @@ subroutine FEM_mechanical_init(fieldBC) character(len=*), parameter :: prefix = 'mechFE_' PetscErrorCode :: ierr - + real(pReal), dimension(3,3) :: devNull class(tNode), pointer :: & num_mesh @@ -258,6 +258,7 @@ subroutine FEM_mechanical_init(fieldBC) call DMPlexVecSetClosure(mechanical_mesh,section,solution_local,cell,px_scal,5,ierr) CHKERRQ(ierr) enddo + call utilities_constitutiveResponse(0.0_pReal,devNull,.true.) end subroutine FEM_mechanical_init @@ -397,7 +398,7 @@ subroutine FEM_mechanical_formResidual(dm_local,xx_local,f_local,dummy,ierr) !-------------------------------------------------------------------------------------------------- ! evaluate constitutive response - call Utilities_constitutiveResponse(params%timeinc,P_av,ForwardData) + call utilities_constitutiveResponse(params%timeinc,P_av,ForwardData) call MPI_Allreduce(MPI_IN_PLACE,terminallyIll,1,MPI_LOGICAL,MPI_LOR,MPI_COMM_WORLD,ierr) ForwardData = .false. From 5ac592eb9e331de2988f97946237092882ad5795 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Thu, 22 Jul 2021 12:31:10 +0200 Subject: [PATCH 24/63] initialize L_(i/p), L_(i,p)0 --- src/mesh/mesh_mech_FEM.f90 | 4 ++-- src/phase_mechanical.f90 | 28 ++++++++++++++-------------- 2 files changed, 16 insertions(+), 16 deletions(-) diff --git a/src/mesh/mesh_mech_FEM.f90 b/src/mesh/mesh_mech_FEM.f90 index 6bcf44c3a..d6d314a42 100644 --- a/src/mesh/mesh_mech_FEM.f90 +++ b/src/mesh/mesh_mech_FEM.f90 @@ -289,8 +289,8 @@ type(tSolutionState) function FEM_mechanical_solution( & params%timeinc = timeinc params%fieldBC = fieldBC - call SNESSolve(mechanical_snes,PETSC_NULL_VEC,solution,ierr); CHKERRQ(ierr) ! solve mechanical_snes based on solution guess (result in solution) - call SNESGetConvergedReason(mechanical_snes,reason,ierr); CHKERRQ(ierr) ! solution converged? + call SNESSolve(mechanical_snes,PETSC_NULL_VEC,solution,ierr); CHKERRQ(ierr) ! solve mechanical_snes based on solution guess (result in solution) + call SNESGetConvergedReason(mechanical_snes,reason,ierr); CHKERRQ(ierr) ! solution converged? terminallyIll = .false. if (reason < 1) then ! 0: still iterating (will not occur), negative -> convergence error diff --git a/src/phase_mechanical.f90 b/src/phase_mechanical.f90 index 3c14accb5..a14102513 100644 --- a/src/phase_mechanical.f90 +++ b/src/phase_mechanical.f90 @@ -215,8 +215,8 @@ module subroutine mechanical_init(materials,phases) allocate(phase_mechanical_F0(phases%length)) allocate(phase_mechanical_Li(phases%length)) allocate(phase_mechanical_Li0(phases%length)) - allocate(phase_mechanical_Lp0(phases%length)) allocate(phase_mechanical_Lp(phases%length)) + allocate(phase_mechanical_Lp0(phases%length)) allocate(phase_mechanical_S(phases%length)) allocate(phase_mechanical_P(phases%length)) allocate(phase_mechanical_S0(phases%length)) @@ -224,20 +224,20 @@ module subroutine mechanical_init(materials,phases) do ph = 1, phases%length Nmembers = count(material_phaseID == ph) - allocate(phase_mechanical_Fi(ph)%data(3,3,Nmembers)) allocate(phase_mechanical_Fe(ph)%data(3,3,Nmembers)) + allocate(phase_mechanical_Fi(ph)%data(3,3,Nmembers)) allocate(phase_mechanical_Fi0(ph)%data(3,3,Nmembers)) allocate(phase_mechanical_Fp(ph)%data(3,3,Nmembers)) allocate(phase_mechanical_Fp0(ph)%data(3,3,Nmembers)) - allocate(phase_mechanical_Li(ph)%data(3,3,Nmembers)) - allocate(phase_mechanical_Li0(ph)%data(3,3,Nmembers)) - allocate(phase_mechanical_Lp0(ph)%data(3,3,Nmembers)) - allocate(phase_mechanical_Lp(ph)%data(3,3,Nmembers)) + allocate(phase_mechanical_F(ph)%data(3,3,Nmembers)) + allocate(phase_mechanical_F0(ph)%data(3,3,Nmembers)) + allocate(phase_mechanical_Li(ph)%data(3,3,Nmembers),source=0.0_pReal) + allocate(phase_mechanical_Li0(ph)%data(3,3,Nmembers),source=0.0_pReal) + allocate(phase_mechanical_Lp(ph)%data(3,3,Nmembers),source=0.0_pReal) + allocate(phase_mechanical_Lp0(ph)%data(3,3,Nmembers),source=0.0_pReal) allocate(phase_mechanical_S(ph)%data(3,3,Nmembers),source=0.0_pReal) allocate(phase_mechanical_P(ph)%data(3,3,Nmembers),source=0.0_pReal) allocate(phase_mechanical_S0(ph)%data(3,3,Nmembers),source=0.0_pReal) - allocate(phase_mechanical_F(ph)%data(3,3,Nmembers)) - allocate(phase_mechanical_F0(ph)%data(3,3,Nmembers)) phase => phases%get(ph) mech => phase%get('mechanical') @@ -494,7 +494,7 @@ function integrateStress(F,subFp0,subFi0,Delta_t,co,ip,el) result(broken) enddo LpLoop call phase_LiAndItsTangents(Li_constitutive, dLi_dS, dLi_dFi, & - S, Fi_new, ph,en) + S, Fi_new, ph,en) !* update current residuum and check for convergence of loop atol_Li = max(num%rtol_crystalliteStress * max(norm2(Liguess),norm2(Li_constitutive)), & ! absolute tolerance from largest acceptable relative error @@ -1130,12 +1130,12 @@ module function phase_mechanical_dPdF(Delta_t,co,ce) result(dPdF) en = material_phaseEntry(co,ce) call phase_hooke_SandItsTangents(devNull,dSdFe,dSdFi, & - phase_mechanical_Fe(ph)%data(1:3,1:3,en), & - phase_mechanical_Fi(ph)%data(1:3,1:3,en),ph,en) + phase_mechanical_Fe(ph)%data(1:3,1:3,en), & + phase_mechanical_Fi(ph)%data(1:3,1:3,en),ph,en) call phase_LiAndItsTangents(devNull,dLidS,dLidFi, & - phase_mechanical_S(ph)%data(1:3,1:3,en), & - phase_mechanical_Fi(ph)%data(1:3,1:3,en), & - ph,en) + phase_mechanical_S(ph)%data(1:3,1:3,en), & + phase_mechanical_Fi(ph)%data(1:3,1:3,en), & + ph,en) invFp = math_inv33(phase_mechanical_Fp(ph)%data(1:3,1:3,en)) invFi = math_inv33(phase_mechanical_Fi(ph)%data(1:3,1:3,en)) From 11cdc63a2768013ec5b5e7e283431f65b6f5cab1 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Thu, 22 Jul 2021 13:50:41 +0200 Subject: [PATCH 25/63] use reasonable load case, even for compile test --- PRIVATE | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/PRIVATE b/PRIVATE index d92b030e5..f2e2d6c71 160000 --- a/PRIVATE +++ b/PRIVATE @@ -1 +1 @@ -Subproject commit d92b030e5777e718c77edc2e1e93abfa0981b024 +Subproject commit f2e2d6c71ea798bfc63230a756b7cf9748599bec From 6936ecb0910066fbfa59fcd19331a93fae4caa50 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Fri, 23 Jul 2021 21:36:26 +0200 Subject: [PATCH 26/63] sl is the subscript for slip --- .../Phase_Dislotwin_TWIP-Steel-FeMnC.yaml | 2 +- .../plastic/dislotwin_IF-steel.yaml | 2 +- src/phase_mechanical_plastic_dislotwin.f90 | 48 +++++++++---------- 3 files changed, 25 insertions(+), 27 deletions(-) diff --git a/examples/config/Phase_Dislotwin_TWIP-Steel-FeMnC.yaml b/examples/config/Phase_Dislotwin_TWIP-Steel-FeMnC.yaml index fe51585be..6c19b9f93 100644 --- a/examples/config/Phase_Dislotwin_TWIP-Steel-FeMnC.yaml +++ b/examples/config/Phase_Dislotwin_TWIP-Steel-FeMnC.yaml @@ -6,7 +6,7 @@ b_sl: [2.56e-10] rho_mob_0: [1.0e+12] rho_dip_0: [1.0] v_0: [1.0e+4] -Q_s: [3.7e-19] +Q_sl: [3.7e-19] p_sl: [1.0] q_sl: [1.0] tau_0: [1.5e+8] diff --git a/examples/config/phase/mechanical/plastic/dislotwin_IF-steel.yaml b/examples/config/phase/mechanical/plastic/dislotwin_IF-steel.yaml index 5478aa869..e48d5437f 100644 --- a/examples/config/phase/mechanical/plastic/dislotwin_IF-steel.yaml +++ b/examples/config/phase/mechanical/plastic/dislotwin_IF-steel.yaml @@ -11,7 +11,7 @@ b_sl: [2.49e-10, 2.49e-10] rho_mob_0: [2.81e12, 2.8e+12] rho_dip_0: [1.0, 1.0] # not given v_0: [1.4e+3, 1.4e+3] -Q_s: [1.57e-19, 1.57e-19] # Delta_F +Q_sl: [1.57e-19, 1.57e-19] # Delta_F tau_0: [454.e+6, 454.e+6] p_sl: [0.325, 0.325] q_sl: [1.55, 1.55] diff --git a/src/phase_mechanical_plastic_dislotwin.f90 b/src/phase_mechanical_plastic_dislotwin.f90 index 6ecf9ee7a..4aa62bcb8 100644 --- a/src/phase_mechanical_plastic_dislotwin.f90 +++ b/src/phase_mechanical_plastic_dislotwin.f90 @@ -22,7 +22,6 @@ submodule(phase:plastic) dislotwin D = 1.0_pReal, & !< grain size p_sb = 1.0_pReal, & !< p-exponent in shear band velocity q_sb = 1.0_pReal, & !< q-exponent in shear band velocity - D_a = 1.0_pReal, & !< adjustment parameter to calculate minimum dipole distance i_tw = 1.0_pReal, & !< adjustment parameter to calculate MFP for twinning L_tw = 1.0_pReal, & !< Length of twin nuclei in Burgers vectors L_tr = 1.0_pReal, & !< Length of trans nuclei in Burgers vectors @@ -42,7 +41,7 @@ submodule(phase:plastic) dislotwin b_sl, & !< absolute length of Burgers vector [m] for each slip system b_tw, & !< absolute length of Burgers vector [m] for each twin system b_tr, & !< absolute length of Burgers vector [m] for each transformation system - Q_s,& !< activation energy for glide [J] for each slip system + Q_sl,& !< activation energy for glide [J] for each slip system v_0, & !< dislocation velocity prefactor [m/s] for each slip system dot_N_0_tw, & !< twin nucleation rate [1/m³s] for each twin system dot_N_0_tr, & !< trans nucleation rate [1/m³s] for each trans system @@ -55,7 +54,8 @@ submodule(phase:plastic) dislotwin s, & !< s-exponent in trans nucleation rate tau_0, & !< strength due to elements in solid solution gamma_char, & !< characteristic shear for twins - B !< drag coefficient + B, & !< drag coefficient + d_caron !< distance of spontaneous annhihilation real(pReal), allocatable, dimension(:,:) :: & h_sl_sl, & !< components of slip-slip interaction matrix h_sl_tw, & !< components of slip-twin interaction matrix @@ -206,7 +206,7 @@ module function plastic_dislotwin_init() result(myPlasticity) rho_dip_0 = pl%get_as1dFloat('rho_dip_0', requiredSize=size(N_sl)) prm%v_0 = pl%get_as1dFloat('v_0', requiredSize=size(N_sl)) prm%b_sl = pl%get_as1dFloat('b_sl', requiredSize=size(N_sl)) - prm%Q_s = pl%get_as1dFloat('Q_s', requiredSize=size(N_sl)) + prm%Q_sl = pl%get_as1dFloat('Q_sl', requiredSize=size(N_sl)) prm%i_sl = pl%get_as1dFloat('i_sl', requiredSize=size(N_sl)) prm%p = pl%get_as1dFloat('p_sl', requiredSize=size(N_sl)) prm%q = pl%get_as1dFloat('q_sl', requiredSize=size(N_sl)) @@ -214,9 +214,8 @@ module function plastic_dislotwin_init() result(myPlasticity) prm%B = pl%get_as1dFloat('B', requiredSize=size(N_sl), & defaultVal=[(0.0_pReal, i=1,size(N_sl))]) - prm%D_a = pl%get_asFloat('D_a') - prm%D_0 = pl%get_asFloat('D_0') - prm%Q_cl = pl%get_asFloat('Q_cl') + prm%D_0 = pl%get_asFloat('D_0') + prm%Q_cl = pl%get_asFloat('Q_cl') prm%ExtendedDislocations = pl%get_asBool('extend_dislocations',defaultVal = .false.) prm%omitDipoles = pl%get_asBool('omit_dipoles',defaultVal = .false.) @@ -231,12 +230,13 @@ module function plastic_dislotwin_init() result(myPlasticity) rho_dip_0 = math_expand(rho_dip_0, N_sl) prm%v_0 = math_expand(prm%v_0, N_sl) prm%b_sl = math_expand(prm%b_sl, N_sl) - prm%Q_s = math_expand(prm%Q_s, N_sl) + prm%Q_sl = math_expand(prm%Q_sl, N_sl) prm%i_sl = math_expand(prm%i_sl, N_sl) prm%p = math_expand(prm%p, N_sl) prm%q = math_expand(prm%q, N_sl) prm%tau_0 = math_expand(prm%tau_0, N_sl) prm%B = math_expand(prm%B, N_sl) + prm%d_caron = pl%get_asFloat('D_a') * prm%b_sl ! sanity checks if ( prm%D_0 <= 0.0_pReal) extmsg = trim(extmsg)//' D_0' @@ -245,14 +245,15 @@ module function plastic_dislotwin_init() result(myPlasticity) if (any(rho_dip_0 < 0.0_pReal)) extmsg = trim(extmsg)//' rho_dip_0' if (any(prm%v_0 < 0.0_pReal)) extmsg = trim(extmsg)//' v_0' if (any(prm%b_sl <= 0.0_pReal)) extmsg = trim(extmsg)//' b_sl' - if (any(prm%Q_s <= 0.0_pReal)) extmsg = trim(extmsg)//' Q_s' + if (any(prm%Q_sl <= 0.0_pReal)) extmsg = trim(extmsg)//' Q_sl' if (any(prm%i_sl <= 0.0_pReal)) extmsg = trim(extmsg)//' i_sl' if (any(prm%B < 0.0_pReal)) extmsg = trim(extmsg)//' B' + if (any(prm%d_caron < 0.0_pReal)) extmsg = trim(extmsg)//' d_caron(D_a,b_sl)' if (any(prm%p<=0.0_pReal .or. prm%p>1.0_pReal)) extmsg = trim(extmsg)//' p_sl' if (any(prm%q< 1.0_pReal .or. prm%q>2.0_pReal)) extmsg = trim(extmsg)//' q_sl' else slipActive rho_mob_0 = emptyRealArray; rho_dip_0 = emptyRealArray - allocate(prm%b_sl,prm%Q_s,prm%v_0,prm%i_sl,prm%p,prm%q,prm%B,source=emptyRealArray) + allocate(prm%b_sl,prm%Q_sl,prm%v_0,prm%i_sl,prm%p,prm%q,prm%B,source=emptyRealArray) allocate(prm%forestProjection(0,0),prm%h_sl_sl(0,0)) endif slipActive @@ -631,7 +632,7 @@ module subroutine dislotwin_dotState(Mp,T,ph,en) integer :: i real(pReal) :: & f_unrotated, & - rho_dip_distance, & + d_hat, & v_cl, & !< climb velocity tau, & sigma_cl, & !< climb stress @@ -639,15 +640,14 @@ module subroutine dislotwin_dotState(Mp,T,ph,en) real(pReal), dimension(param(ph)%sum_N_sl) :: & dot_rho_dip_formation, & dot_rho_dip_climb, & - rho_dip_distance_min, & dot_gamma_sl real(pReal), dimension(param(ph)%sum_N_tw) :: & dot_gamma_tw real(pReal), dimension(param(ph)%sum_N_tr) :: & dot_gamma_tr - associate(prm => param(ph), stt => state(ph), & - dot => dotState(ph), dst => dependentState(ph)) + + associate(prm => param(ph), stt => state(ph), dot => dotState(ph), dst => dependentState(ph)) f_unrotated = 1.0_pReal & - sum(stt%f_tw(1:prm%sum_N_tw,en)) & @@ -656,8 +656,6 @@ module subroutine dislotwin_dotState(Mp,T,ph,en) call kinetics_sl(Mp,T,ph,en,dot_gamma_sl) dot%gamma_sl(:,en) = abs(dot_gamma_sl) - rho_dip_distance_min = prm%D_a*prm%b_sl - slipState: do i = 1, prm%sum_N_sl tau = math_tensordot(Mp,prm%P_sl(1:3,1:3,i)) @@ -665,14 +663,14 @@ module subroutine dislotwin_dotState(Mp,T,ph,en) dot_rho_dip_formation(i) = 0.0_pReal dot_rho_dip_climb(i) = 0.0_pReal else significantSlipStress - rho_dip_distance = 3.0_pReal*prm%mu*prm%b_sl(i)/(16.0_pReal*PI*abs(tau)) - rho_dip_distance = math_clip(rho_dip_distance, right = dst%Lambda_sl(i,en)) - rho_dip_distance = math_clip(rho_dip_distance, left = rho_dip_distance_min(i)) + d_hat = 3.0_pReal*prm%mu*prm%b_sl(i)/(16.0_pReal*PI*abs(tau)) + d_hat = math_clip(d_hat, right = dst%Lambda_sl(i,en)) + d_hat = math_clip(d_hat, left = prm%d_caron(i)) - dot_rho_dip_formation(i) = 2.0_pReal*(rho_dip_distance-rho_dip_distance_min(i))/prm%b_sl(i) & + dot_rho_dip_formation(i) = 2.0_pReal*(d_hat-prm%d_caron(i))/prm%b_sl(i) & * stt%rho_mob(i,en)*abs(dot_gamma_sl(i)) - if (dEq(rho_dip_distance,rho_dip_distance_min(i))) then + if (dEq(d_hat,prm%d_caron(i))) then dot_rho_dip_climb(i) = 0.0_pReal else ! Argon & Moffat, Acta Metallurgica, Vol. 29, pg 293 to 299, 1981 @@ -685,17 +683,17 @@ module subroutine dislotwin_dotState(Mp,T,ph,en) * (exp(abs(sigma_cl)*prm%b_sl(i)**3.0_pReal/(kB*T)) - 1.0_pReal) dot_rho_dip_climb(i) = 4.0_pReal*v_cl*stt%rho_dip(i,en) & - / (rho_dip_distance-rho_dip_distance_min(i)) + / (d_hat-prm%d_caron(i)) endif endif significantSlipStress enddo slipState dot%rho_mob(:,en) = abs(dot_gamma_sl)/(prm%b_sl*dst%Lambda_sl(:,en)) & - dot_rho_dip_formation & - - 2.0_pReal*rho_dip_distance_min/prm%b_sl * stt%rho_mob(:,en)*abs(dot_gamma_sl) + - 2.0_pReal*prm%d_caron/prm%b_sl * stt%rho_mob(:,en)*abs(dot_gamma_sl) dot%rho_dip(:,en) = dot_rho_dip_formation & - - 2.0_pReal*rho_dip_distance_min/prm%b_sl * stt%rho_dip(:,en)*abs(dot_gamma_sl) & + - 2.0_pReal*prm%d_caron/prm%b_sl * stt%rho_dip(:,en)*abs(dot_gamma_sl) & - dot_rho_dip_climb call kinetics_tw(Mp,T,dot_gamma_sl,ph,en,dot_gamma_tw) @@ -885,7 +883,7 @@ pure subroutine kinetics_sl(Mp,T,ph,en, & significantStress: where(tau_eff > tol_math_check) stressRatio = tau_eff/prm%tau_0 StressRatio_p = stressRatio** prm%p - BoltzmannRatio = prm%Q_s/(kB*T) + BoltzmannRatio = prm%Q_sl/(kB*T) v_wait_inverse = prm%v_0**(-1.0_pReal) * exp(BoltzmannRatio*(1.0_pReal-StressRatio_p)** prm%q) v_run_inverse = prm%B/(tau_eff*prm%b_sl) From f96b0371ac27ee38331613cbbf3a4d18b80b3405 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Fri, 23 Jul 2021 21:48:40 +0200 Subject: [PATCH 27/63] proper output name/format --- PRIVATE | 2 +- src/phase_mechanical_plastic_kinehardening.f90 | 4 ++-- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/PRIVATE b/PRIVATE index 174ecac2d..347932db3 160000 --- a/PRIVATE +++ b/PRIVATE @@ -1 +1 @@ -Subproject commit 174ecac2d3ab7596bdb60184d6bb9e1a52cb7378 +Subproject commit 347932db368ec5c0d9e9092ef106c225605d830a diff --git a/src/phase_mechanical_plastic_kinehardening.f90 b/src/phase_mechanical_plastic_kinehardening.f90 index 880b128a9..f1c9b0332 100644 --- a/src/phase_mechanical_plastic_kinehardening.f90 +++ b/src/phase_mechanical_plastic_kinehardening.f90 @@ -359,11 +359,11 @@ module subroutine plastic_kinehardening_results(ph,group) case ('xi') if(prm%sum_N_sl>0) call results_writeDataset(stt%xi,group,trim(prm%output(o)), & 'resistance against plastic slip','Pa') - case ('tau_b') !ToDo: chi + case ('chi') if(prm%sum_N_sl>0) call results_writeDataset(stt%chi,group,trim(prm%output(o)), & 'back stress','Pa') case ('sgn(gamma)') - if(prm%sum_N_sl>0) call results_writeDataset(stt%sgn_gamma,group,trim(prm%output(o)), & ! ToDo: could be int + if(prm%sum_N_sl>0) call results_writeDataset(int(stt%sgn_gamma),group,trim(prm%output(o)), & 'sense of shear','1') case ('chi_0') if(prm%sum_N_sl>0) call results_writeDataset(stt%chi_0,group,trim(prm%output(o)), & From 2d1d99542c4ab6e341c3b66bf8370e7224088f89 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Fri, 23 Jul 2021 21:51:51 +0200 Subject: [PATCH 28/63] takeover from original paper --- .../config/Phase_DisloUCLA_Tungsten.config | 42 ------------------- examples/config/Phase_Dislotungsten_W.yaml | 26 ++++++++++++ 2 files changed, 26 insertions(+), 42 deletions(-) delete mode 100644 examples/config/Phase_DisloUCLA_Tungsten.config create mode 100644 examples/config/Phase_Dislotungsten_W.yaml diff --git a/examples/config/Phase_DisloUCLA_Tungsten.config b/examples/config/Phase_DisloUCLA_Tungsten.config deleted file mode 100644 index 1eef32577..000000000 --- a/examples/config/Phase_DisloUCLA_Tungsten.config +++ /dev/null @@ -1,42 +0,0 @@ -[Tungsten] -elasticity hooke -plasticity disloucla - -(output) edge_density -(output) dipole_density -(output) shear_rate_slip -(output) accumulated_shear_slip -(output) resolved_stress_slip -(output) threshold_stress_slip - -grainsize 2.7e-5 # Average grain size [m] 2.0e-5 -SolidSolutionStrength 0.0 # Strength due to elements in solid solution - -### Dislocation glide parameters ### -#per family -Nslip 12 -slipburgers 2.72e-10 # Burgers vector of slip system [m] -rhoedge0 1.0e12 # Initial edge dislocation density [m/m**3] -rhoedgedip0 1.0 # Initial edged dipole dislocation density [m/m**3] -Qedge 2.61154e-19 # Activation energy for dislocation glide [J], 1.63 eV -v0 1 # Initial glide velocity [m/s] -p_slip 0.86 # p-exponent in glide velocity -q_slip 1.69 # q-exponent in glide velocity -tau_peierls 2.03e9 # peierls stress [Pa] - -#mobility law -kink_height 2.567e-10 # kink height sqrt(6)/3*lattice_parameter [m] -omega 9.1e11 # attemp frequency (from kMC paper) [s^(-1)] -kink_width 29.95e-10 # kink pair width ~ 11 b (kMC paper) [m] -dislolength 78e-10 # dislocation length (ideally lambda) [m] initial value 11b -friction_coeff 8.3e-5 # friction coeff. B [Pa*s] - -#hardening -dipoleformationfactor 0 # to have hardening due to dipole formation off -CLambdaSlip 10.0 # Adj. parameter controlling dislocation mean free path -D0 4.0e-5 # Vacancy diffusion prefactor [m**2/s] -Qsd 4.5e-19 # Activation energy for climb [J] -Catomicvolume 1.0 # Adj. parameter controlling the atomic volume [in b] -Cedgedipmindistance 1.0 # Adj. parameter controlling the minimum dipole distance [in b] -interaction_slipslip 0.009 0.72 0.009 0.05 0.05 0.06 0.09 -nonschmid_coefficients 0.938 0.71 4.43 0.0 0.0 0.0 diff --git a/examples/config/Phase_Dislotungsten_W.yaml b/examples/config/Phase_Dislotungsten_W.yaml new file mode 100644 index 000000000..38cf2d4e4 --- /dev/null +++ b/examples/config/Phase_Dislotungsten_W.yaml @@ -0,0 +1,26 @@ +type: dislotungsten + +N_sl: [12] + +rho_mob_0: [1.0e+9] +rho_dip_0: [1.0] + +nu_a: [9.1e+11] +b_sl: [2.72e-10] +Delta_H_kp,0: [2.61154e-19] # 1.63 eV, Delta_H0 + +tau_Peierls: [2.03e+9] +p_sl: [0.86] +q_sl: [1.69] +h: [2.566e-10] +w: [2.992e-09] +B: [8.3e-5] +D_a: 1.0 # d_edge + +# climb (disabled) +D_0: 0.0 +Q_cl: 0.0 +V_cl: [0.0] + +h_sl-sl: [0.009, 0.72, 0.009, 0.05, 0.05, 0.06, 0.09] +a_nonSchmid: [0.938, 0.71, 4.43] From c46813657f38e0a898bba78eb2d78fc71dfc1d60 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Fri, 23 Jul 2021 21:56:34 +0200 Subject: [PATCH 29/63] better to understand --- ...phase_mechanical_plastic_dislotungsten.f90 | 33 +++++++++---------- 1 file changed, 15 insertions(+), 18 deletions(-) diff --git a/src/phase_mechanical_plastic_dislotungsten.f90 b/src/phase_mechanical_plastic_dislotungsten.f90 index 6ca8e440e..0fc0076fc 100644 --- a/src/phase_mechanical_plastic_dislotungsten.f90 +++ b/src/phase_mechanical_plastic_dislotungsten.f90 @@ -56,7 +56,7 @@ submodule(phase:plastic) dislotungsten type :: tDisloTungstendependentState real(pReal), dimension(:,:), allocatable :: & Lambda_sl, & - threshold_stress + tau_pass end type tDisloTungstendependentState !-------------------------------------------------------------------------------------------------- @@ -246,8 +246,8 @@ module function plastic_dislotungsten_init() result(myPlasticity) plasticState(ph)%atol(startIndex:endIndex) = pl%get_asFloat('atol_gamma',defaultVal=1.0e-6_pReal) if(any(plasticState(ph)%atol(startIndex:endIndex) < 0.0_pReal)) extmsg = trim(extmsg)//' atol_gamma' - allocate(dst%Lambda_sl(prm%sum_N_sl,Nmembers), source=0.0_pReal) - allocate(dst%threshold_stress(prm%sum_N_sl,Nmembers), source=0.0_pReal) + allocate(dst%Lambda_sl(prm%sum_N_sl,Nmembers), source=0.0_pReal) + allocate(dst%tau_pass(prm%sum_N_sl,Nmembers), source=0.0_pReal) end associate @@ -380,8 +380,8 @@ module subroutine dislotungsten_dependentState(ph,en) associate(prm => param(ph), stt => state(ph),dst => dependentState(ph)) dislocationSpacing = sqrt(matmul(prm%forestProjection,stt%rho_mob(:,en)+stt%rho_dip(:,en))) - dst%threshold_stress(:,en) = prm%mu*prm%b_sl & - * sqrt(matmul(prm%h_sl_sl,stt%rho_mob(:,en)+stt%rho_dip(:,en))) + dst%tau_pass(:,en) = prm%mu*prm%b_sl & + * sqrt(matmul(prm%h_sl_sl,stt%rho_mob(:,en)+stt%rho_dip(:,en))) dst%Lambda_sl(:,en) = prm%D/(1.0_pReal+prm%D*dislocationSpacing/prm%i_sl) @@ -416,7 +416,7 @@ module subroutine plastic_dislotungsten_results(ph,group) if(prm%sum_N_sl>0) call results_writeDataset(dst%Lambda_sl,group,trim(prm%output(o)), & 'mean free path for slip','m') case('tau_pass') - if(prm%sum_N_sl>0) call results_writeDataset(dst%threshold_stress,group,trim(prm%output(o)), & + if(prm%sum_N_sl>0) call results_writeDataset(dst%tau_pass,group,trim(prm%output(o)), & 'threshold stress for slip','Pa') end select enddo outputsLoop @@ -456,8 +456,7 @@ pure subroutine kinetics(Mp,T,ph,en, & StressRatio_p,StressRatio_pminus1, & dvel, vel, & tau_pos,tau_neg, & - t_n, t_k, dtk,dtn, & - needsGoodName ! ToDo: @Karo: any idea? + t_n, t_k, dtk,dtn integer :: j associate(prm => param(ph), stt => state(ph), dst => dependentState(ph)) @@ -475,13 +474,12 @@ pure subroutine kinetics(Mp,T,ph,en, & dot_gamma_0 => stt%rho_mob(:,en)*prm%b_sl*prm%v_0, & effectiveLength => dst%Lambda_sl(:,en) - prm%w) - significantPositiveTau: where(abs(tau_pos)-dst%threshold_stress(:,en) > tol_math_check) - StressRatio = (abs(tau_pos)-dst%threshold_stress(:,en))/prm%tau_Peierls + significantPositiveTau: where(abs(tau_pos)-dst%tau_pass(:,en) > tol_math_check) + StressRatio = (abs(tau_pos)-dst%tau_pass(:,en))/prm%tau_Peierls StressRatio_p = StressRatio** prm%p StressRatio_pminus1 = StressRatio**(prm%p-1.0_pReal) - needsGoodName = exp(-BoltzmannRatio*(1-StressRatio_p) ** prm%q) - t_n = prm%b_sl/(needsGoodName*prm%omega*effectiveLength) + t_n = prm%b_sl/(exp(-BoltzmannRatio*(1-StressRatio_p) ** prm%q)*prm%omega*effectiveLength) t_k = effectiveLength * prm%B /(2.0_pReal*prm%b_sl*tau_pos) vel = prm%h/(t_n + t_k) @@ -492,7 +490,7 @@ pure subroutine kinetics(Mp,T,ph,en, & end where significantPositiveTau if (present(ddot_gamma_dtau_pos)) then - significantPositiveTau2: where(abs(tau_pos)-dst%threshold_stress(:,en) > tol_math_check) + significantPositiveTau2: where(abs(tau_pos)-dst%tau_pass(:,en) > tol_math_check) dtn = -1.0_pReal * t_n * BoltzmannRatio * prm%p * prm%q * (1.0_pReal-StressRatio_p)**(prm%q - 1.0_pReal) & * (StressRatio)**(prm%p - 1.0_pReal) / prm%tau_Peierls dtk = -1.0_pReal * t_k / tau_pos @@ -505,13 +503,12 @@ pure subroutine kinetics(Mp,T,ph,en, & end where significantPositiveTau2 endif - significantNegativeTau: where(abs(tau_neg)-dst%threshold_stress(:,en) > tol_math_check) - StressRatio = (abs(tau_neg)-dst%threshold_stress(:,en))/prm%tau_Peierls + significantNegativeTau: where(abs(tau_neg)-dst%tau_pass(:,en) > tol_math_check) + StressRatio = (abs(tau_neg)-dst%tau_pass(:,en))/prm%tau_Peierls StressRatio_p = StressRatio** prm%p StressRatio_pminus1 = StressRatio**(prm%p-1.0_pReal) - needsGoodName = exp(-BoltzmannRatio*(1-StressRatio_p) ** prm%q) - t_n = prm%b_sl/(needsGoodName*prm%omega*effectiveLength) + t_n = prm%b_sl/(exp(-BoltzmannRatio*(1-StressRatio_p) ** prm%q)*prm%omega*effectiveLength) t_k = effectiveLength * prm%B /(2.0_pReal*prm%b_sl*tau_pos) vel = prm%h/(t_n + t_k) @@ -522,7 +519,7 @@ pure subroutine kinetics(Mp,T,ph,en, & end where significantNegativeTau if (present(ddot_gamma_dtau_neg)) then - significantNegativeTau2: where(abs(tau_neg)-dst%threshold_stress(:,en) > tol_math_check) + significantNegativeTau2: where(abs(tau_neg)-dst%tau_pass(:,en) > tol_math_check) dtn = -1.0_pReal * t_n * BoltzmannRatio * prm%p * prm%q * (1.0_pReal-StressRatio_p)**(prm%q - 1.0_pReal) & * (StressRatio)**(prm%p - 1.0_pReal) / prm%tau_Peierls dtk = -1.0_pReal * t_k / tau_neg From 2a78174547b473f1d96ff5916ad5e805c2f53c88 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Fri, 23 Jul 2021 22:04:39 +0200 Subject: [PATCH 30/63] standard names in source code --- examples/config/Phase_Dislotungsten_W.yaml | 2 +- ...phase_mechanical_plastic_dislotungsten.f90 | 36 +++++++++---------- 2 files changed, 17 insertions(+), 21 deletions(-) diff --git a/examples/config/Phase_Dislotungsten_W.yaml b/examples/config/Phase_Dislotungsten_W.yaml index 38cf2d4e4..bf8796cfa 100644 --- a/examples/config/Phase_Dislotungsten_W.yaml +++ b/examples/config/Phase_Dislotungsten_W.yaml @@ -16,7 +16,7 @@ h: [2.566e-10] w: [2.992e-09] B: [8.3e-5] D_a: 1.0 # d_edge - + # climb (disabled) D_0: 0.0 Q_cl: 0.0 diff --git a/src/phase_mechanical_plastic_dislotungsten.f90 b/src/phase_mechanical_plastic_dislotungsten.f90 index 0fc0076fc..e7cf32080 100644 --- a/src/phase_mechanical_plastic_dislotungsten.f90 +++ b/src/phase_mechanical_plastic_dislotungsten.f90 @@ -18,7 +18,7 @@ submodule(phase:plastic) dislotungsten Q_cl = 1.0_pReal !< activation energy for dislocation climb real(pReal), allocatable, dimension(:) :: & b_sl, & !< magnitude of Burgers vector [m] - D_a, & + d_caron, & !< distance of spontaneous annhihilation i_sl, & !< Adj. parameter for distance between 2 forest dislocations f_at, & !< factor to calculate atomic volume tau_Peierls, & !< Peierls stress @@ -172,7 +172,6 @@ module function plastic_dislotungsten_init() result(myPlasticity) prm%D_0 = pl%get_asFloat('D_0') prm%Q_cl = pl%get_asFloat('Q_cl') prm%f_at = pl%get_asFloat('f_at') * prm%b_sl**3.0_pReal - prm%D_a = pl%get_asFloat('D_a') * prm%b_sl prm%dipoleformation = .not. pl%get_asBool('no_dipole_formation', defaultVal = .false.) @@ -191,7 +190,7 @@ module function plastic_dislotungsten_init() result(myPlasticity) prm%B = math_expand(prm%B, N_sl) prm%i_sl = math_expand(prm%i_sl, N_sl) prm%f_at = math_expand(prm%f_at, N_sl) - prm%D_a = math_expand(prm%D_a, N_sl) + prm%d_caron = pl%get_asFloat('D_a') * prm%b_sl ! sanity checks if ( prm%D_0 <= 0.0_pReal) extmsg = trim(extmsg)//' D_0' @@ -202,12 +201,13 @@ module function plastic_dislotungsten_init() result(myPlasticity) if (any(prm%b_sl <= 0.0_pReal)) extmsg = trim(extmsg)//' b_sl' if (any(prm%Q_s <= 0.0_pReal)) extmsg = trim(extmsg)//' Q_s' if (any(prm%tau_Peierls < 0.0_pReal)) extmsg = trim(extmsg)//' tau_Peierls' - if (any(prm%D_a <= 0.0_pReal)) extmsg = trim(extmsg)//' D_a or b_sl' + if (any(prm%B < 0.0_pReal)) extmsg = trim(extmsg)//' B' + if (any(prm%d_caron < 0.0_pReal)) extmsg = trim(extmsg)//' d_caron(D_a,b_sl)' if (any(prm%f_at <= 0.0_pReal)) extmsg = trim(extmsg)//' f_at or b_sl' else slipActive rho_mob_0= emptyRealArray; rho_dip_0 = emptyRealArray - allocate(prm%b_sl,prm%D_a,prm%i_sl,prm%f_at,prm%tau_Peierls, & + allocate(prm%b_sl,prm%d_caron,prm%i_sl,prm%f_at,prm%tau_Peierls, & prm%Q_s,prm%v_0,prm%p,prm%q,prm%B,prm%h,prm%w,prm%omega, & source = emptyRealArray) allocate(prm%forestProjection(0,0)) @@ -316,8 +316,6 @@ module subroutine dislotungsten_dotState(Mp,T,ph,en) ph, & en - real(pReal) :: & - VacancyDiffusion real(pReal), dimension(param(ph)%sum_N_sl) :: & dot_gamma_pos, dot_gamma_neg,& tau_pos,& @@ -325,38 +323,36 @@ module subroutine dislotungsten_dotState(Mp,T,ph,en) v_cl, & dot_rho_dip_formation, & dot_rho_dip_climb, & - dip_distance + d_hat - associate(prm => param(ph), stt => state(ph),& - dot => dotState(ph), dst => dependentState(ph)) + associate(prm => param(ph), stt => state(ph), dot => dotState(ph), dst => dependentState(ph)) call kinetics(Mp,T,ph,en,& dot_gamma_pos,dot_gamma_neg, & tau_pos_out = tau_pos,tau_neg_out = tau_neg) dot%gamma_sl(:,en) = (dot_gamma_pos+dot_gamma_neg) ! ToDo: needs to be abs - VacancyDiffusion = prm%D_0*exp(-prm%Q_cl/(kB*T)) where(dEq0(tau_pos)) ! ToDo: use avg of pos and neg dot_rho_dip_formation = 0.0_pReal dot_rho_dip_climb = 0.0_pReal else where - dip_distance = math_clip(3.0_pReal*prm%mu*prm%b_sl/(16.0_pReal*PI*abs(tau_pos)), & - prm%D_a, & ! lower limit - dst%Lambda_sl(:,en)) ! upper limit - dot_rho_dip_formation = merge(2.0_pReal*dip_distance* stt%rho_mob(:,en)*abs(dot%gamma_sl(:,en))/prm%b_sl, & ! ToDo: ignore region of spontaneous annihilation + d_hat = math_clip(3.0_pReal*prm%mu*prm%b_sl/(16.0_pReal*PI*abs(tau_pos)), & + prm%d_caron, & ! lower limit + dst%Lambda_sl(:,en)) ! upper limit + dot_rho_dip_formation = merge(2.0_pReal*d_hat* stt%rho_mob(:,en)*abs(dot%gamma_sl(:,en))/prm%b_sl, & ! ToDo: ignore region of spontaneous annihilation 0.0_pReal, & prm%dipoleformation) - v_cl = (3.0_pReal*prm%mu*VacancyDiffusion*prm%f_at/(2.0_pReal*pi*kB*T)) & - * (1.0_pReal/(dip_distance+prm%D_a)) - dot_rho_dip_climb = (4.0_pReal*v_cl*stt%rho_dip(:,en))/(dip_distance-prm%D_a) ! ToDo: Discuss with Franz: Stress dependency? + v_cl = (3.0_pReal*prm%mu*prm%D_0*exp(-prm%Q_cl/(kB*T))*prm%f_at/(2.0_pReal*PI*kB*T)) & + * (1.0_pReal/(d_hat+prm%d_caron)) + dot_rho_dip_climb = (4.0_pReal*v_cl*stt%rho_dip(:,en))/(d_hat-prm%d_caron) ! ToDo: Discuss with Franz: Stress dependency? end where dot%rho_mob(:,en) = abs(dot%gamma_sl(:,en))/(prm%b_sl*dst%Lambda_sl(:,en)) & ! multiplication - dot_rho_dip_formation & - - (2.0_pReal*prm%D_a)/prm%b_sl*stt%rho_mob(:,en)*abs(dot%gamma_sl(:,en)) ! Spontaneous annihilation of 2 single edge dislocations + - (2.0_pReal*prm%d_caron)/prm%b_sl*stt%rho_mob(:,en)*abs(dot%gamma_sl(:,en)) ! Spontaneous annihilation of 2 edges dot%rho_dip(:,en) = dot_rho_dip_formation & - - (2.0_pReal*prm%D_a)/prm%b_sl*stt%rho_dip(:,en)*abs(dot%gamma_sl(:,en)) & ! Spontaneous annihilation of a single edge dislocation with a dipole constituent + - (2.0_pReal*prm%d_caron)/prm%b_sl*stt%rho_dip(:,en)*abs(dot%gamma_sl(:,en)) & ! Spontaneous annihilation of an edge with a dipole - dot_rho_dip_climb end associate From 527fd306e27b4f7731474ce9d945644f9a3ca2f2 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Fri, 23 Jul 2021 22:25:06 +0200 Subject: [PATCH 31/63] change of behavior - gamma (state) increases monotoneously - region of spontaneous annihilation is ignored in dipole formation --- src/phase_mechanical_plastic_dislotungsten.f90 | 16 ++++++++-------- 1 file changed, 8 insertions(+), 8 deletions(-) diff --git a/src/phase_mechanical_plastic_dislotungsten.f90 b/src/phase_mechanical_plastic_dislotungsten.f90 index e7cf32080..4f86b23b1 100644 --- a/src/phase_mechanical_plastic_dislotungsten.f90 +++ b/src/phase_mechanical_plastic_dislotungsten.f90 @@ -331,16 +331,16 @@ module subroutine dislotungsten_dotState(Mp,T,ph,en) dot_gamma_pos,dot_gamma_neg, & tau_pos_out = tau_pos,tau_neg_out = tau_neg) - dot%gamma_sl(:,en) = (dot_gamma_pos+dot_gamma_neg) ! ToDo: needs to be abs + dot%gamma_sl(:,en) = abs(dot_gamma_pos+dot_gamma_neg) - where(dEq0(tau_pos)) ! ToDo: use avg of pos and neg + where(dEq0(tau_pos)) ! ToDo: use avg of +/- dot_rho_dip_formation = 0.0_pReal dot_rho_dip_climb = 0.0_pReal else where - d_hat = math_clip(3.0_pReal*prm%mu*prm%b_sl/(16.0_pReal*PI*abs(tau_pos)), & + d_hat = math_clip(3.0_pReal*prm%mu*prm%b_sl/(16.0_pReal*PI*abs(tau_pos)), & ! ToDo: use avg of +/- prm%d_caron, & ! lower limit dst%Lambda_sl(:,en)) ! upper limit - dot_rho_dip_formation = merge(2.0_pReal*d_hat* stt%rho_mob(:,en)*abs(dot%gamma_sl(:,en))/prm%b_sl, & ! ToDo: ignore region of spontaneous annihilation + dot_rho_dip_formation = merge(2.0_pReal*(d_hat-prm%d_caron)*stt%rho_mob(:,en)*dot%gamma_sl(:,en)/prm%b_sl, & 0.0_pReal, & prm%dipoleformation) v_cl = (3.0_pReal*prm%mu*prm%D_0*exp(-prm%Q_cl/(kB*T))*prm%f_at/(2.0_pReal*PI*kB*T)) & @@ -348,11 +348,11 @@ module subroutine dislotungsten_dotState(Mp,T,ph,en) dot_rho_dip_climb = (4.0_pReal*v_cl*stt%rho_dip(:,en))/(d_hat-prm%d_caron) ! ToDo: Discuss with Franz: Stress dependency? end where - dot%rho_mob(:,en) = abs(dot%gamma_sl(:,en))/(prm%b_sl*dst%Lambda_sl(:,en)) & ! multiplication + dot%rho_mob(:,en) = dot%gamma_sl(:,en)/(prm%b_sl*dst%Lambda_sl(:,en)) & ! multiplication - dot_rho_dip_formation & - - (2.0_pReal*prm%d_caron)/prm%b_sl*stt%rho_mob(:,en)*abs(dot%gamma_sl(:,en)) ! Spontaneous annihilation of 2 edges + - (2.0_pReal*prm%d_caron)/prm%b_sl*stt%rho_mob(:,en)*dot%gamma_sl(:,en) ! Spontaneous annihilation of 2 edges dot%rho_dip(:,en) = dot_rho_dip_formation & - - (2.0_pReal*prm%d_caron)/prm%b_sl*stt%rho_dip(:,en)*abs(dot%gamma_sl(:,en)) & ! Spontaneous annihilation of an edge with a dipole + - (2.0_pReal*prm%d_caron)/prm%b_sl*stt%rho_dip(:,en)*dot%gamma_sl(:,en) & ! Spontaneous annihilation of an edge with a dipole - dot_rho_dip_climb end associate @@ -373,7 +373,7 @@ module subroutine dislotungsten_dependentState(ph,en) dislocationSpacing - associate(prm => param(ph), stt => state(ph),dst => dependentState(ph)) + associate(prm => param(ph), stt => state(ph), dst => dependentState(ph)) dislocationSpacing = sqrt(matmul(prm%forestProjection,stt%rho_mob(:,en)+stt%rho_dip(:,en))) dst%tau_pass(:,en) = prm%mu*prm%b_sl & From d4ffc778c2d491f7930c488197bda23f56c2d071 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sat, 24 Jul 2021 07:09:35 +0200 Subject: [PATCH 32/63] easier to read removed comment regarding use of dot_state in kinetics_t(w/r). Data stored in dotState is not reliable, FPI integrator for writes to it and Runge-Kutta calls the dot state function at different time steps --- src/phase_mechanical_plastic_dislotwin.f90 | 243 +++++++++++---------- 1 file changed, 123 insertions(+), 120 deletions(-) diff --git a/src/phase_mechanical_plastic_dislotwin.f90 b/src/phase_mechanical_plastic_dislotwin.f90 index 4aa62bcb8..1d6fadd13 100644 --- a/src/phase_mechanical_plastic_dislotwin.f90 +++ b/src/phase_mechanical_plastic_dislotwin.f90 @@ -517,7 +517,7 @@ module subroutine dislotwin_LpAndItsTangent(Lp,dLp_dMp,Mp,T,ph,en) integer :: i,k,l,m,n real(pReal) :: & f_unrotated,StressRatio_p,& - BoltzmannRatio, & + E_kB_T, & ddot_gamma_dtau, & tau real(pReal), dimension(param(ph)%sum_N_sl) :: & @@ -587,7 +587,7 @@ module subroutine dislotwin_LpAndItsTangent(Lp,dLp_dMp,Mp,T,ph,en) shearBandingContribution: if(dNeq0(prm%v_sb)) then - BoltzmannRatio = prm%E_sb/(kB*T) + E_kB_T = prm%E_sb/(kB*T) call math_eigh33(eigValues,eigVectors,Mp) ! is Mp symmetric by design? do i = 1,6 @@ -597,8 +597,8 @@ module subroutine dislotwin_LpAndItsTangent(Lp,dLp_dMp,Mp,T,ph,en) significantShearBandStress: if (abs(tau) > tol_math_check) then StressRatio_p = (abs(tau)/prm%xi_sb)**prm%p_sb - dot_gamma_sb = sign(prm%v_sb*exp(-BoltzmannRatio*(1-StressRatio_p)**prm%q_sb), tau) - ddot_gamma_dtau = abs(dot_gamma_sb)*BoltzmannRatio* prm%p_sb*prm%q_sb/ prm%xi_sb & + dot_gamma_sb = sign(prm%v_sb*exp(-E_kB_T*(1-StressRatio_p)**prm%q_sb), tau) + ddot_gamma_dtau = abs(dot_gamma_sb)*E_kB_T*prm%p_sb*prm%q_sb/prm%xi_sb & * (abs(tau)/prm%xi_sb)**(prm%p_sb-1.0_pReal) & * (1.0_pReal-StressRatio_p)**(prm%q_sb-1.0_pReal) @@ -649,58 +649,58 @@ module subroutine dislotwin_dotState(Mp,T,ph,en) associate(prm => param(ph), stt => state(ph), dot => dotState(ph), dst => dependentState(ph)) - f_unrotated = 1.0_pReal & - - sum(stt%f_tw(1:prm%sum_N_tw,en)) & - - sum(stt%f_tr(1:prm%sum_N_tr,en)) + f_unrotated = 1.0_pReal & + - sum(stt%f_tw(1:prm%sum_N_tw,en)) & + - sum(stt%f_tr(1:prm%sum_N_tr,en)) - call kinetics_sl(Mp,T,ph,en,dot_gamma_sl) - dot%gamma_sl(:,en) = abs(dot_gamma_sl) + call kinetics_sl(Mp,T,ph,en,dot_gamma_sl) + dot%gamma_sl(:,en) = abs(dot_gamma_sl) - slipState: do i = 1, prm%sum_N_sl - tau = math_tensordot(Mp,prm%P_sl(1:3,1:3,i)) + slipState: do i = 1, prm%sum_N_sl + tau = math_tensordot(Mp,prm%P_sl(1:3,1:3,i)) - significantSlipStress: if (dEq0(tau) .or. prm%omitDipoles) then - dot_rho_dip_formation(i) = 0.0_pReal - dot_rho_dip_climb(i) = 0.0_pReal - else significantSlipStress - d_hat = 3.0_pReal*prm%mu*prm%b_sl(i)/(16.0_pReal*PI*abs(tau)) - d_hat = math_clip(d_hat, right = dst%Lambda_sl(i,en)) - d_hat = math_clip(d_hat, left = prm%d_caron(i)) - - dot_rho_dip_formation(i) = 2.0_pReal*(d_hat-prm%d_caron(i))/prm%b_sl(i) & - * stt%rho_mob(i,en)*abs(dot_gamma_sl(i)) - - if (dEq(d_hat,prm%d_caron(i))) then + significantSlipStress: if (dEq0(tau) .or. prm%omitDipoles) then + dot_rho_dip_formation(i) = 0.0_pReal dot_rho_dip_climb(i) = 0.0_pReal - else - ! Argon & Moffat, Acta Metallurgica, Vol. 29, pg 293 to 299, 1981 - sigma_cl = dot_product(prm%n0_sl(1:3,i),matmul(Mp,prm%n0_sl(1:3,i))) - b_d = merge(24.0_pReal*PI*(1.0_pReal - prm%nu)/(2.0_pReal + prm%nu) & - * (prm%Gamma_sf(1) + prm%Gamma_sf(2) * T) / (prm%mu*prm%b_sl(i)), & - 1.0_pReal, & - prm%ExtendedDislocations) - v_cl = 2.0_pReal*prm%omega*b_d**2.0_pReal*exp(-prm%Q_cl/(kB*T)) & - * (exp(abs(sigma_cl)*prm%b_sl(i)**3.0_pReal/(kB*T)) - 1.0_pReal) + else significantSlipStress + d_hat = 3.0_pReal*prm%mu*prm%b_sl(i)/(16.0_pReal*PI*abs(tau)) + d_hat = math_clip(d_hat, right = dst%Lambda_sl(i,en)) + d_hat = math_clip(d_hat, left = prm%d_caron(i)) - dot_rho_dip_climb(i) = 4.0_pReal*v_cl*stt%rho_dip(i,en) & - / (d_hat-prm%d_caron(i)) - endif - endif significantSlipStress - enddo slipState + dot_rho_dip_formation(i) = 2.0_pReal*(d_hat-prm%d_caron(i))/prm%b_sl(i) & + * stt%rho_mob(i,en)*abs(dot_gamma_sl(i)) - dot%rho_mob(:,en) = abs(dot_gamma_sl)/(prm%b_sl*dst%Lambda_sl(:,en)) & - - dot_rho_dip_formation & - - 2.0_pReal*prm%d_caron/prm%b_sl * stt%rho_mob(:,en)*abs(dot_gamma_sl) + if (dEq(d_hat,prm%d_caron(i))) then + dot_rho_dip_climb(i) = 0.0_pReal + else + ! Argon & Moffat, Acta Metallurgica, Vol. 29, pg 293 to 299, 1981 + sigma_cl = dot_product(prm%n0_sl(1:3,i),matmul(Mp,prm%n0_sl(1:3,i))) + b_d = merge(24.0_pReal*PI*(1.0_pReal - prm%nu)/(2.0_pReal + prm%nu) & + * (prm%Gamma_sf(1) + prm%Gamma_sf(2) * T) / (prm%mu*prm%b_sl(i)), & + 1.0_pReal, & + prm%ExtendedDislocations) + v_cl = 2.0_pReal*prm%omega*b_d**2.0_pReal*exp(-prm%Q_cl/(kB*T)) & + * (exp(abs(sigma_cl)*prm%b_sl(i)**3.0_pReal/(kB*T)) - 1.0_pReal) - dot%rho_dip(:,en) = dot_rho_dip_formation & - - 2.0_pReal*prm%d_caron/prm%b_sl * stt%rho_dip(:,en)*abs(dot_gamma_sl) & - - dot_rho_dip_climb + dot_rho_dip_climb(i) = 4.0_pReal*v_cl*stt%rho_dip(i,en) & + / (d_hat-prm%d_caron(i)) + endif + endif significantSlipStress + enddo slipState - call kinetics_tw(Mp,T,dot_gamma_sl,ph,en,dot_gamma_tw) - dot%f_tw(:,en) = f_unrotated*dot_gamma_tw/prm%gamma_char + dot%rho_mob(:,en) = abs(dot_gamma_sl)/(prm%b_sl*dst%Lambda_sl(:,en)) & + - dot_rho_dip_formation & + - 2.0_pReal*prm%d_caron/prm%b_sl * stt%rho_mob(:,en)*abs(dot_gamma_sl) - call kinetics_tr(Mp,T,dot_gamma_sl,ph,en,dot_gamma_tr) - dot%f_tr(:,en) = f_unrotated*dot_gamma_tr + dot%rho_dip(:,en) = dot_rho_dip_formation & + - 2.0_pReal*prm%d_caron/prm%b_sl * stt%rho_dip(:,en)*abs(dot_gamma_sl) & + - dot_rho_dip_climb + + call kinetics_tw(Mp,T,dot_gamma_sl,ph,en,dot_gamma_tw) + dot%f_tw(:,en) = f_unrotated*dot_gamma_tw/prm%gamma_char + + call kinetics_tr(Mp,T,dot_gamma_sl,ph,en,dot_gamma_tr) + dot%f_tr(:,en) = f_unrotated*dot_gamma_tr end associate @@ -865,7 +865,7 @@ pure subroutine kinetics_sl(Mp,T,ph,en, & tau, & stressRatio, & StressRatio_p, & - BoltzmannRatio, & + Q_kB_T, & v_wait_inverse, & !< inverse of the effective velocity of a dislocation waiting at obstacles (unsigned) v_run_inverse, & !< inverse of the velocity of a free moving dislocation (unsigned) dV_wait_inverse_dTau, & @@ -874,33 +874,34 @@ pure subroutine kinetics_sl(Mp,T,ph,en, & tau_eff !< effective resolved stress integer :: i + associate(prm => param(ph), stt => state(ph), dst => dependentState(ph)) - tau = [(math_tensordot(Mp,prm%P_sl(1:3,1:3,i)),i = 1, prm%sum_N_sl)] + tau = [(math_tensordot(Mp,prm%P_sl(1:3,1:3,i)),i = 1, prm%sum_N_sl)] - tau_eff = abs(tau)-dst%tau_pass(:,en) + tau_eff = abs(tau)-dst%tau_pass(:,en) - significantStress: where(tau_eff > tol_math_check) - stressRatio = tau_eff/prm%tau_0 - StressRatio_p = stressRatio** prm%p - BoltzmannRatio = prm%Q_sl/(kB*T) - v_wait_inverse = prm%v_0**(-1.0_pReal) * exp(BoltzmannRatio*(1.0_pReal-StressRatio_p)** prm%q) - v_run_inverse = prm%B/(tau_eff*prm%b_sl) + significantStress: where(tau_eff > tol_math_check) + stressRatio = tau_eff/prm%tau_0 + StressRatio_p = stressRatio** prm%p + Q_kB_T = prm%Q_sl/(kB*T) + v_wait_inverse = prm%v_0**(-1.0_pReal) * exp(Q_kB_T*(1.0_pReal-StressRatio_p)** prm%q) + v_run_inverse = prm%B/(tau_eff*prm%b_sl) - dot_gamma_sl = sign(stt%rho_mob(:,en)*prm%b_sl/(v_wait_inverse+v_run_inverse),tau) + dot_gamma_sl = sign(stt%rho_mob(:,en)*prm%b_sl/(v_wait_inverse+v_run_inverse),tau) - dV_wait_inverse_dTau = -1.0_pReal * v_wait_inverse * prm%p * prm%q * BoltzmannRatio & - * (stressRatio**(prm%p-1.0_pReal)) & - * (1.0_pReal-StressRatio_p)**(prm%q-1.0_pReal) & - / prm%tau_0 - dV_run_inverse_dTau = -1.0_pReal * v_run_inverse/tau_eff - dV_dTau = -1.0_pReal * (dV_wait_inverse_dTau+dV_run_inverse_dTau) & - / (v_wait_inverse+v_run_inverse)**2.0_pReal - ddot_gamma_dtau = dV_dTau*stt%rho_mob(:,en)*prm%b_sl - else where significantStress - dot_gamma_sl = 0.0_pReal - ddot_gamma_dtau = 0.0_pReal - end where significantStress + dV_wait_inverse_dTau = -1.0_pReal * v_wait_inverse * prm%p * prm%q * Q_kB_T & + * (stressRatio**(prm%p-1.0_pReal)) & + * (1.0_pReal-StressRatio_p)**(prm%q-1.0_pReal) & + / prm%tau_0 + dV_run_inverse_dTau = -1.0_pReal * v_run_inverse/tau_eff + dV_dTau = -1.0_pReal * (dV_wait_inverse_dTau+dV_run_inverse_dTau) & + / (v_wait_inverse+v_run_inverse)**2.0_pReal + ddot_gamma_dtau = dV_dTau*stt%rho_mob(:,en)*prm%b_sl + else where significantStress + dot_gamma_sl = 0.0_pReal + ddot_gamma_dtau = 0.0_pReal + end where significantStress end associate @@ -943,34 +944,35 @@ pure subroutine kinetics_tw(Mp,T,dot_gamma_sl,ph,en,& integer :: i,s1,s2 + associate(prm => param(ph), stt => state(ph), dst => dependentState(ph)) - do i = 1, prm%sum_N_tw - tau(i) = math_tensordot(Mp,prm%P_tw(1:3,1:3,i)) - isFCC: if (prm%fccTwinTransNucleation) then - s1=prm%fcc_twinNucleationSlipPair(1,i) - s2=prm%fcc_twinNucleationSlipPair(2,i) - if (tau(i) < dst%tau_r_tw(i,en)) then ! ToDo: correct? - Ndot0=(abs(dot_gamma_sl(s1))*(stt%rho_mob(s2,en)+stt%rho_dip(s2,en))+& - abs(dot_gamma_sl(s2))*(stt%rho_mob(s1,en)+stt%rho_dip(s1,en)))/& ! ToDo: MD: it would be more consistent to use shearrates from state - (prm%L_tw*prm%b_sl(i))*& - (1.0_pReal-exp(-prm%V_cs/(kB*T)*(dst%tau_r_tw(i,en)-tau(i)))) ! P_ncs - else - Ndot0=0.0_pReal - end if - else isFCC - Ndot0=prm%dot_N_0_tw(i) - endif isFCC - enddo + do i = 1, prm%sum_N_tw + tau(i) = math_tensordot(Mp,prm%P_tw(1:3,1:3,i)) + isFCC: if (prm%fccTwinTransNucleation) then + s1=prm%fcc_twinNucleationSlipPair(1,i) + s2=prm%fcc_twinNucleationSlipPair(2,i) + if (tau(i) < dst%tau_r_tw(i,en)) then ! ToDo: correct? + Ndot0=(abs(dot_gamma_sl(s1))*(stt%rho_mob(s2,en)+stt%rho_dip(s2,en))+& + abs(dot_gamma_sl(s2))*(stt%rho_mob(s1,en)+stt%rho_dip(s1,en)))/& + (prm%L_tw*prm%b_sl(i))*& + (1.0_pReal-exp(-prm%V_cs/(kB*T)*(dst%tau_r_tw(i,en)-tau(i)))) + else + Ndot0=0.0_pReal + end if + else isFCC + Ndot0=prm%dot_N_0_tw(i) + endif isFCC + enddo - significantStress: where(tau > tol_math_check) - StressRatio_r = (dst%tau_hat_tw(:,en)/tau)**prm%r - dot_gamma_tw = prm%gamma_char * dst%V_tw(:,en) * Ndot0*exp(-StressRatio_r) - ddot_gamma_dtau = (dot_gamma_tw*prm%r/tau)*StressRatio_r - else where significantStress - dot_gamma_tw = 0.0_pReal - ddot_gamma_dtau = 0.0_pReal - end where significantStress + significantStress: where(tau > tol_math_check) + StressRatio_r = (dst%tau_hat_tw(:,en)/tau)**prm%r + dot_gamma_tw = prm%gamma_char * dst%V_tw(:,en) * Ndot0*exp(-StressRatio_r) + ddot_gamma_dtau = (dot_gamma_tw*prm%r/tau)*StressRatio_r + else where significantStress + dot_gamma_tw = 0.0_pReal + ddot_gamma_dtau = 0.0_pReal + end where significantStress end associate @@ -1009,36 +1011,37 @@ pure subroutine kinetics_tr(Mp,T,dot_gamma_sl,ph,en,& Ndot0, & stressRatio_s, & ddot_gamma_dtau - integer :: i,s1,s2 + + associate(prm => param(ph), stt => state(ph), dst => dependentState(ph)) - do i = 1, prm%sum_N_tr - tau(i) = math_tensordot(Mp,prm%P_tr(1:3,1:3,i)) - isFCC: if (prm%fccTwinTransNucleation) then - s1=prm%fcc_twinNucleationSlipPair(1,i) - s2=prm%fcc_twinNucleationSlipPair(2,i) - if (tau(i) < dst%tau_r_tr(i,en)) then ! ToDo: correct? - Ndot0=(abs(dot_gamma_sl(s1))*(stt%rho_mob(s2,en)+stt%rho_dip(s2,en))+& - abs(dot_gamma_sl(s2))*(stt%rho_mob(s1,en)+stt%rho_dip(s1,en)))/& ! ToDo: MD: it would be more consistent to use shearrates from state - (prm%L_tr*prm%b_sl(i))*& - (1.0_pReal-exp(-prm%V_cs/(kB*T)*(dst%tau_r_tr(i,en)-tau(i)))) ! P_ncs - else - Ndot0=0.0_pReal - end if - else isFCC - Ndot0=prm%dot_N_0_tr(i) - endif isFCC - enddo + do i = 1, prm%sum_N_tr + tau(i) = math_tensordot(Mp,prm%P_tr(1:3,1:3,i)) + isFCC: if (prm%fccTwinTransNucleation) then + s1=prm%fcc_twinNucleationSlipPair(1,i) + s2=prm%fcc_twinNucleationSlipPair(2,i) + if (tau(i) < dst%tau_r_tr(i,en)) then ! ToDo: correct? + Ndot0=(abs(dot_gamma_sl(s1))*(stt%rho_mob(s2,en)+stt%rho_dip(s2,en))+& + abs(dot_gamma_sl(s2))*(stt%rho_mob(s1,en)+stt%rho_dip(s1,en)))/& + (prm%L_tr*prm%b_sl(i))*& + (1.0_pReal-exp(-prm%V_cs/(kB*T)*(dst%tau_r_tr(i,en)-tau(i)))) + else + Ndot0=0.0_pReal + end if + else isFCC + Ndot0=prm%dot_N_0_tr(i) + endif isFCC + enddo - significantStress: where(tau > tol_math_check) - StressRatio_s = (dst%tau_hat_tr(:,en)/tau)**prm%s - dot_gamma_tr = dst%V_tr(:,en) * Ndot0*exp(-StressRatio_s) - ddot_gamma_dtau = (dot_gamma_tr*prm%s/tau)*StressRatio_s - else where significantStress - dot_gamma_tr = 0.0_pReal - ddot_gamma_dtau = 0.0_pReal - end where significantStress + significantStress: where(tau > tol_math_check) + StressRatio_s = (dst%tau_hat_tr(:,en)/tau)**prm%s + dot_gamma_tr = dst%V_tr(:,en) * Ndot0*exp(-StressRatio_s) + ddot_gamma_dtau = (dot_gamma_tr*prm%s/tau)*StressRatio_s + else where significantStress + dot_gamma_tr = 0.0_pReal + ddot_gamma_dtau = 0.0_pReal + end where significantStress end associate From 6f906647eb1198d6563851089f274d01b73a2daa Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sat, 24 Jul 2021 09:10:16 +0200 Subject: [PATCH 33/63] documenting hex --- PRIVATE | 2 +- .../mechanical/plastic/phenopowerlaw_Mg.yaml} | 33 ++++++++++++------- .../mechanical/plastic/phenopowerlaw_Ti.yaml | 5 ++- 3 files changed, 26 insertions(+), 14 deletions(-) rename examples/config/{Phase_Phenopowerlaw_Magnesium.yaml => phase/mechanical/plastic/phenopowerlaw_Mg.yaml} (65%) diff --git a/PRIVATE b/PRIVATE index 347932db3..020cf0621 160000 --- a/PRIVATE +++ b/PRIVATE @@ -1 +1 @@ -Subproject commit 347932db368ec5c0d9e9092ef106c225605d830a +Subproject commit 020cf06218fc4f14cb55ae0527972d12a57343bc diff --git a/examples/config/Phase_Phenopowerlaw_Magnesium.yaml b/examples/config/phase/mechanical/plastic/phenopowerlaw_Mg.yaml similarity index 65% rename from examples/config/Phase_Phenopowerlaw_Magnesium.yaml rename to examples/config/phase/mechanical/plastic/phenopowerlaw_Mg.yaml index 4d4b9eb8e..b29075173 100644 --- a/examples/config/Phase_Phenopowerlaw_Magnesium.yaml +++ b/examples/config/phase/mechanical/plastic/phenopowerlaw_Mg.yaml @@ -1,20 +1,29 @@ -N_sl: [3, 3, 0, 6, 0, 6] -N_tw: [6, 0, 0, 6] -h_0_tw-tw: 50.0e+6 -h_0_sl-sl: 500.0e+6 -h_0_tw-sl: 150.0e+6 -h_sl-sl: [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1] -h_tw-tw: [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1] -h_sl-tw: [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1] -h_tw-sl: [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1] -output: [xi_sl, xi_tw] type: phenopowerlaw -xi_0_sl: [10.e+6, 55.e+6, 0., 60.e+6, 0., 60.e+6] +references: + - F. Wang et al., + Acta Materialia 80:77-93, 2014, + https://doi.org/10.1016/j.actamat.2014.07.048 + +output: [xi_sl, xi_tw] + +N_sl: [3, 3, 0, 6, 0, 6] # basal, 1. prism, -, 1. pyr, -, 2. pyr +N_tw: [6, 0, 6] # tension, -, compression + +xi_0_sl: [10.e+6, 55.e+6, 0., 60.e+6, 0., 60.e+6] xi_inf_sl: [40.e+6, 135.e+6, 0., 150.e+6, 0., 150.e+6] -xi_0_tw: [40.e+6, 0., 0., 60.e+6] +xi_0_tw: [40.e+6, 0., 60.e+6] + a_sl: 2.25 dot_gamma_0_sl: 0.001 dot_gamma_0_tw: 0.001 n_sl: 20 n_tw: 20 f_sat_sl-tw: 10.0 + +h_0_sl-sl: 500.0e+6 +h_0_tw-tw: 50.0e+6 +h_0_tw-sl: 150.0e+6 +h_sl-sl: [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1] +h_tw-tw: [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1] +h_tw-sl: [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1] +h_sl-tw: [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1] diff --git a/examples/config/phase/mechanical/plastic/phenopowerlaw_Ti.yaml b/examples/config/phase/mechanical/plastic/phenopowerlaw_Ti.yaml index 5e05206df..804cf4541 100644 --- a/examples/config/phase/mechanical/plastic/phenopowerlaw_Ti.yaml +++ b/examples/config/phase/mechanical/plastic/phenopowerlaw_Ti.yaml @@ -7,14 +7,17 @@ references: Acta Materialia 132:598-610, 2017, https://doi.org/10.1016/j.actamat.2017.05.015 output: [gamma_sl] -N_sl: [3, 3, 0, 0, 12] + +N_sl: [3, 3, 0, 0, 12] # basal, 1. prism, -, -, 2. pyr n_sl: 20 a_sl: 2.0 dot_gamma_0_sl: 0.001 h_0_sl-sl: 200.e+6 + # C. Zambaldi et al.: xi_0_sl: [349.e+6, 150.e+6, 0.0, 0.0, 1107.e+6] xi_inf_sl: [568.e+6, 150.e+7, 0.0, 0.0, 3420.e+6] # L. Wang et al. : # xi_0_sl: [127.e+6, 96.e+6, 0.0, 0.0, 240.e+6] + h_sl-sl: [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1] From 931dc9955752ec99759c2aa36f9c137ee5ca0b89 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sat, 24 Jul 2021 11:21:07 +0200 Subject: [PATCH 34/63] not used climb formulation was updated a while ago --- PRIVATE | 2 +- .../config/phase/mechanical/plastic/dislotwin_IF-steel.yaml | 1 - src/phase_mechanical_plastic_dislotwin.f90 | 5 +---- 3 files changed, 2 insertions(+), 6 deletions(-) diff --git a/PRIVATE b/PRIVATE index 020cf0621..c656b6f08 160000 --- a/PRIVATE +++ b/PRIVATE @@ -1 +1 @@ -Subproject commit 020cf06218fc4f14cb55ae0527972d12a57343bc +Subproject commit c656b6f08489756c9ee6a6e1a62858c8b7836f10 diff --git a/examples/config/phase/mechanical/plastic/dislotwin_IF-steel.yaml b/examples/config/phase/mechanical/plastic/dislotwin_IF-steel.yaml index e48d5437f..362797d0b 100644 --- a/examples/config/phase/mechanical/plastic/dislotwin_IF-steel.yaml +++ b/examples/config/phase/mechanical/plastic/dislotwin_IF-steel.yaml @@ -19,6 +19,5 @@ i_sl: [23.3, 23.3] D_a: 7.4 # C_anni B: [0.001, 0.001] h_sl-sl: [0.1, 0.72, 0.1, 0.053, 0.053, 0.073, 0.137, 0.72, 0.72, 0.053, 0.053, 0.053, 0.053, 0.073, 0.073, 0.073, 0.073, 0.073, 0.073, 0.137, 0.073, 0.073, 0.137, 0.073] -D_0: 4.0e-05 Q_cl: 5.4e-19 # no recovery! D: 40.e-6 # estimated diff --git a/src/phase_mechanical_plastic_dislotwin.f90 b/src/phase_mechanical_plastic_dislotwin.f90 index 1d6fadd13..2cc6f4af8 100644 --- a/src/phase_mechanical_plastic_dislotwin.f90 +++ b/src/phase_mechanical_plastic_dislotwin.f90 @@ -16,7 +16,6 @@ submodule(phase:plastic) dislotwin real(pReal) :: & mu = 1.0_pReal, & !< equivalent shear modulus nu = 1.0_pReal, & !< equivalent shear Poisson's ratio - D_0 = 1.0_pReal, & !< prefactor for self-diffusion coefficient Q_cl = 1.0_pReal, & !< activation energy for dislocation climb omega = 1.0_pReal, & !< frequency factor for dislocation climb D = 1.0_pReal, & !< grain size @@ -214,7 +213,6 @@ module function plastic_dislotwin_init() result(myPlasticity) prm%B = pl%get_as1dFloat('B', requiredSize=size(N_sl), & defaultVal=[(0.0_pReal, i=1,size(N_sl))]) - prm%D_0 = pl%get_asFloat('D_0') prm%Q_cl = pl%get_asFloat('Q_cl') prm%ExtendedDislocations = pl%get_asBool('extend_dislocations',defaultVal = .false.) @@ -230,7 +228,7 @@ module function plastic_dislotwin_init() result(myPlasticity) rho_dip_0 = math_expand(rho_dip_0, N_sl) prm%v_0 = math_expand(prm%v_0, N_sl) prm%b_sl = math_expand(prm%b_sl, N_sl) - prm%Q_sl = math_expand(prm%Q_sl, N_sl) + prm%Q_sl = math_expand(prm%Q_sl, N_sl) prm%i_sl = math_expand(prm%i_sl, N_sl) prm%p = math_expand(prm%p, N_sl) prm%q = math_expand(prm%q, N_sl) @@ -239,7 +237,6 @@ module function plastic_dislotwin_init() result(myPlasticity) prm%d_caron = pl%get_asFloat('D_a') * prm%b_sl ! sanity checks - if ( prm%D_0 <= 0.0_pReal) extmsg = trim(extmsg)//' D_0' if ( prm%Q_cl <= 0.0_pReal) extmsg = trim(extmsg)//' Q_cl' if (any(rho_mob_0 < 0.0_pReal)) extmsg = trim(extmsg)//' rho_mob_0' if (any(rho_dip_0 < 0.0_pReal)) extmsg = trim(extmsg)//' rho_dip_0' From d6ce721a2516d1bbb8aa7f3253ae49bc0d061aff Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sat, 24 Jul 2021 12:10:59 +0200 Subject: [PATCH 35/63] need to check type of leaf otherwise, using a [list] where a scalar is expected results in a crash, not in a meaningful error message --- src/YAML_types.f90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/YAML_types.f90 b/src/YAML_types.f90 index 45907fbae..e11be1275 100644 --- a/src/YAML_types.f90 +++ b/src/YAML_types.f90 @@ -849,7 +849,7 @@ function tNode_get_byKey_as1dFloat(self,k,defaultVal,requiredSize) result(nodeAs if (self%contains(k)) then node => self%get(k) - select type(self) + select type(node) class is(tList) list => node%asList() nodeAs1dFloat = list%as1dFloat() From f22f30e05da9008998a6c9316c470c8ac621c93e Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sat, 24 Jul 2021 14:27:00 +0200 Subject: [PATCH 36/63] same functionality as for 1D --- src/YAML_types.f90 | 7 ++++++- 1 file changed, 6 insertions(+), 1 deletion(-) diff --git a/src/YAML_types.f90 b/src/YAML_types.f90 index e11be1275..8f81923a8 100644 --- a/src/YAML_types.f90 +++ b/src/YAML_types.f90 @@ -872,11 +872,12 @@ end function tNode_get_byKey_as1dFloat !-------------------------------------------------------------------------------------------------- !> @brief Access by key and convert to float array (2D) !-------------------------------------------------------------------------------------------------- -function tNode_get_byKey_as2dFloat(self,k,defaultVal) result(nodeAs2dFloat) +function tNode_get_byKey_as2dFloat(self,k,defaultVal,requiredShape) result(nodeAs2dFloat) class(tNode), intent(in), target :: self character(len=*), intent(in) :: k real(pReal), intent(in), dimension(:,:), optional :: defaultVal + integer, intent(in), dimension(2), optional :: requiredShape real(pReal), dimension(:,:), allocatable :: nodeAs2dFloat @@ -898,6 +899,10 @@ function tNode_get_byKey_as2dFloat(self,k,defaultVal) result(nodeAs2dFloat) call IO_error(143,ext_msg=k) endif + if (present(requiredShape)) then + if (any(requiredShape /= shape(nodeAs2dFloat))) call IO_error(146,ext_msg=k) + endif + end function tNode_get_byKey_as2dFloat From b3f5e1223235043c9c736130058b1b564d19f194 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sat, 24 Jul 2021 11:43:14 +0200 Subject: [PATCH 37/63] using burgers vector of tw/tr system seems to make more sense here otherwise twinning will not work in many cases. Matching number is only required for nucleation of tw/tr --- PRIVATE | 2 +- src/phase_mechanical_plastic_dislotwin.f90 | 16 +++++++--------- 2 files changed, 8 insertions(+), 10 deletions(-) diff --git a/PRIVATE b/PRIVATE index c656b6f08..5049664c5 160000 --- a/PRIVATE +++ b/PRIVATE @@ -1 +1 @@ -Subproject commit c656b6f08489756c9ee6a6e1a62858c8b7836f10 +Subproject commit 5049664c5cac3cc1571c7b61f3345f1ba8d627f6 diff --git a/src/phase_mechanical_plastic_dislotwin.f90 b/src/phase_mechanical_plastic_dislotwin.f90 index 2cc6f4af8..4242d0e7a 100644 --- a/src/phase_mechanical_plastic_dislotwin.f90 +++ b/src/phase_mechanical_plastic_dislotwin.f90 @@ -205,7 +205,7 @@ module function plastic_dislotwin_init() result(myPlasticity) rho_dip_0 = pl%get_as1dFloat('rho_dip_0', requiredSize=size(N_sl)) prm%v_0 = pl%get_as1dFloat('v_0', requiredSize=size(N_sl)) prm%b_sl = pl%get_as1dFloat('b_sl', requiredSize=size(N_sl)) - prm%Q_sl = pl%get_as1dFloat('Q_sl', requiredSize=size(N_sl)) + prm%Q_sl = pl%get_as1dFloat('Q_sl', requiredSize=size(N_sl)) prm%i_sl = pl%get_as1dFloat('i_sl', requiredSize=size(N_sl)) prm%p = pl%get_as1dFloat('p_sl', requiredSize=size(N_sl)) prm%q = pl%get_as1dFloat('q_sl', requiredSize=size(N_sl)) @@ -758,19 +758,17 @@ module subroutine dislotwin_dependentState(T,ph,en) dst%tau_pass(:,en) = prm%mu*prm%b_sl* sqrt(matmul(prm%h_sl_sl,stt%rho_mob(:,en)+stt%rho_dip(:,en))) !* threshold stress for growing twin/martensite - if(prm%sum_N_tw == prm%sum_N_sl) & - dst%tau_hat_tw(:,en) = Gamma/(3.0_pReal*prm%b_tw) & - + 3.0_pReal*prm%b_tw*prm%mu/(prm%L_tw*prm%b_sl) ! slip Burgers here correct? - if(prm%sum_N_tr == prm%sum_N_sl) & - dst%tau_hat_tr(:,en) = Gamma/(3.0_pReal*prm%b_tr) & - + 3.0_pReal*prm%b_tr*prm%mu/(prm%L_tr*prm%b_sl) & ! slip Burgers here correct? - + prm%h*prm%delta_G/ (3.0_pReal*prm%b_tr) + dst%tau_hat_tw(:,en) = Gamma/(3.0_pReal*prm%b_tw) & + + 3.0_pReal*prm%b_tw*prm%mu/(prm%L_tw*prm%b_tw) + dst%tau_hat_tr(:,en) = Gamma/(3.0_pReal*prm%b_tr) & + + 3.0_pReal*prm%b_tr*prm%mu/(prm%L_tr*prm%b_tr) & + + prm%h*prm%delta_G/(3.0_pReal*prm%b_tr) dst%V_tw(:,en) = (PI/4.0_pReal)*prm%t_tw*dst%Lambda_tw(:,en)**2.0_pReal dst%V_tr(:,en) = (PI/4.0_pReal)*prm%t_tr*dst%Lambda_tr(:,en)**2.0_pReal - x0 = prm%mu*prm%b_tw**2.0_pReal/(Gamma*8.0_pReal*PI)*(2.0_pReal+prm%nu)/(1.0_pReal-prm%nu) ! ToDo: In the paper, this is the Burgers vector for slip and is the same for twin and trans + x0 = prm%mu*prm%b_tw**2.0_pReal/(Gamma*8.0_pReal*PI)*(2.0_pReal+prm%nu)/(1.0_pReal-prm%nu) ! ToDo: In the paper, this is the Burgers vector for slip dst%tau_r_tw(:,en) = prm%mu*prm%b_tw/(2.0_pReal*PI)*(1.0_pReal/(x0+prm%x_c_tw)+cos(pi/3.0_pReal)/x0) x0 = prm%mu*prm%b_tr**2.0_pReal/(Gamma*8.0_pReal*PI)*(2.0_pReal+prm%nu)/(1.0_pReal-prm%nu) ! ToDo: In the paper, this is the Burgers vector for slip From 3bb5ae3d9ec03059e7e18b96544ab431984e1883 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sun, 25 Jul 2021 13:20:39 +0200 Subject: [PATCH 38/63] concurrent/parallel execution is possible here --- src/math.f90 | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/math.f90 b/src/math.f90 index 36fe91bac..aeab47093 100644 --- a/src/math.f90 +++ b/src/math.f90 @@ -736,9 +736,9 @@ pure function math_3333to99(m3333) integer :: i,j - do i=1,9; do j=1,9 + do concurrent(i=1:9, j=1:9) math_3333to99(i,j) = m3333(MAPPLAIN(1,i),MAPPLAIN(2,i),MAPPLAIN(1,j),MAPPLAIN(2,j)) - enddo; enddo + enddo end function math_3333to99 @@ -753,9 +753,9 @@ pure function math_99to3333(m99) integer :: i,j - do i=1,9; do j=1,9 + do concurrent(i=1:9, j=1:9) math_99to3333(MAPPLAIN(1,i),MAPPLAIN(2,i),MAPPLAIN(1,j),MAPPLAIN(2,j)) = m99(i,j) - enddo; enddo + enddo end function math_99to3333 From 31a40006553b87cea43090ac43702606dd3ddbad Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sun, 25 Jul 2021 13:32:11 +0200 Subject: [PATCH 39/63] no need for arbitrary dimension not sure if correct, not used at all (even identity4th for 3 dim is not used, but now at least tested) --- src/math.f90 | 21 +++++++++++---------- 1 file changed, 11 insertions(+), 10 deletions(-) diff --git a/src/math.f90 b/src/math.f90 index aeab47093..110af5128 100644 --- a/src/math.f90 +++ b/src/math.f90 @@ -233,7 +233,7 @@ end function math_range !-------------------------------------------------------------------------------------------------- -!> @brief second rank identity tensor of specified dimension +!> @brief Rank two identity tensor of specified dimension. !-------------------------------------------------------------------------------------------------- pure function math_eye(d) @@ -250,20 +250,18 @@ end function math_eye !-------------------------------------------------------------------------------------------------- -!> @brief symmetric fourth rank identity tensor of specified dimension +!> @brief Symmetric rank four identity tensor. ! from http://en.wikipedia.org/wiki/Tensor_derivative_(continuum_mechanics)#Derivative_of_a_second-order_tensor_with_respect_to_itself !-------------------------------------------------------------------------------------------------- -pure function math_identity4th(d) +pure function math_identity4th() + + real(pReal), dimension(3,3,3,3) :: math_identity4th - integer, intent(in) :: d !< tensor dimension integer :: i,j,k,l - real(pReal), dimension(d,d,d,d) :: math_identity4th - real(pReal), dimension(d,d) :: identity2nd - identity2nd = math_eye(d) - do i=1,d; do j=1,d; do k=1,d; do l=1,d - math_identity4th(i,j,k,l) = 0.5_pReal & - *(identity2nd(i,k)*identity2nd(j,l)+identity2nd(i,l)*identity2nd(j,k)) + + do i=1,3; do j=1,3; do k=1,3; do l=1,3 + math_identity4th(i,j,k,l) = 0.5_pReal*(math_I3(i,k)*math_I3(j,l)+math_I3(i,l)*math_I3(j,k)) enddo; enddo; enddo; enddo end function math_identity4th @@ -1242,6 +1240,9 @@ subroutine selfTest if(dNeq(math_det33(math_symmetric33(t33)),math_detSym33(math_symmetric33(t33)),tol=1.0e-12_pReal)) & error stop 'math_det33/math_detSym33' + if(any(dNeq0(t33+transpose(t33)-math_mul3333xx33(math_identity4th(),t33+transpose(t33))))) & + error stop 'math_mul3333xx33/math_identity4th' + if(any(dNeq0(math_eye(3),math_inv33(math_I3)))) & error stop 'math_inv33(math_I3)' From d19ab4c4f62977bcfafccda895656f35774ba126 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sun, 25 Jul 2021 13:38:56 +0200 Subject: [PATCH 40/63] co concurrent possible here note: do concurrent cannot be used for the double loops for 66 to 3333! --- PRIVATE | 2 +- src/math.f90 | 30 +++++++++++++++++++----------- 2 files changed, 20 insertions(+), 12 deletions(-) diff --git a/PRIVATE b/PRIVATE index 174ecac2d..d92b030e5 160000 --- a/PRIVATE +++ b/PRIVATE @@ -1 +1 @@ -Subproject commit 174ecac2d3ab7596bdb60184d6bb9e1a52cb7378 +Subproject commit d92b030e5777e718c77edc2e1e93abfa0981b024 diff --git a/src/math.f90 b/src/math.f90 index 110af5128..89ee7abfe 100644 --- a/src/math.f90 +++ b/src/math.f90 @@ -260,9 +260,9 @@ pure function math_identity4th() integer :: i,j,k,l - do i=1,3; do j=1,3; do k=1,3; do l=1,3 + do concurrent(i=1:3, j=1:3, k=1:3, l=1:3) math_identity4th(i,j,k,l) = 0.5_pReal*(math_I3(i,k)*math_I3(j,l)+math_I3(i,l)*math_I3(j,k)) - enddo; enddo; enddo; enddo + enddo end function math_identity4th @@ -329,9 +329,10 @@ pure function math_outer(A,B) real(pReal), dimension(size(A,1),size(B,1)) :: math_outer integer :: i,j - do i=1,size(A,1); do j=1,size(B,1) + + do concurrent(i=1:size(A,1), j=1:size(B,1)) math_outer(i,j) = A(i)*B(j) - enddo; enddo + enddo end function math_outer @@ -371,9 +372,10 @@ pure function math_mul3333xx33(A,B) real(pReal), dimension(3,3) :: math_mul3333xx33 integer :: i,j - do i=1,3; do j=1,3 + + do concurrent(i=1:3, j=1:3) math_mul3333xx33(i,j) = sum(A(i,j,1:3,1:3)*B(1:3,1:3)) - enddo; enddo + enddo end function math_mul3333xx33 @@ -388,9 +390,9 @@ pure function math_mul3333xx3333(A,B) real(pReal), dimension(3,3,3,3), intent(in) :: B real(pReal), dimension(3,3,3,3) :: math_mul3333xx3333 - do i=1,3; do j=1,3; do k=1,3; do l=1,3 + do concurrent(i=1:3, j=1:3, k=1:3, l=1:3) math_mul3333xx3333(i,j,k,l) = sum(A(i,j,1:3,1:3)*B(1:3,1:3,k,l)) - enddo; enddo; enddo; enddo + enddo end function math_mul3333xx3333 @@ -710,6 +712,7 @@ pure function math_6toSym33(v6,weighted) real(pReal), dimension(6) :: w integer :: i + if(present(weighted)) then w = merge(INVNRMMANDEL,1.0_pReal,weighted) else @@ -734,6 +737,7 @@ pure function math_3333to99(m3333) integer :: i,j + do concurrent(i=1:9, j=1:9) math_3333to99(i,j) = m3333(MAPPLAIN(1,i),MAPPLAIN(2,i),MAPPLAIN(1,j),MAPPLAIN(2,j)) enddo @@ -751,6 +755,7 @@ pure function math_99to3333(m99) integer :: i,j + do concurrent(i=1:9, j=1:9) math_99to3333(MAPPLAIN(1,i),MAPPLAIN(2,i),MAPPLAIN(1,j),MAPPLAIN(2,j)) = m99(i,j) enddo @@ -773,15 +778,16 @@ pure function math_sym3333to66(m3333,weighted) real(pReal), dimension(6) :: w integer :: i,j + if(present(weighted)) then w = merge(NRMMANDEL,1.0_pReal,weighted) else w = NRMMANDEL endif - do i=1,6; do j=1,6 + do concurrent(i=1:6, j=1:6) math_sym3333to66(i,j) = w(i)*w(j)*m3333(MAPNYE(1,i),MAPNYE(2,i),MAPNYE(1,j),MAPNYE(2,j)) - enddo; enddo + enddo end function math_sym3333to66 @@ -801,6 +807,7 @@ pure function math_66toSym3333(m66,weighted) real(pReal), dimension(6) :: w integer :: i,j + if(present(weighted)) then w = merge(INVNRMMANDEL,1.0_pReal,weighted) else @@ -826,6 +833,7 @@ pure function math_Voigt66to3333(m66) real(pReal), dimension(6,6), intent(in) :: m66 !< 6x6 matrix integer :: i,j + do i=1,6; do j=1, 6 math_Voigt66to3333(MAPVOIGT(1,i),MAPVOIGT(2,i),MAPVOIGT(1,j),MAPVOIGT(2,j)) = m66(i,j) math_Voigt66to3333(MAPVOIGT(2,i),MAPVOIGT(1,i),MAPVOIGT(1,j),MAPVOIGT(2,j)) = m66(i,j) @@ -1240,7 +1248,7 @@ subroutine selfTest if(dNeq(math_det33(math_symmetric33(t33)),math_detSym33(math_symmetric33(t33)),tol=1.0e-12_pReal)) & error stop 'math_det33/math_detSym33' - if(any(dNeq0(t33+transpose(t33)-math_mul3333xx33(math_identity4th(),t33+transpose(t33))))) & + if(any(dNeq(t33+transpose(t33),math_mul3333xx33(math_identity4th(),t33+transpose(t33))))) & error stop 'math_mul3333xx33/math_identity4th' if(any(dNeq0(math_eye(3),math_inv33(math_I3)))) & From e56f2e09a46ff105524946ac0ca8e9232d707ec7 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Mon, 26 Jul 2021 10:22:53 +0200 Subject: [PATCH 41/63] trivial and never used even the more commoly used fucntions for 3x3 matrices are not all in use --- src/math.f90 | 13 ------------- 1 file changed, 13 deletions(-) diff --git a/src/math.f90 b/src/math.f90 index 36fe91bac..a48f44b93 100644 --- a/src/math.f90 +++ b/src/math.f90 @@ -545,19 +545,6 @@ pure function math_symmetric33(m) end function math_symmetric33 -!-------------------------------------------------------------------------------------------------- -!> @brief symmetrize a 6x6 matrix -!-------------------------------------------------------------------------------------------------- -pure function math_symmetric66(m) - - real(pReal), dimension(6,6) :: math_symmetric66 - real(pReal), dimension(6,6), intent(in) :: m - - math_symmetric66 = 0.5_pReal * (m + transpose(m)) - -end function math_symmetric66 - - !-------------------------------------------------------------------------------------------------- !> @brief skew part of a 3x3 matrix !-------------------------------------------------------------------------------------------------- From ce4052a49d9da97e75c68d78083beed8a81166b7 Mon Sep 17 00:00:00 2001 From: Test User Date: Mon, 26 Jul 2021 13:30:08 +0200 Subject: [PATCH 42/63] [skip ci] updated version information after successful test of v3.0.0-alpha4-141-g18dbda784 --- python/damask/VERSION | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/python/damask/VERSION b/python/damask/VERSION index 026bbcf9f..e1c11a09b 100644 --- a/python/damask/VERSION +++ b/python/damask/VERSION @@ -1 +1 @@ -v3.0.0-alpha4-137-gb69b85754 +v3.0.0-alpha4-141-g18dbda784 From 93814054d51f6ccc465f51208a874ca2216e79e9 Mon Sep 17 00:00:00 2001 From: Test User Date: Mon, 26 Jul 2021 15:02:14 +0200 Subject: [PATCH 43/63] [skip ci] updated version information after successful test of v3.0.0-alpha4-171-g6a9c892d7 --- python/damask/VERSION | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/python/damask/VERSION b/python/damask/VERSION index 026bbcf9f..c4c59cb81 100644 --- a/python/damask/VERSION +++ b/python/damask/VERSION @@ -1 +1 @@ -v3.0.0-alpha4-137-gb69b85754 +v3.0.0-alpha4-171-g6a9c892d7 From 04ab83a11fa0d56cb689b8d7f1a1b2d6481abdab Mon Sep 17 00:00:00 2001 From: Test User Date: Mon, 26 Jul 2021 18:16:19 +0200 Subject: [PATCH 44/63] [skip ci] updated version information after successful test of v3.0.0-alpha4-182-gac6d31b1f --- python/damask/VERSION | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/python/damask/VERSION b/python/damask/VERSION index e1c11a09b..13d57229b 100644 --- a/python/damask/VERSION +++ b/python/damask/VERSION @@ -1 +1 @@ -v3.0.0-alpha4-141-g18dbda784 +v3.0.0-alpha4-182-gac6d31b1f From 835c9b80024ac6878c878f6efaa1afc7e7c01888 Mon Sep 17 00:00:00 2001 From: Test User Date: Mon, 26 Jul 2021 19:54:00 +0200 Subject: [PATCH 45/63] [skip ci] updated version information after successful test of v3.0.0-alpha4-194-g012f09caa --- python/damask/VERSION | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/python/damask/VERSION b/python/damask/VERSION index e1c11a09b..c76dacadf 100644 --- a/python/damask/VERSION +++ b/python/damask/VERSION @@ -1 +1 @@ -v3.0.0-alpha4-141-g18dbda784 +v3.0.0-alpha4-194-g012f09caa From 58c1f5fdc1f4ac58b5456b9282376e2a0b174acc Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Mon, 26 Jul 2021 22:39:23 +0200 Subject: [PATCH 46/63] write string, only from one MPI process needed to store simulation setup. I don't see a situation where strings are distributed over multiple processes --- src/HDF5_utilities.f90 | 34 +++++++++++++++++++++++++++++++++- 1 file changed, 33 insertions(+), 1 deletion(-) diff --git a/src/HDF5_utilities.f90 b/src/HDF5_utilities.f90 index c981fad53..7de6be431 100644 --- a/src/HDF5_utilities.f90 +++ b/src/HDF5_utilities.f90 @@ -85,6 +85,7 @@ module HDF5_utilities HDF5_utilities_init, & HDF5_read, & HDF5_write, & + HDF5_write_str, & HDF5_addAttribute, & HDF5_addGroup, & HDF5_openGroup, & @@ -1467,6 +1468,37 @@ subroutine HDF5_write_real7(dataset,loc_id,datasetName,parallel) end subroutine HDF5_write_real7 +subroutine HDF5_write_str(dataset,loc_id,datasetName) + + character(len=*), intent(in) :: dataset + integer(HID_T), intent(in) :: loc_id + character(len=*), intent(in) :: datasetName !< name of the dataset in the file + + INTEGER(HID_T) :: filetype_id, space_id, dataset_id + INTEGER :: hdferr + + character(len=len_trim(dataset)+1,kind=C_CHAR), dimension(1), target :: dataset_ + type(C_PTR), target, dimension(1) :: ptr + + + dataset_(1) = trim(dataset)//C_NULL_CHAR + ptr(1) = c_loc(dataset_(1)) + + call h5tcopy_f(H5T_STRING, filetype_id, hdferr) + call h5tset_size_f(filetype_id, int(len(dataset_),HSIZE_T), hdferr) + + call h5screate_f(H5S_SCALAR_F, space_id, hdferr) + call h5dcreate_f(loc_id, datasetName, H5T_STRING, space_id, dataset_id, hdferr) + + call h5dwrite_f(dataset_id, H5T_STRING, c_loc(ptr), hdferr); + + call h5dclose_f(dataset_id, hdferr) + call h5sclose_f(space_id, hdferr) + call h5tclose_f(filetype_id, hdferr) + +end subroutine HDF5_write_str + + !-------------------------------------------------------------------------------------------------- !> @brief write dataset of type integer with 1 dimension !-------------------------------------------------------------------------------------------------- @@ -1872,7 +1904,7 @@ subroutine initialize_write(dset_id, filespace_id, memspace_id, plist_id, & integer(HSIZE_T), parameter :: chunkSize = 1024_HSIZE_T**2/8_HSIZE_T !------------------------------------------------------------------------------------------------- -! creating a property list for transfer properties (is collective when reading in parallel) +! creating a property list for transfer properties (is collective when writing in parallel) call h5pcreate_f(H5P_DATASET_XFER_F, plist_id, hdferr) if(hdferr < 0) error stop 'HDF5 error' #ifdef PETSC From 820d590b6b8a5096226f65029111f2475089119a Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Tue, 27 Jul 2021 07:25:44 +0200 Subject: [PATCH 47/63] write string to result requires to open file without MPI support --- src/HDF5_utilities.f90 | 9 +++++++-- src/results.f90 | 27 +++++++++++++++++++++++++-- 2 files changed, 32 insertions(+), 4 deletions(-) diff --git a/src/HDF5_utilities.f90 b/src/HDF5_utilities.f90 index 7de6be431..bd88feee2 100644 --- a/src/HDF5_utilities.f90 +++ b/src/HDF5_utilities.f90 @@ -129,10 +129,11 @@ end subroutine HDF5_utilities_init !-------------------------------------------------------------------------------------------------- !> @brief open and initializes HDF5 output file !-------------------------------------------------------------------------------------------------- -integer(HID_T) function HDF5_openFile(fileName,mode) +integer(HID_T) function HDF5_openFile(fileName,mode,parallel) character(len=*), intent(in) :: fileName character, intent(in), optional :: mode + logical, intent(in), optional :: parallel character :: m integer(HID_T) :: plist_id @@ -149,7 +150,11 @@ integer(HID_T) function HDF5_openFile(fileName,mode) if(hdferr < 0) error stop 'HDF5 error' #ifdef PETSC - call h5pset_fapl_mpio_f(plist_id, PETSC_COMM_WORLD, MPI_INFO_NULL, hdferr) + if (present(parallel)) then + if (parallel) call h5pset_fapl_mpio_f(plist_id, PETSC_COMM_WORLD, MPI_INFO_NULL, hdferr) + else + call h5pset_fapl_mpio_f(plist_id, PETSC_COMM_WORLD, MPI_INFO_NULL, hdferr) + endif if(hdferr < 0) error stop 'HDF5 error' #endif diff --git a/src/results.f90 b/src/results.f90 index 94625a4b9..ef6298151 100644 --- a/src/results.f90 +++ b/src/results.f90 @@ -52,6 +52,7 @@ module results results_openGroup, & results_closeGroup, & results_writeDataset, & + results_writeDataset_str, & results_setLink, & results_addAttribute, & results_removeLink, & @@ -90,9 +91,12 @@ end subroutine results_init !-------------------------------------------------------------------------------------------------- !> @brief opens the results file to append data !-------------------------------------------------------------------------------------------------- -subroutine results_openJobFile +subroutine results_openJobFile(parallel) - resultsFile = HDF5_openFile(getSolverJobName()//'.hdf5','a') + logical, intent(in), optional :: parallel + + + resultsFile = HDF5_openFile(getSolverJobName()//'.hdf5','a',parallel) end subroutine results_openJobFile @@ -297,6 +301,25 @@ subroutine results_removeLink(link) end subroutine results_removeLink +!-------------------------------------------------------------------------------------------------- +!> @brief Store string dataset. +!> @details Not collective, must be called by one process at at time. +!-------------------------------------------------------------------------------------------------- +subroutine results_writeDataset_str(dataset,group,label,description) + + character(len=*), intent(in) :: label,group,description,dataset + + integer(HID_T) :: groupHandle + + + groupHandle = results_openGroup(group) + call HDF5_write_str(dataset,groupHandle,label) + call executionStamp(group//'/'//label,description) + call HDF5_closeGroup(groupHandle) + +end subroutine results_writeDataset_str + + !-------------------------------------------------------------------------------------------------- !> @brief Store real scalar dataset with associated metadata. !-------------------------------------------------------------------------------------------------- From f20302e0eed91f63043a8f97971ee26365879d56 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Tue, 27 Jul 2021 07:27:18 +0200 Subject: [PATCH 48/63] early initialization of results need to be done before config if config writes to it --- src/CPFEM.f90 | 27 ++++++++++++++------------- src/CPFEM2.f90 | 34 ++++++++++++++++++---------------- 2 files changed, 32 insertions(+), 29 deletions(-) diff --git a/src/CPFEM.f90 b/src/CPFEM.f90 index aa532859a..4123af37e 100644 --- a/src/CPFEM.f90 +++ b/src/CPFEM.f90 @@ -4,22 +4,23 @@ !> @brief CPFEM engine !-------------------------------------------------------------------------------------------------- module CPFEM + use DAMASK_interface use prec - use math - use rotations + use IO use YAML_types use YAML_parse - use discretization_marc - use material - use config - use homogenization - use IO - use discretization - use DAMASK_interface use HDF5_utilities use results + use config + use math + use rotations use lattice + use material use phase + use homogenization + + use discretization + use discretization_marc implicit none private @@ -68,7 +69,7 @@ contains !-------------------------------------------------------------------------------------------------- -!> @brief call all module initializations +!> @brief Initialize all modules. !-------------------------------------------------------------------------------------------------- subroutine CPFEM_initAll @@ -77,13 +78,13 @@ subroutine CPFEM_initAll call IO_init call YAML_types_init call YAML_parse_init + call HDF5_utilities_init + call results_init(.false.) call config_init call math_init call rotations_init - call HDF5_utilities_init - call results_init(.false.) - call discretization_marc_init call lattice_init + call discretization_marc_init call material_init(.false.) call phase_init call homogenization_init diff --git a/src/CPFEM2.f90 b/src/CPFEM2.f90 index ea5820852..2bb8420b9 100644 --- a/src/CPFEM2.f90 +++ b/src/CPFEM2.f90 @@ -4,28 +4,29 @@ !> @brief needs a good name and description !-------------------------------------------------------------------------------------------------- module CPFEM2 - use prec use parallelization - use config - use math - use rotations + use DAMASK_interface + use prec + use IO use YAML_types use YAML_parse - use material - use lattice - use IO - use base64 - use DAMASK_interface - use discretization use HDF5 use HDF5_utilities use results - use homogenization + use config + use math + use rotations + use lattice + use material use phase + use homogenization + + use discretization #if defined(MESH) use FEM_quadrature use discretization_mesh #elif defined(GRID) + use base64 use discretization_grid #endif @@ -36,7 +37,7 @@ contains !-------------------------------------------------------------------------------------------------- -!> @brief call all module initializations +!> @brief Initialize all modules. !-------------------------------------------------------------------------------------------------- subroutine CPFEM_initAll @@ -44,18 +45,19 @@ subroutine CPFEM_initAll call DAMASK_interface_init ! Spectral and FEM interface to commandline call prec_init call IO_init - call base64_init -#ifdef MESH +#if defined(MESH) call FEM_quadrature_init +#elif defined(GRID) + call base64_init #endif call YAML_types_init call YAML_parse_init + call HDF5_utilities_init + call results_init(restart=interface_restartInc>0) call config_init call math_init call rotations_init call lattice_init - call HDF5_utilities_init - call results_init(restart=interface_restartInc>0) #if defined(MESH) call discretization_mesh_init(restart=interface_restartInc>0) #elif defined(GRID) From 3b06498c2f1b3d9d563e7675a7f25ab9032b87bf Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Tue, 27 Jul 2021 07:28:35 +0200 Subject: [PATCH 49/63] make results self-contained for reproducibility ToDo: same functionality for load and geom --- src/config.f90 | 45 ++++++++++++++++++++++++++++++++++----------- src/results.f90 | 2 ++ 2 files changed, 36 insertions(+), 11 deletions(-) diff --git a/src/config.f90 b/src/config.f90 index 7d75cb444..8ee25182f 100644 --- a/src/config.f90 +++ b/src/config.f90 @@ -8,7 +8,8 @@ module config use IO use YAML_parse use YAML_types - + use results + use parallelization implicit none private @@ -31,6 +32,7 @@ subroutine config_init print'(/,a)', ' <<<+- config init -+>>>'; flush(IO_STDOUT) + call parse_material call parse_numerics call parse_debug @@ -41,14 +43,21 @@ end subroutine config_init !-------------------------------------------------------------------------------------------------- !> @brief Read material.yaml or .yaml. !-------------------------------------------------------------------------------------------------- -subroutine parse_material +subroutine parse_material() logical :: fileExists + character(len=:), allocatable :: fileContent inquire(file='material.yaml',exist=fileExists) if(.not. fileExists) call IO_error(100,ext_msg='material.yaml') print*, 'reading material.yaml'; flush(IO_STDOUT) + fileContent = IO_read('material.yaml') + if (worldrank == 0) then + call results_openJobFile(parallel=.false.) + call results_writeDataset_str(fileContent,'setup','material.yaml','DAMASK main configuration') + call results_closeJobFile + endif config_material => YAML_parse_file('material.yaml') end subroutine parse_material @@ -57,15 +66,22 @@ end subroutine parse_material !-------------------------------------------------------------------------------------------------- !> @brief Read numerics.yaml. !-------------------------------------------------------------------------------------------------- -subroutine parse_numerics +subroutine parse_numerics() - logical :: fexist + logical :: fileExists + character(len=:), allocatable :: fileContent config_numerics => emptyDict - inquire(file='numerics.yaml', exist=fexist) - if (fexist) then + inquire(file='numerics.yaml', exist=fileExists) + if (fileExists) then print*, 'reading numerics.yaml'; flush(IO_STDOUT) + fileContent = IO_read('numerics.yaml') + if (worldrank == 0) then + call results_openJobFile(parallel=.false.) + call results_writeDataset_str(fileContent,'setup','numerics.yaml','numerics configuration (optional)') + call results_closeJobFile + endif config_numerics => YAML_parse_file('numerics.yaml') endif @@ -75,17 +91,24 @@ end subroutine parse_numerics !-------------------------------------------------------------------------------------------------- !> @brief Read debug.yaml. !-------------------------------------------------------------------------------------------------- -subroutine parse_debug +subroutine parse_debug() - logical :: fexist + logical :: fileExists + character(len=:), allocatable :: fileContent config_debug => emptyDict - inquire(file='debug.yaml', exist=fexist) - fileExists: if (fexist) then + inquire(file='debug.yaml', exist=fileExists) + if (fileExists) then print*, 'reading debug.yaml'; flush(IO_STDOUT) + fileContent = IO_read('debug.yaml') + if (worldrank == 0) then + call results_openJobFile(parallel=.false.) + call results_writeDataset_str(fileContent,'setup','debug.yaml','debug configuration (optional)') + call results_closeJobFile + endif config_debug => YAML_parse_file('debug.yaml') - endif fileExists + endif end subroutine parse_debug diff --git a/src/results.f90 b/src/results.f90 index ef6298151..8b04800ef 100644 --- a/src/results.f90 +++ b/src/results.f90 @@ -82,6 +82,8 @@ subroutine results_init(restart) call results_addAttribute('call',trim(commandLine)) call results_closeGroup(results_addGroup('cell_to')) call results_addAttribute('description','mappings to place data in space','cell_to') + call results_closeGroup(results_addGroup('setup')) + call results_addAttribute('description','input data used to run the simulation','setup') call results_closeJobFile endif From e01d271ee4c668fc87c58cfcebf3e8f3d35676c1 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Tue, 27 Jul 2021 07:43:35 +0200 Subject: [PATCH 50/63] store geometry (for full reproducibility) --- src/grid/discretization_grid.f90 | 39 ++++++++++++++++++++------------ 1 file changed, 24 insertions(+), 15 deletions(-) diff --git a/src/grid/discretization_grid.f90 b/src/grid/discretization_grid.f90 index 049b3b60d..f94c9921b 100644 --- a/src/grid/discretization_grid.f90 +++ b/src/grid/discretization_grid.f90 @@ -68,11 +68,28 @@ subroutine discretization_grid_init(restart) devNull, z, z_offset integer, dimension(worldsize) :: & displs, sendcounts + integer :: fileUnit, myStat + integer(pI64) :: & + fileLength + character(len=:), allocatable :: & + fileContent + print'(/,a)', ' <<<+- discretization_grid init -+>>>'; flush(IO_STDOUT) + if(worldrank == 0) then - call readVTI(grid,geomSize,origin,materialAt_global) + inquire(file = trim(interface_geomFile), size=fileLength) + open(newunit=fileUnit, file=trim(interface_geomFile), access='stream',& + status='old', position='rewind', action='read',iostat=myStat) + if(myStat /= 0) call IO_error(100,ext_msg=trim(interface_geomFile)) + allocate(character(len=fileLength)::fileContent) + read(fileUnit) fileContent + close(fileUnit) + call readVTI(grid,geomSize,origin,materialAt_global,fileContent) + call results_openJobFile(parallel=.false.) + call results_writeDataset_str(fileContent,'setup',interface_geomFile,'geometry definition (grid solver)') + call results_closeJobFile else allocate(materialAt_global(0)) ! needed for IntelMPI endif @@ -157,7 +174,8 @@ end subroutine discretization_grid_init !> @brief Parse vtk image data (.vti) !> @details https://vtk.org/Wiki/VTK_XML_Formats !-------------------------------------------------------------------------------------------------- -subroutine readVTI(grid,geomSize,origin,material) +subroutine readVTI(grid,geomSize,origin,material, & + fileContent) integer, dimension(3), intent(out) :: & grid ! grid (across all processes!) @@ -166,28 +184,19 @@ subroutine readVTI(grid,geomSize,origin,material) origin ! origin (across all processes!) integer, dimension(:), intent(out), allocatable :: & material + character(len=*), intent(in) :: & + fileContent - character(len=:), allocatable :: fileContent, dataType, headerType + character(len=:), allocatable :: dataType, headerType logical :: inFile,inImage,gotCellData,compressed - integer :: fileUnit, myStat integer(pI64) :: & - fileLength, & !< length of the geom file (in characters) startPos, endPos, & s + grid = -1 geomSize = -1.0_pReal -!-------------------------------------------------------------------------------------------------- -! read raw data as stream - inquire(file = trim(interface_geomFile), size=fileLength) - open(newunit=fileUnit, file=trim(interface_geomFile), access='stream',& - status='old', position='rewind', action='read',iostat=myStat) - if(myStat /= 0) call IO_error(100,ext_msg=trim(interface_geomFile)) - allocate(character(len=fileLength)::fileContent) - read(fileUnit) fileContent - close(fileUnit) - inFile = .false. inImage = .false. gotCelldata = .false. From ddb0429a1d5451237e69202c0e312bc8e8c984d8 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Tue, 27 Jul 2021 07:57:04 +0200 Subject: [PATCH 51/63] store load case (full reproducibility for grid solver) --- src/IO.f90 | 3 ++- src/grid/DAMASK_grid.f90 | 8 ++++++++ src/grid/discretization_grid.f90 | 11 +---------- 3 files changed, 11 insertions(+), 11 deletions(-) diff --git a/src/IO.f90 b/src/IO.f90 index 399e2e6df..961e92d63 100644 --- a/src/IO.f90 +++ b/src/IO.f90 @@ -119,9 +119,10 @@ function IO_read(fileName) result(fileContent) character(len=:), allocatable :: fileContent integer :: & - fileLength, & fileUnit, & myStat + integer(pI64) :: & + fileLength inquire(file = fileName, size=fileLength) diff --git a/src/grid/DAMASK_grid.f90 b/src/grid/DAMASK_grid.f90 index 9c31b6f26..800003552 100644 --- a/src/grid/DAMASK_grid.f90 +++ b/src/grid/DAMASK_grid.f90 @@ -107,6 +107,8 @@ program DAMASK_grid step_bc, & step_mech, & step_discretization + character(len=:), allocatable :: & + fileContent !-------------------------------------------------------------------------------------------------- ! init DAMASK (all modules) @@ -127,6 +129,12 @@ program DAMASK_grid if (stagItMax < 0) call IO_error(301,ext_msg='maxStaggeredIter') if (maxCutBack < 0) call IO_error(301,ext_msg='maxCutBack') + if (worldrank == 0) then + fileContent = IO_read(interface_loadFile) + call results_openJobFile(parallel=.false.) + call results_writeDataset_str(fileContent,'setup',interface_loadFile,'load case definition (grid solver)') + call results_closeJobFile + endif config_load => YAML_parse_file(trim(interface_loadFile)) solver => config_load%get('solver') diff --git a/src/grid/discretization_grid.f90 b/src/grid/discretization_grid.f90 index f94c9921b..8592f7963 100644 --- a/src/grid/discretization_grid.f90 +++ b/src/grid/discretization_grid.f90 @@ -68,9 +68,6 @@ subroutine discretization_grid_init(restart) devNull, z, z_offset integer, dimension(worldsize) :: & displs, sendcounts - integer :: fileUnit, myStat - integer(pI64) :: & - fileLength character(len=:), allocatable :: & fileContent @@ -79,13 +76,7 @@ subroutine discretization_grid_init(restart) if(worldrank == 0) then - inquire(file = trim(interface_geomFile), size=fileLength) - open(newunit=fileUnit, file=trim(interface_geomFile), access='stream',& - status='old', position='rewind', action='read',iostat=myStat) - if(myStat /= 0) call IO_error(100,ext_msg=trim(interface_geomFile)) - allocate(character(len=fileLength)::fileContent) - read(fileUnit) fileContent - close(fileUnit) + fileContent = IO_read(interface_geomFile) call readVTI(grid,geomSize,origin,materialAt_global,fileContent) call results_openJobFile(parallel=.false.) call results_writeDataset_str(fileContent,'setup',interface_geomFile,'geometry definition (grid solver)') From 367c088b1661a817d6fe95aa39ddcf294b968804 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Tue, 27 Jul 2021 08:00:24 +0200 Subject: [PATCH 52/63] sorting according to dependency (no Makefile) --- src/commercialFEM_fileList.f90 | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/commercialFEM_fileList.f90 b/src/commercialFEM_fileList.f90 index 5fe754ce3..e1d53ca83 100644 --- a/src/commercialFEM_fileList.f90 +++ b/src/commercialFEM_fileList.f90 @@ -7,18 +7,18 @@ #include "IO.f90" #include "YAML_types.f90" #include "YAML_parse.f90" +#include "HDF5_utilities.f90" +#include "results.f90" #include "config.f90" #include "LAPACK_interface.f90" #include "math.f90" #include "rotations.f90" +#include "lattice.f90" #include "element.f90" -#include "HDF5_utilities.f90" -#include "results.f90" #include "geometry_plastic_nonlocal.f90" #include "discretization.f90" #include "Marc/discretization_Marc.f90" #include "material.f90" -#include "lattice.f90" #include "phase.f90" #include "phase_mechanical.f90" #include "phase_mechanical_elastic.f90" From 812b0f07f5fc8037df99cf62f26c9c7bccb9ed3e Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Tue, 27 Jul 2021 08:35:52 +0200 Subject: [PATCH 53/63] read file only once (per process) --- src/YAML_parse.f90 | 12 ++++++------ src/config.f90 | 19 ++++++++++++++++--- src/grid/DAMASK_grid.f90 | 4 ++-- 3 files changed, 24 insertions(+), 11 deletions(-) diff --git a/src/YAML_parse.f90 b/src/YAML_parse.f90 index d7cfa016f..2df09936d 100644 --- a/src/YAML_parse.f90 +++ b/src/YAML_parse.f90 @@ -14,7 +14,7 @@ module YAML_parse public :: & YAML_parse_init, & - YAML_parse_file + YAML_parse_str contains @@ -29,16 +29,16 @@ end subroutine YAML_parse_init !-------------------------------------------------------------------------------------------------- -!> @brief Parse a YAML file into a a structure of nodes. +!> @brief Parse a YAML string into a a structure of nodes. !-------------------------------------------------------------------------------------------------- -function YAML_parse_file(fname) result(node) +function YAML_parse_str(str) result(node) - character(len=*), intent(in) :: fname + character(len=*), intent(in) :: str class (tNode), pointer :: node - node => parse_flow(to_flow(IO_read(fname))) + node => parse_flow(to_flow(str)) -end function YAML_parse_file +end function YAML_parse_str !-------------------------------------------------------------------------------------------------- diff --git a/src/config.f90 b/src/config.f90 index 8ee25182f..b80658213 100644 --- a/src/config.f90 +++ b/src/config.f90 @@ -51,14 +51,17 @@ subroutine parse_material() inquire(file='material.yaml',exist=fileExists) if(.not. fileExists) call IO_error(100,ext_msg='material.yaml') + print*, 'reading material.yaml'; flush(IO_STDOUT) fileContent = IO_read('material.yaml') + if (worldrank == 0) then call results_openJobFile(parallel=.false.) call results_writeDataset_str(fileContent,'setup','material.yaml','DAMASK main configuration') call results_closeJobFile endif - config_material => YAML_parse_file('material.yaml') + + config_material => YAML_parse_str(fileContent) end subroutine parse_material @@ -73,16 +76,21 @@ subroutine parse_numerics() config_numerics => emptyDict + inquire(file='numerics.yaml', exist=fileExists) if (fileExists) then + print*, 'reading numerics.yaml'; flush(IO_STDOUT) fileContent = IO_read('numerics.yaml') + if (worldrank == 0) then call results_openJobFile(parallel=.false.) call results_writeDataset_str(fileContent,'setup','numerics.yaml','numerics configuration (optional)') call results_closeJobFile endif - config_numerics => YAML_parse_file('numerics.yaml') + + config_numerics => YAML_parse_str(fileContent) + endif end subroutine parse_numerics @@ -98,16 +106,21 @@ subroutine parse_debug() config_debug => emptyDict + inquire(file='debug.yaml', exist=fileExists) if (fileExists) then + print*, 'reading debug.yaml'; flush(IO_STDOUT) fileContent = IO_read('debug.yaml') + if (worldrank == 0) then call results_openJobFile(parallel=.false.) call results_writeDataset_str(fileContent,'setup','debug.yaml','debug configuration (optional)') call results_closeJobFile endif - config_debug => YAML_parse_file('debug.yaml') + + config_debug => YAML_parse_str(fileContent) + endif end subroutine parse_debug diff --git a/src/grid/DAMASK_grid.f90 b/src/grid/DAMASK_grid.f90 index 800003552..65c050b0c 100644 --- a/src/grid/DAMASK_grid.f90 +++ b/src/grid/DAMASK_grid.f90 @@ -129,13 +129,13 @@ program DAMASK_grid if (stagItMax < 0) call IO_error(301,ext_msg='maxStaggeredIter') if (maxCutBack < 0) call IO_error(301,ext_msg='maxCutBack') + fileContent = IO_read(interface_loadFile) if (worldrank == 0) then - fileContent = IO_read(interface_loadFile) call results_openJobFile(parallel=.false.) call results_writeDataset_str(fileContent,'setup',interface_loadFile,'load case definition (grid solver)') call results_closeJobFile endif - config_load => YAML_parse_file(trim(interface_loadFile)) + config_load => YAML_parse_str(fileContent) solver => config_load%get('solver') !-------------------------------------------------------------------------------------------------- From b9d4eb23ccebc8d6d4ecbfc02b1c8166743db7e3 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Tue, 27 Jul 2021 08:54:17 +0200 Subject: [PATCH 54/63] only rank 0 reads file for MPI --- src/config.f90 | 18 +++++++++--------- src/grid/DAMASK_grid.f90 | 4 +++- src/parallelization.f90 | 34 ++++++++++++++++++++++++++++++++-- 3 files changed, 44 insertions(+), 12 deletions(-) diff --git a/src/config.f90 b/src/config.f90 index b80658213..8f5680158 100644 --- a/src/config.f90 +++ b/src/config.f90 @@ -52,14 +52,14 @@ subroutine parse_material() inquire(file='material.yaml',exist=fileExists) if(.not. fileExists) call IO_error(100,ext_msg='material.yaml') - print*, 'reading material.yaml'; flush(IO_STDOUT) - fileContent = IO_read('material.yaml') - if (worldrank == 0) then + print*, 'reading material.yaml'; flush(IO_STDOUT) + fileContent = IO_read('material.yaml') call results_openJobFile(parallel=.false.) call results_writeDataset_str(fileContent,'setup','material.yaml','DAMASK main configuration') call results_closeJobFile endif + call parallelization_bcast_str(fileContent) config_material => YAML_parse_str(fileContent) @@ -80,14 +80,14 @@ subroutine parse_numerics() inquire(file='numerics.yaml', exist=fileExists) if (fileExists) then - print*, 'reading numerics.yaml'; flush(IO_STDOUT) - fileContent = IO_read('numerics.yaml') - if (worldrank == 0) then + print*, 'reading numerics.yaml'; flush(IO_STDOUT) + fileContent = IO_read('numerics.yaml') call results_openJobFile(parallel=.false.) call results_writeDataset_str(fileContent,'setup','numerics.yaml','numerics configuration (optional)') call results_closeJobFile endif + call parallelization_bcast_str(fileContent) config_numerics => YAML_parse_str(fileContent) @@ -110,14 +110,14 @@ subroutine parse_debug() inquire(file='debug.yaml', exist=fileExists) if (fileExists) then - print*, 'reading debug.yaml'; flush(IO_STDOUT) - fileContent = IO_read('debug.yaml') - if (worldrank == 0) then + print*, 'reading debug.yaml'; flush(IO_STDOUT) + fileContent = IO_read('debug.yaml') call results_openJobFile(parallel=.false.) call results_writeDataset_str(fileContent,'setup','debug.yaml','debug configuration (optional)') call results_closeJobFile endif + call parallelization_bcast_str(fileContent) config_debug => YAML_parse_str(fileContent) diff --git a/src/grid/DAMASK_grid.f90 b/src/grid/DAMASK_grid.f90 index 65c050b0c..89521a223 100644 --- a/src/grid/DAMASK_grid.f90 +++ b/src/grid/DAMASK_grid.f90 @@ -129,12 +129,14 @@ program DAMASK_grid if (stagItMax < 0) call IO_error(301,ext_msg='maxStaggeredIter') if (maxCutBack < 0) call IO_error(301,ext_msg='maxCutBack') - fileContent = IO_read(interface_loadFile) if (worldrank == 0) then + fileContent = IO_read(interface_loadFile) call results_openJobFile(parallel=.false.) call results_writeDataset_str(fileContent,'setup',interface_loadFile,'load case definition (grid solver)') call results_closeJobFile endif + + call parallelization_bcast_str(fileContent) config_load => YAML_parse_str(fileContent) solver => config_load%get('solver') diff --git a/src/parallelization.f90 b/src/parallelization.f90 index 534478cef..31da0f969 100644 --- a/src/parallelization.f90 +++ b/src/parallelization.f90 @@ -24,9 +24,18 @@ module parallelization worldrank = 0, & !< MPI worldrank (/=0 for MPI simulations only) worldsize = 1 !< MPI worldsize (/=1 for MPI simulations only) -#ifdef PETSC +#ifndef PETSC +public :: parallelization_broadcast_str + +contains +subroutine parallelization_bcast_str(string) + character(len=*), allocatable, intent(inout) :: string +end subroutine parallelization_bcast_str + +#else public :: & - parallelization_init + parallelization_init, & + parallelization_bcast_str contains @@ -101,6 +110,27 @@ subroutine parallelization_init !$ call omp_set_num_threads(OMP_NUM_THREADS) end subroutine parallelization_init + + +!-------------------------------------------------------------------------------------------------- +!> @brief Broadcast a string from process 0. +!-------------------------------------------------------------------------------------------------- +subroutine parallelization_bcast_str(string) + + character(len=:), allocatable, intent(inout) :: string + + integer :: strlen, ierr ! pI64 for strlen not supported by MPI + + + if (worldrank == 0) strlen = len(string,pI64) + call MPI_Bcast(strlen,1,MPI_INTEGER,0,MPI_COMM_WORLD, ierr) + if (worldrank /= 0) allocate(character(len=strlen)::string) + + call MPI_Bcast(string,strlen,MPI_CHARACTER,0,MPI_COMM_WORLD, ierr) + + +end subroutine parallelization_bcast_str + #endif end module parallelization From 26ae352a4c35c52edf99ea384f8deed2f808f85f Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Tue, 27 Jul 2021 09:31:57 +0200 Subject: [PATCH 55/63] make setup data easily accessible to the user --- PRIVATE | 2 +- python/damask/_result.py | 32 ++++++++++++++++++++++++++++++-- src/parallelization.f90 | 4 ++-- src/results.f90 | 2 +- 4 files changed, 34 insertions(+), 6 deletions(-) diff --git a/PRIVATE b/PRIVATE index 72c581038..7d48d8158 160000 --- a/PRIVATE +++ b/PRIVATE @@ -1 +1 @@ -Subproject commit 72c58103860e127d37ccf3a06827331de29406ca +Subproject commit 7d48d8158c1b2c0ab8bdd1c5d2baa3d020ae006d diff --git a/python/damask/_result.py b/python/damask/_result.py index 1fa376f63..0f69c4c5f 100644 --- a/python/damask/_result.py +++ b/python/damask/_result.py @@ -99,8 +99,10 @@ class Result: self.version_major = f.attrs['DADF5_version_major'] self.version_minor = f.attrs['DADF5_version_minor'] - if self.version_major != 0 or not 12 <= self.version_minor <= 13: + if self.version_major != 0 or not 12 <= self.version_minor <= 14: raise TypeError(f'Unsupported DADF5 version {self.version_major}.{self.version_minor}') + if self.version_major == 0 and self.version_minor < 14: + self.export_setup = None self.structured = 'cells' in f['geometry'].attrs.keys() @@ -1395,7 +1397,7 @@ class Result: def export_XDMF(self,output='*'): """ - Write XDMF file to directly visualize data in DADF5 file. + Write XDMF file to directly visualize data from DADF5 file. The XDMF format is only supported for structured grids with single phase and single constituent. @@ -1748,3 +1750,29 @@ class Result: if flatten: r = util.dict_flatten(r) return None if (type(r) == dict and r == {}) else r + + + def export_setup(self,output='*',overwrite=False): + """ + Export configuration files. + + Parameters + ---------- + output : (list of) str, optional + Names of the datasets to export to the file. + Defaults to '*', in which case all datasets are exported. + overwrite : boolean, optional + Overwrite existing configuration files. + Defaults to False. + + """ + with h5py.File(self.fname,'r') as f_in: + for out in _match(output,f_in['setup'].keys()): + description = f_in['/'.join(('setup',out))].attrs['description'] + if not h5py3: description = description.decode() + if not Path(out).exists() or overwrite: + with open(out,'w') as f_out: + f_out.write(f_in['/'.join(('setup',out))][()].decode()) + print(f"exported {description} to '{out}'") + else: + print(f"'{out}' exists, {description} not exported") diff --git a/src/parallelization.f90 b/src/parallelization.f90 index 31da0f969..a5e89561c 100644 --- a/src/parallelization.f90 +++ b/src/parallelization.f90 @@ -25,11 +25,11 @@ module parallelization worldsize = 1 !< MPI worldsize (/=1 for MPI simulations only) #ifndef PETSC -public :: parallelization_broadcast_str +public :: parallelization_bcast_str contains subroutine parallelization_bcast_str(string) - character(len=*), allocatable, intent(inout) :: string + character(len=:), allocatable, intent(inout) :: string end subroutine parallelization_bcast_str #else diff --git a/src/results.f90 b/src/results.f90 index 8b04800ef..9c4507938 100644 --- a/src/results.f90 +++ b/src/results.f90 @@ -74,7 +74,7 @@ subroutine results_init(restart) if(.not. restart) then resultsFile = HDF5_openFile(getSolverJobName()//'.hdf5','w') call results_addAttribute('DADF5_version_major',0) - call results_addAttribute('DADF5_version_minor',13) + call results_addAttribute('DADF5_version_minor',14) call get_command_argument(0,commandLine) call results_addAttribute('creator',trim(commandLine)//' '//DAMASKVERSION) call results_addAttribute('created',now()) From 64d52dbbf47632891f08b61a686b18be2de7fc0c Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Tue, 27 Jul 2021 09:46:51 +0200 Subject: [PATCH 56/63] store input deck, MSC.Marc now 100% reproducible --- src/Marc/discretization_Marc.f90 | 8 +++++++- src/parallelization.f90 | 2 +- 2 files changed, 8 insertions(+), 2 deletions(-) diff --git a/src/Marc/discretization_Marc.f90 b/src/Marc/discretization_Marc.f90 index d93eea2d5..fd6b8699d 100644 --- a/src/Marc/discretization_Marc.f90 +++ b/src/Marc/discretization_Marc.f90 @@ -216,7 +216,13 @@ subroutine inputRead(elem,node0_elem,connectivity_elem,materialAt) mapElemSet !< list of elements in elementSet - inputFile = IO_readlines(trim(getSolverJobName())//trim(InputFileExtension)) + call results_openJobFile + call results_writeDataset_str(IO_read(trim(getSolverJobName())//InputFileExtension), 'setup', & + trim(getSolverJobName())//InputFileExtension, & + 'MSC.Marc input deck') + call results_closeJobFile + + inputFile = IO_readlines(trim(getSolverJobName())//InputFileExtension) call inputRead_fileFormat(fileFormatVersion, & inputFile) call inputRead_tableStyles(initialcondTableStyle,hypoelasticTableStyle, & diff --git a/src/parallelization.f90 b/src/parallelization.f90 index a5e89561c..0f8b1b8f8 100644 --- a/src/parallelization.f90 +++ b/src/parallelization.f90 @@ -122,7 +122,7 @@ subroutine parallelization_bcast_str(string) integer :: strlen, ierr ! pI64 for strlen not supported by MPI - if (worldrank == 0) strlen = len(string,pI64) + if (worldrank == 0) strlen = string call MPI_Bcast(strlen,1,MPI_INTEGER,0,MPI_COMM_WORLD, ierr) if (worldrank /= 0) allocate(character(len=strlen)::string) From a6f7e4f1a6cc2759450d2ef7fccd84c482fe072c Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Tue, 27 Jul 2021 09:50:23 +0200 Subject: [PATCH 57/63] play it safe --- src/HDF5_utilities.f90 | 31 +++++++++++++++++++------------ src/parallelization.f90 | 2 +- 2 files changed, 20 insertions(+), 13 deletions(-) diff --git a/src/HDF5_utilities.f90 b/src/HDF5_utilities.f90 index bd88feee2..2d2632e13 100644 --- a/src/HDF5_utilities.f90 +++ b/src/HDF5_utilities.f90 @@ -1479,27 +1479,34 @@ subroutine HDF5_write_str(dataset,loc_id,datasetName) integer(HID_T), intent(in) :: loc_id character(len=*), intent(in) :: datasetName !< name of the dataset in the file - INTEGER(HID_T) :: filetype_id, space_id, dataset_id - INTEGER :: hdferr - + integer(HID_T) :: filetype_id, space_id, dataset_id + integer :: hdferr character(len=len_trim(dataset)+1,kind=C_CHAR), dimension(1), target :: dataset_ type(C_PTR), target, dimension(1) :: ptr - + dataset_(1) = trim(dataset)//C_NULL_CHAR ptr(1) = c_loc(dataset_(1)) - - call h5tcopy_f(H5T_STRING, filetype_id, hdferr) - call h5tset_size_f(filetype_id, int(len(dataset_),HSIZE_T), hdferr) - - call h5screate_f(H5S_SCALAR_F, space_id, hdferr) - call h5dcreate_f(loc_id, datasetName, H5T_STRING, space_id, dataset_id, hdferr) - call h5dwrite_f(dataset_id, H5T_STRING, c_loc(ptr), hdferr); - + call h5tcopy_f(H5T_STRING, filetype_id, hdferr) + if(hdferr < 0) error stop 'HDF5 error' + call h5tset_size_f(filetype_id, int(len(dataset_),HSIZE_T), hdferr) + if(hdferr < 0) error stop 'HDF5 error' + + call h5screate_f(H5S_SCALAR_F, space_id, hdferr) + if(hdferr < 0) error stop 'HDF5 error' + call h5dcreate_f(loc_id, datasetName, H5T_STRING, space_id, dataset_id, hdferr) + if(hdferr < 0) error stop 'HDF5 error' + + call h5dwrite_f(dataset_id, H5T_STRING, c_loc(ptr), hdferr) + if(hdferr < 0) error stop 'HDF5 error' + call h5dclose_f(dataset_id, hdferr) + if(hdferr < 0) error stop 'HDF5 error' call h5sclose_f(space_id, hdferr) + if(hdferr < 0) error stop 'HDF5 error' call h5tclose_f(filetype_id, hdferr) + if(hdferr < 0) error stop 'HDF5 error' end subroutine HDF5_write_str diff --git a/src/parallelization.f90 b/src/parallelization.f90 index 0f8b1b8f8..8ea22ef0a 100644 --- a/src/parallelization.f90 +++ b/src/parallelization.f90 @@ -122,7 +122,7 @@ subroutine parallelization_bcast_str(string) integer :: strlen, ierr ! pI64 for strlen not supported by MPI - if (worldrank == 0) strlen = string + if (worldrank == 0) strlen = len(string) call MPI_Bcast(strlen,1,MPI_INTEGER,0,MPI_COMM_WORLD, ierr) if (worldrank /= 0) allocate(character(len=strlen)::string) From af714cfd5b2a2e1d12c28cb0c3b83bce108f4469 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Tue, 27 Jul 2021 13:56:07 +0200 Subject: [PATCH 58/63] broken Intel compiler (do concurrent) errors occur under unclear conditions, better avoid at the moment: https://community.intel.com/t5/Intel-Fortran-Compiler/do-concurrent-broken-with-openmp/m-p/1301528#M156942 https://community.intel.com/t5/Intel-Fortran-Compiler/Bug-with-do-concurrent-and-openmp/m-p/1173473m --- src/element.f90 | 1 + src/math.f90 | 51 +++++++++++++++++++++++++++++++++++++++++++++---- 2 files changed, 48 insertions(+), 4 deletions(-) diff --git a/src/element.f90 b/src/element.f90 index 5b7af36bd..b061f942d 100644 --- a/src/element.f90 +++ b/src/element.f90 @@ -151,6 +151,7 @@ module element integer, dimension(NIPNEIGHBOR(CELLTYPE(1)),NIP(1)), parameter :: IPNEIGHBOR1 = & reshape([& -2,-3,-1 & +! Note: This fix is for gfortran 9 only. gfortran 8 supports neither, gfortran > 9 both variants #if !defined(__GFORTRAN__) ],shape(IPNEIGHBOR1)) #else diff --git a/src/math.f90 b/src/math.f90 index cbc924f08..71cbcffec 100644 --- a/src/math.f90 +++ b/src/math.f90 @@ -90,12 +90,12 @@ subroutine math_init integer, dimension(:), allocatable :: randInit class(tNode), pointer :: & num_generic - + print'(/,a)', ' <<<+- math init -+>>>'; flush(IO_STDOUT) num_generic => config_numerics%get('generic',defaultVal=emptyDict) randomSeed = num_generic%get_asInt('random_seed', defaultVal = 0) - + call random_seed(size=randSize) allocate(randInit(randSize)) if (randomSeed > 0) then @@ -260,9 +260,15 @@ pure function math_identity4th() integer :: i,j,k,l +#ifndef __INTEL_COMPILER do concurrent(i=1:3, j=1:3, k=1:3, l=1:3) math_identity4th(i,j,k,l) = 0.5_pReal*(math_I3(i,k)*math_I3(j,l)+math_I3(i,l)*math_I3(j,k)) enddo +#else + do i=1,3; do j=1,3; do k=1,3; do l=1,3 + math_identity4th(i,j,k,l) = 0.5_pReal*(math_I3(i,k)*math_I3(j,l)+math_I3(i,l)*math_I3(j,k)) + enddo; enddo; enddo; enddo +#endif end function math_identity4th @@ -330,9 +336,15 @@ pure function math_outer(A,B) integer :: i,j +#ifndef __INTEL_COMPILER do concurrent(i=1:size(A,1), j=1:size(B,1)) math_outer(i,j) = A(i)*B(j) enddo +#else + do i=1,size(A,1); do j=1,size(B,1) + math_outer(i,j) = A(i)*B(j) + enddo; enddo +#endif end function math_outer @@ -373,9 +385,15 @@ pure function math_mul3333xx33(A,B) integer :: i,j +#ifndef __INTEL_COMPILER do concurrent(i=1:3, j=1:3) math_mul3333xx33(i,j) = sum(A(i,j,1:3,1:3)*B(1:3,1:3)) enddo +#else + do i=1,3; do j=1,3 + math_mul3333xx33(i,j) = sum(A(i,j,1:3,1:3)*B(1:3,1:3)) + enddo; enddo +#endif end function math_mul3333xx33 @@ -390,9 +408,16 @@ pure function math_mul3333xx3333(A,B) real(pReal), dimension(3,3,3,3), intent(in) :: B real(pReal), dimension(3,3,3,3) :: math_mul3333xx3333 + +#ifndef __INTEL_COMPILER do concurrent(i=1:3, j=1:3, k=1:3, l=1:3) math_mul3333xx3333(i,j,k,l) = sum(A(i,j,1:3,1:3)*B(1:3,1:3,k,l)) enddo +#else + do i=1,3; do j=1,3; do k=1,3; do l=1,3 + math_mul3333xx3333(i,j,k,l) = sum(A(i,j,1:3,1:3)*B(1:3,1:3,k,l)) + enddo; enddo; enddo; enddo +#endif end function math_mul3333xx3333 @@ -725,9 +750,15 @@ pure function math_3333to99(m3333) integer :: i,j +#ifndef __INTEL_COMPILER do concurrent(i=1:9, j=1:9) math_3333to99(i,j) = m3333(MAPPLAIN(1,i),MAPPLAIN(2,i),MAPPLAIN(1,j),MAPPLAIN(2,j)) enddo +#else + do i=1,9; do j=1,9 + math_3333to99(i,j) = m3333(MAPPLAIN(1,i),MAPPLAIN(2,i),MAPPLAIN(1,j),MAPPLAIN(2,j)) + enddo; enddo +#endif end function math_3333to99 @@ -742,10 +773,15 @@ pure function math_99to3333(m99) integer :: i,j - +#ifndef __INTEL_COMPILER do concurrent(i=1:9, j=1:9) math_99to3333(MAPPLAIN(1,i),MAPPLAIN(2,i),MAPPLAIN(1,j),MAPPLAIN(2,j)) = m99(i,j) enddo +#else + do i=1,9; do j=1,9 + math_99to3333(MAPPLAIN(1,i),MAPPLAIN(2,i),MAPPLAIN(1,j),MAPPLAIN(2,j)) = m99(i,j) + enddo; enddo +#endif end function math_99to3333 @@ -772,9 +808,15 @@ pure function math_sym3333to66(m3333,weighted) w = NRMMANDEL endif +#ifndef __INTEL_COMPILER do concurrent(i=1:6, j=1:6) math_sym3333to66(i,j) = w(i)*w(j)*m3333(MAPNYE(1,i),MAPNYE(2,i),MAPNYE(1,j),MAPNYE(2,j)) enddo +#else + do i=1,6; do j=1,6 + math_sym3333to66(i,j) = w(i)*w(j)*m3333(MAPNYE(1,i),MAPNYE(2,i),MAPNYE(1,j),MAPNYE(2,j)) + enddo; enddo +#endif end function math_sym3333to66 @@ -874,10 +916,11 @@ subroutine math_eigh(w,v,error,m) real(pReal), dimension(size(m,1)), intent(out) :: w !< eigenvalues real(pReal), dimension(size(m,1),size(m,1)), intent(out) :: v !< eigenvectors logical, intent(out) :: error - + integer :: ierr real(pReal), dimension(size(m,1)**2) :: work + v = m ! copy matrix to input (doubles as output) array call dsyev('V','U',size(m,1),v,size(m,1),w,work,size(work,1),ierr) error = (ierr /= 0) From 5d0c0f42831fd740c8dcf34651d738a2d8d24c02 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Tue, 27 Jul 2021 16:44:52 +0200 Subject: [PATCH 59/63] outdated instructions --- src/DAMASK_interface.f90 | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/DAMASK_interface.f90 b/src/DAMASK_interface.f90 index 52f58d363..ca3179afc 100644 --- a/src/DAMASK_interface.f90 +++ b/src/DAMASK_interface.f90 @@ -139,10 +139,10 @@ subroutine DAMASK_interface_init print'(a)', ' Optional arguments:' print'(/,a)',' --workingdirectory PathToWorkingDirectory' print'(a)', ' Specifies the working directory and overwrites the default ./' - print'(a)', ' Make sure the file "material.config" exists in the working' + print'(a)', ' Make sure the file "material.yaml" exists in the working' print'(a)', ' directory.' - print'(a)', ' For further configuration place "numerics.config"' - print'(a)',' and "debug.config" in that directory.' + print'(a)', ' For further configuration place "numerics.yaml"' + print'(a)',' and "debug.yaml" in that directory.' print'(/,a)',' --restart N' print'(a)', ' Reads in increment N and continues with calculating' print'(a)', ' increment N+1 based on this.' From 29155e325f9bc9c8a31d5c338438f5e8424cea5f Mon Sep 17 00:00:00 2001 From: Test User Date: Thu, 29 Jul 2021 12:58:11 +0200 Subject: [PATCH 60/63] [skip ci] updated version information after successful test of v3.0.0-alpha4-205-g5d0c0f428 --- python/damask/VERSION | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/python/damask/VERSION b/python/damask/VERSION index 13d57229b..2a51b3669 100644 --- a/python/damask/VERSION +++ b/python/damask/VERSION @@ -1 +1 @@ -v3.0.0-alpha4-182-gac6d31b1f +v3.0.0-alpha4-205-g5d0c0f428 From 044a0489446f2b18051c528e7cfac304d03ff7d7 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Tue, 27 Jul 2021 11:59:46 +0200 Subject: [PATCH 61/63] taking care of corner cases (e.g. restart) adjusting tests to take care of new 'setup' group --- PRIVATE | 2 +- src/HDF5_utilities.f90 | 4 ++++ src/IO.f90 | 8 ++++---- src/config.f90 | 2 +- src/grid/DAMASK_grid.f90 | 6 ++++-- src/grid/discretization_grid.f90 | 6 ++++-- src/math.f90 | 2 +- src/prec.f90 | 1 + src/results.f90 | 19 +++++++++++++++++-- 9 files changed, 37 insertions(+), 13 deletions(-) diff --git a/PRIVATE b/PRIVATE index 7d48d8158..9699f20f2 160000 --- a/PRIVATE +++ b/PRIVATE @@ -1 +1 @@ -Subproject commit 7d48d8158c1b2c0ab8bdd1c5d2baa3d020ae006d +Subproject commit 9699f20f21f8a5f532c735a1aa9daeba395da94d diff --git a/src/HDF5_utilities.f90 b/src/HDF5_utilities.f90 index 2d2632e13..fa62e4840 100644 --- a/src/HDF5_utilities.f90 +++ b/src/HDF5_utilities.f90 @@ -1473,6 +1473,10 @@ subroutine HDF5_write_real7(dataset,loc_id,datasetName,parallel) end subroutine HDF5_write_real7 +!-------------------------------------------------------------------------------------------------- +!> @brief Write dataset of type string (scalar). +!> @details Not collective, must be called by one process at at time. +!-------------------------------------------------------------------------------------------------- subroutine HDF5_write_str(dataset,loc_id,datasetName) character(len=*), intent(in) :: dataset diff --git a/src/IO.f90 b/src/IO.f90 index 961e92d63..a32bd5ef0 100644 --- a/src/IO.f90 +++ b/src/IO.f90 @@ -128,19 +128,19 @@ function IO_read(fileName) result(fileContent) inquire(file = fileName, size=fileLength) open(newunit=fileUnit, file=fileName, access='stream',& status='old', position='rewind', action='read',iostat=myStat) - if(myStat /= 0) call IO_error(100,ext_msg=trim(fileName)) + if (myStat /= 0) call IO_error(100,ext_msg=trim(fileName)) allocate(character(len=fileLength)::fileContent) - if(fileLength==0) then + if (fileLength==0) then close(fileUnit) return endif read(fileUnit,iostat=myStat) fileContent - if(myStat /= 0) call IO_error(102,ext_msg=trim(fileName)) + if (myStat /= 0) call IO_error(102,ext_msg=trim(fileName)) close(fileUnit) if (scan(fileContent(:index(fileContent,LF)),CR//LF) /= 0) fileContent = CRLF2LF(fileContent) - if(fileContent(fileLength:fileLength) /= IO_EOL) fileContent = fileContent//IO_EOL ! ensure EOL@EOF + if (fileContent(fileLength:fileLength) /= IO_EOL) fileContent = fileContent//IO_EOL ! ensure EOL@EOF end function IO_read diff --git a/src/config.f90 b/src/config.f90 index 8f5680158..c460ce7c5 100644 --- a/src/config.f90 +++ b/src/config.f90 @@ -56,7 +56,7 @@ subroutine parse_material() print*, 'reading material.yaml'; flush(IO_STDOUT) fileContent = IO_read('material.yaml') call results_openJobFile(parallel=.false.) - call results_writeDataset_str(fileContent,'setup','material.yaml','DAMASK main configuration') + call results_writeDataset_str(fileContent,'setup','material.yaml','main configuration') call results_closeJobFile endif call parallelization_bcast_str(fileContent) diff --git a/src/grid/DAMASK_grid.f90 b/src/grid/DAMASK_grid.f90 index 89521a223..b6e45a75b 100644 --- a/src/grid/DAMASK_grid.f90 +++ b/src/grid/DAMASK_grid.f90 @@ -108,7 +108,7 @@ program DAMASK_grid step_mech, & step_discretization character(len=:), allocatable :: & - fileContent + fileContent, fname !-------------------------------------------------------------------------------------------------- ! init DAMASK (all modules) @@ -131,8 +131,10 @@ program DAMASK_grid if (worldrank == 0) then fileContent = IO_read(interface_loadFile) + fname = interface_loadFile + if (scan(fname,'/') /= 0) fname = fname(scan(fname,'/',.true.)+1:) call results_openJobFile(parallel=.false.) - call results_writeDataset_str(fileContent,'setup',interface_loadFile,'load case definition (grid solver)') + call results_writeDataset_str(fileContent,'setup',fname,'load case definition (grid solver)') call results_closeJobFile endif diff --git a/src/grid/discretization_grid.f90 b/src/grid/discretization_grid.f90 index 8592f7963..1c88e5817 100644 --- a/src/grid/discretization_grid.f90 +++ b/src/grid/discretization_grid.f90 @@ -69,7 +69,7 @@ subroutine discretization_grid_init(restart) integer, dimension(worldsize) :: & displs, sendcounts character(len=:), allocatable :: & - fileContent + fileContent, fname print'(/,a)', ' <<<+- discretization_grid init -+>>>'; flush(IO_STDOUT) @@ -78,8 +78,10 @@ subroutine discretization_grid_init(restart) if(worldrank == 0) then fileContent = IO_read(interface_geomFile) call readVTI(grid,geomSize,origin,materialAt_global,fileContent) + fname = interface_geomFile + if (scan(fname,'/') /= 0) fname = fname(scan(fname,'/',.true.)+1:) call results_openJobFile(parallel=.false.) - call results_writeDataset_str(fileContent,'setup',interface_geomFile,'geometry definition (grid solver)') + call results_writeDataset_str(fileContent,'setup',fname,'geometry definition (grid solver)') call results_closeJobFile else allocate(materialAt_global(0)) ! needed for IntelMPI diff --git a/src/math.f90 b/src/math.f90 index a48f44b93..052ff7fa0 100644 --- a/src/math.f90 +++ b/src/math.f90 @@ -1209,7 +1209,7 @@ subroutine selfTest error stop 'math_sym33to6/math_6toSym33' call random_number(t66) - if(any(dNeq(math_sym3333to66(math_66toSym3333(t66)),t66))) & + if(any(dNeq(math_sym3333to66(math_66toSym3333(t66)),t66,1.0e-15_pReal))) & error stop 'math_sym3333to66/math_66toSym3333' call random_number(v6) diff --git a/src/prec.f90 b/src/prec.f90 index 3090eb67c..4b70ffbfa 100644 --- a/src/prec.f90 +++ b/src/prec.f90 @@ -89,6 +89,7 @@ end subroutine prec_init ! replaces "==" but for certain (relative) tolerance. Counterpart to dNeq ! https://randomascii.wordpress.com/2012/02/25/comparing-floating-point-numbers-2012-edition/ ! AlmostEqualRelative +! ToDo: Use 'spacing': https://gcc.gnu.org/onlinedocs/gfortran/SPACING.html#SPACING !-------------------------------------------------------------------------------------------------- logical elemental pure function dEq(a,b,tol) diff --git a/src/results.f90 b/src/results.f90 index 9c4507938..5553a883a 100644 --- a/src/results.f90 +++ b/src/results.f90 @@ -65,13 +65,17 @@ subroutine results_init(restart) logical, intent(in) :: restart character(len=pPathLen) :: commandLine + integer :: hdferr + integer(HID_T) :: group_id + character(len=:), allocatable :: date + print'(/,a)', ' <<<+- results init -+>>>'; flush(IO_STDOUT) print*, 'M. Diehl et al., Integrating Materials and Manufacturing Innovation 6(1):83–91, 2017' print*, 'https://doi.org/10.1007/s40192-017-0084-5'//IO_EOL - if(.not. restart) then + if (.not. restart) then resultsFile = HDF5_openFile(getSolverJobName()//'.hdf5','w') call results_addAttribute('DADF5_version_major',0) call results_addAttribute('DADF5_version_minor',14) @@ -84,9 +88,20 @@ subroutine results_init(restart) call results_addAttribute('description','mappings to place data in space','cell_to') call results_closeGroup(results_addGroup('setup')) call results_addAttribute('description','input data used to run the simulation','setup') - call results_closeJobFile + else + date = now() + call results_openJobFile + call get_command(commandLine) + call results_addAttribute('call (restart at '//date//')',trim(commandLine)) + call h5gmove_f(resultsFile,'setup','tmp',hdferr) + call results_addAttribute('description','input data used to run the simulation (backup from restart at '//date//')','tmp') + call results_closeGroup(results_addGroup('setup')) + call results_addAttribute('description','input data used to run the simulation','setup') + call h5gmove_f(resultsFile,'tmp','setup/backup',hdferr) endif + call results_closeJobFile + end subroutine results_init From e8312a49ed3a509e8995e808e590c676344d6287 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Mon, 2 Aug 2021 10:38:59 +0200 Subject: [PATCH 62/63] polishing --- python/damask/_result.py | 21 ++++++++++++--------- src/config.f90 | 4 ++-- src/results.f90 | 4 ++-- 3 files changed, 16 insertions(+), 13 deletions(-) diff --git a/python/damask/_result.py b/python/damask/_result.py index 0f69c4c5f..d502969a1 100644 --- a/python/damask/_result.py +++ b/python/damask/_result.py @@ -1766,13 +1766,16 @@ class Result: Defaults to False. """ - with h5py.File(self.fname,'r') as f_in: - for out in _match(output,f_in['setup'].keys()): - description = f_in['/'.join(('setup',out))].attrs['description'] - if not h5py3: description = description.decode() - if not Path(out).exists() or overwrite: - with open(out,'w') as f_out: - f_out.write(f_in['/'.join(('setup',out))][()].decode()) - print(f"exported {description} to '{out}'") + def export(name,obj,output,overwrite): + if type(obj) == h5py.Dataset and _match(output,[name]): + d = obj.attrs['description'] if h5py3 else obj.attrs['description'].decode() + if not Path(name).exists() or overwrite: + with open(name,'w') as f_out: f_out.write(obj[()].decode()) + print(f"Exported {d} to '{name}'.") else: - print(f"'{out}' exists, {description} not exported") + print(f"'{name}' exists, {d} not exported.") + elif type(obj) == h5py.Group: + os.makedirs(name, exist_ok=True) + + with h5py.File(self.fname,'r') as f_in: + f_in['setup'].visititems(partial(export,output=output,overwrite=overwrite)) diff --git a/src/config.f90 b/src/config.f90 index c460ce7c5..ecde0831c 100644 --- a/src/config.f90 +++ b/src/config.f90 @@ -84,7 +84,7 @@ subroutine parse_numerics() print*, 'reading numerics.yaml'; flush(IO_STDOUT) fileContent = IO_read('numerics.yaml') call results_openJobFile(parallel=.false.) - call results_writeDataset_str(fileContent,'setup','numerics.yaml','numerics configuration (optional)') + call results_writeDataset_str(fileContent,'setup','numerics.yaml','numerics configuration') call results_closeJobFile endif call parallelization_bcast_str(fileContent) @@ -114,7 +114,7 @@ subroutine parse_debug() print*, 'reading debug.yaml'; flush(IO_STDOUT) fileContent = IO_read('debug.yaml') call results_openJobFile(parallel=.false.) - call results_writeDataset_str(fileContent,'setup','debug.yaml','debug configuration (optional)') + call results_writeDataset_str(fileContent,'setup','debug.yaml','debug configuration') call results_closeJobFile endif call parallelization_bcast_str(fileContent) diff --git a/src/results.f90 b/src/results.f90 index 5553a883a..f6bc4a045 100644 --- a/src/results.f90 +++ b/src/results.f90 @@ -94,10 +94,10 @@ subroutine results_init(restart) call get_command(commandLine) call results_addAttribute('call (restart at '//date//')',trim(commandLine)) call h5gmove_f(resultsFile,'setup','tmp',hdferr) - call results_addAttribute('description','input data used to run the simulation (backup from restart at '//date//')','tmp') + call results_addAttribute('description','input data used to run the simulation up to restart at '//date,'tmp') call results_closeGroup(results_addGroup('setup')) call results_addAttribute('description','input data used to run the simulation','setup') - call h5gmove_f(resultsFile,'tmp','setup/backup',hdferr) + call h5gmove_f(resultsFile,'tmp','setup/previous',hdferr) endif call results_closeJobFile From d8c03ad0d55d4cabeb6199a72e51c6852bef2eaa Mon Sep 17 00:00:00 2001 From: Test User Date: Tue, 3 Aug 2021 18:19:37 +0200 Subject: [PATCH 63/63] [skip ci] updated version information after successful test of v3.0.0-alpha4-221-g4a8c83611 --- python/damask/VERSION | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/python/damask/VERSION b/python/damask/VERSION index 2a51b3669..b7d94c61e 100644 --- a/python/damask/VERSION +++ b/python/damask/VERSION @@ -1 +1 @@ -v3.0.0-alpha4-205-g5d0c0f428 +v3.0.0-alpha4-221-g4a8c83611