Merge branch 'development' into vectorize_rotation

PRIVATE

@@ -1 +1 @@
-Subproject commit 038af521a1ef70ed77b132c426bc1a4880db01ef
+Subproject commit 72d526e5750366a9efe4d1fd9d92e0d1ecd2cd38

cmake/Compiler-GNU.cmake

@@ -14,7 +14,7 @@ elseif (OPTIMIZATION STREQUAL "AGGRESSIVE")
   set (OPTIMIZATION_FLAGS "-O3 -ffast-math -funroll-loops -ftree-vectorize")
 endif ()
 
-set (STANDARD_CHECK "-std=f2008ts -pedantic-errors" )
+set (STANDARD_CHECK "-std=f2018 -pedantic-errors" )
 set (LINKER_FLAGS  "${LINKER_FLAGS} -Wl")
 # options parsed directly to the linker
 set (LINKER_FLAGS  "${LINKER_FLAGS},-undefined,dynamic_lookup" )

python/damask/grid_filters.py

@@ -361,7 +361,7 @@ def node_2_cell(node_data):
          + _np.roll(node_data,1,(0,))  + _np.roll(node_data,1,(1,))  + _np.roll(node_data,1,(2,))
          + _np.roll(node_data,1,(0,1)) + _np.roll(node_data,1,(1,2)) + _np.roll(node_data,1,(2,0)))*0.125
 
-    return c[:-1,:-1,:-1]
+    return c[1:,1:,1:]
 
 
 def node_coord0_gridSizeOrigin(coord0,ordered=True):

python/tests/test_Table.py

@@ -18,7 +18,7 @@ def reference_dir(reference_dir_base):
     return os.path.join(reference_dir_base,'Table')
 
 class TestTable:
-    
+
     def test_get_scalar(self,default):
         d = default.get('s')
         assert np.allclose(d,1.0) and d.shape[1:] == (1,)
@@ -29,12 +29,12 @@ class TestTable:
 
     def test_get_tensor(self,default):
         d = default.get('F')
-        assert np.allclose(d,1.0) and d.shape[1:] == (3,3) 
+        assert np.allclose(d,1.0) and d.shape[1:] == (3,3)
 
     def test_get_component(self,default):
         d = default.get('5_F')
         assert np.allclose(d,1.0) and d.shape[1:] == (1,)
-  
+
     def test_write_read_str(self,default,tmpdir):
         default.to_ASCII(str(tmpdir.join('default.txt')))
         new = Table.from_ASCII(str(tmpdir.join('default.txt')))
@@ -69,15 +69,20 @@ class TestTable:
     def test_read_strange(self,reference_dir,fname):
         with open(os.path.join(reference_dir,fname)) as f:
             Table.from_ASCII(f)
-    
+
     def test_set(self,default):
         default.set('F',np.zeros((5,3,3)),'set to zero')
         d=default.get('F')
         assert np.allclose(d,0.0) and d.shape[1:] == (3,3)
 
+    def test_set_component(self,default):
+        default.set('1_F',np.zeros((5)),'set to zero')
+        d=default.get('F')
+        assert np.allclose(d[...,0,0],0.0) and d.shape[1:] == (3,3)
+
     def test_labels(self,default):
         assert default.labels == ['F','v','s']
-        
+
     def test_add(self,default):
         d = np.random.random((5,9))
         default.add('nine',d,'random data')

python/tests/test_grid_filters.py

@@ -42,12 +42,25 @@ class TestGridFilters:
          assert np.allclose(grid_filters.node_displacement_fluct(size,F),
                             grid_filters.cell_2_node(grid_filters.cell_displacement_fluct(size,F)))
 
-    def test_interpolation_nonperiodic(self):
+    def test_interpolation_to_node(self):
          size = np.random.random(3)
          grid = np.random.randint(8,32,(3))
          F    = np.random.random(tuple(grid)+(3,3))
-         assert np.allclose(grid_filters.node_coord(size,F) [1:-1,1:-1,1:-1],grid_filters.cell_2_node(
+         assert np.allclose(grid_filters.node_coord(size,F) [1:-1,1:-1,1:-1],
-                            grid_filters.cell_coord(size,F))[1:-1,1:-1,1:-1])
+                            grid_filters.cell_2_node(grid_filters.cell_coord(size,F))[1:-1,1:-1,1:-1])
+
+    def test_interpolation_to_cell(self):
+         grid = np.random.randint(1,30,(3))
+
+         node_coord_x = np.linspace(0,np.pi*2,num=grid[0]+1)
+         node_field_x = np.cos(node_coord_x)
+         node_field   = np.broadcast_to(node_field_x.reshape(-1,1,1),grid+1)
+
+         cell_coord_x = node_coord_x[:-1]+node_coord_x[1]*.5
+         cell_field_x = np.interp(cell_coord_x,node_coord_x,node_field_x,period=np.pi*2.)
+         cell_field   = np.broadcast_to(cell_field_x.reshape(-1,1,1),grid)
+
+         assert np.allclose(cell_field,grid_filters.node_2_cell(node_field))
 
     @pytest.mark.parametrize('mode',['cell','node'])
     def test_coord0_origin(self,mode):

src/DAMASK_interface.f90

@@ -106,7 +106,7 @@ subroutine DAMASK_interface_init
     typeSize
   integer, dimension(8) :: &
     dateAndTime
-  integer        :: mpi_err
+  integer        :: err
   PetscErrorCode :: petsc_err
   external :: &
     quit
@@ -118,8 +118,8 @@ subroutine DAMASK_interface_init
 #ifdef _OPENMP
   ! If openMP is enabled, check if the MPI libary supports it and initialize accordingly.
   ! Otherwise, the first call to PETSc will do the initialization.
-  call MPI_Init_Thread(MPI_THREAD_FUNNELED,threadLevel,mpi_err)
+  call MPI_Init_Thread(MPI_THREAD_FUNNELED,threadLevel,err)
-  if (mpi_err /= 0) call quit(1)
+  if (err /= 0) call quit(1)
   if (threadLevel<MPI_THREAD_FUNNELED) then
     write(6,'(/,a)') ' ERROR: MPI library does not support OpenMP'
     call quit(1)
@@ -128,10 +128,10 @@ subroutine DAMASK_interface_init
   call PETScInitializeNoArguments(petsc_err)                                                        ! according to PETSc manual, that should be the first line in the code
   CHKERRQ(petsc_err)                                                                                ! this is a macro definition, it is case sensitive
 
-  call MPI_Comm_rank(PETSC_COMM_WORLD,worldrank,mpi_err)
+  call MPI_Comm_rank(PETSC_COMM_WORLD,worldrank,err)
-  if (mpi_err /= 0) call quit(1)
+  if (err /= 0) call quit(1)
-  call MPI_Comm_size(PETSC_COMM_WORLD,worldsize,mpi_err)
+  call MPI_Comm_size(PETSC_COMM_WORLD,worldsize,err)
-  if (mpi_err /= 0) call quit(1)
+  if (err /= 0) call quit(1)
 
   mainProcess: if (worldrank == 0) then
     if (output_unit /= 6) then
@@ -181,22 +181,23 @@ subroutine DAMASK_interface_init
   write(6,'(/,a,2(i2.2,a),i4.4)') ' Date: ',dateAndTime(3),'/',dateAndTime(2),'/', dateAndTime(1)
   write(6,'(a,2(i2.2,a),i2.2)')   ' Time: ',dateAndTime(5),':', dateAndTime(6),':', dateAndTime(7)
 
-  call MPI_Type_size(MPI_INTEGER,typeSize,mpi_err)
+  call MPI_Type_size(MPI_INTEGER,typeSize,err)
-  if (mpi_err /= 0) call quit(1)
+  if (err /= 0) call quit(1)
   if (typeSize*8 /= bit_size(0)) then
     write(6,'(a)') ' Mismatch between MPI and DAMASK integer' 
     call quit(1)
   endif
 
-  call MPI_Type_size(MPI_DOUBLE,typeSize,mpi_err)
+  call MPI_Type_size(MPI_DOUBLE,typeSize,err)
-  if (mpi_err /= 0) call quit(1)
+  if (err /= 0) call quit(1)
   if (typeSize*8 /= storage_size(0.0_pReal)) then
     write(6,'(a)') ' Mismatch between MPI and DAMASK real' 
     call quit(1)
   endif
 
   do i = 1, command_argument_count()
-    call get_command_argument(i,arg)
+    call get_command_argument(i,arg,status=err)
+    if (err /= 0) call quit(1)
     select case(trim(arg))                                                                          ! extract key
       case ('-h','--help')
         write(6,'(a)')  ' #######################################################################'
@@ -236,19 +237,20 @@ subroutine DAMASK_interface_init
         write(6,'(a,/)')'        Prints this message and exits'
         call quit(0)                                                                                ! normal Termination
       case ('-l', '--load', '--loadcase')
-        call get_command_argument(i+1,loadCaseArg)
+        call get_command_argument(i+1,loadCaseArg,status=err)
       case ('-g', '--geom', '--geometry')
-        call get_command_argument(i+1,geometryArg)
+        call get_command_argument(i+1,geometryArg,status=err)
       case ('-w', '--wd', '--workingdir', '--workingdirectory')
-        call get_command_argument(i+1,workingDirArg)
+        call get_command_argument(i+1,workingDirArg,status=err)
       case ('-r', '--rs', '--restart')
-        call get_command_argument(i+1,arg)
+        call get_command_argument(i+1,arg,status=err)
         read(arg,*,iostat=stat) interface_restartInc
         if (interface_restartInc < 0 .or. stat /=0) then
           write(6,'(/,a)') ' ERROR: Could not parse restart increment: '//trim(arg)
           call quit(1)
         endif
     end select
+    if (err /= 0) call quit(1)
   enddo
 
   if (len_trim(loadcaseArg) == 0 .or. len_trim(geometryArg) == 0) then

src/IO.f90

@@ -49,7 +49,7 @@ subroutine IO_init
 
   write(6,'(/,a)') ' <<<+-  IO init  -+>>>'; flush(6)
 
-  call unitTest
+  call selfTest
 
 end subroutine IO_init
 
@@ -696,7 +696,7 @@ end subroutine IO_warning
 !--------------------------------------------------------------------------------------------------
 !> @brief check correctness of some IO functions
 !--------------------------------------------------------------------------------------------------
-subroutine unitTest
+subroutine selfTest
 
   integer, dimension(:), allocatable :: chunkPos
   character(len=:),      allocatable :: str
@@ -745,6 +745,6 @@ subroutine unitTest
   str = IO_rmComment(' ab #')
   if (str /= ' ab'.or. len(str) /= 3)               call IO_error(0,ext_msg='IO_rmComment/7')
 
-end subroutine unitTest
+end subroutine selfTest
 
 end module IO

src/YAML_types.f90

@@ -180,7 +180,7 @@ subroutine YAML_types_init
 
   write(6,'(/,a)') ' <<<+-  YAML_types init  -+>>>'
 
-  call unitTest
+  call selfTest
 
 end subroutine YAML_types_init
 
@@ -188,7 +188,7 @@ end subroutine YAML_types_init
 !--------------------------------------------------------------------------------------------------
 !> @brief check correctness of some type bound procedures
 !--------------------------------------------------------------------------------------------------
-subroutine unitTest
+subroutine selfTest
 
   class(tNode), pointer  :: s1,s2
   allocate(tScalar::s1)
@@ -260,7 +260,7 @@ subroutine unitTest
      if(n%get_asString(1) /= 'True')                call IO_error(0,ext_msg='byIndex_asString')
   end block
 
-end subroutine unitTest
+end subroutine selfTest
 
 
 !---------------------------------------------------------------------------------------------------

src/lattice.f90

@@ -529,7 +529,7 @@ subroutine lattice_init
     lattice_DamageMobility(p) = config_phase(p)%getFloat('damage_mobility',defaultVal=0.0_pReal)
     ! SHOULD NOT BE PART OF LATTICE END
 
-    call unitTest
+    call selfTest
 
   enddo
 
@@ -1436,6 +1436,7 @@ function lattice_SchmidMatrix_slip(Nslip,structure,cOverA) result(SchmidMatrix)
       NslipMax    = BCT_NSLIPSYSTEM
       slipSystems = BCT_SYSTEMSLIP
     case default
+      allocate(NslipMax(0))
       call IO_error(137,ext_msg='lattice_SchmidMatrix_slip: '//trim(structure))
   end select
 
@@ -1485,6 +1486,7 @@ function lattice_SchmidMatrix_twin(Ntwin,structure,cOverA) result(SchmidMatrix)
       NtwinMax    = HEX_NTWINSYSTEM
       twinSystems = HEX_SYSTEMTWIN
     case default
+      allocate(NtwinMax(0))
       call IO_error(137,ext_msg='lattice_SchmidMatrix_twin: '//trim(structure))
   end select
 
@@ -1564,6 +1566,7 @@ function lattice_SchmidMatrix_cleavage(Ncleavage,structure,cOverA) result(Schmid
       NcleavageMax    = BCC_NCLEAVAGESYSTEM
       cleavageSystems = BCC_SYSTEMCLEAVAGE
     case default
+      allocate(NcleavageMax(0))
       call IO_error(137,ext_msg='lattice_SchmidMatrix_cleavage: '//trim(structure))
   end select
 
@@ -1919,6 +1922,7 @@ function coordinateSystem_slip(Nslip,structure,cOverA) result(coordinateSystem)
       NslipMax    = BCT_NSLIPSYSTEM
       slipSystems = BCT_SYSTEMSLIP
     case default
+      allocate(NslipMax(0))
       call IO_error(137,ext_msg='coordinateSystem_slip: '//trim(structure))
   end select
 
@@ -2291,7 +2295,7 @@ end function equivalent_mu
 !--------------------------------------------------------------------------------------------------
 !> @brief check correctness of some lattice functions
 !--------------------------------------------------------------------------------------------------
-subroutine unitTest
+subroutine selfTest
 
   real(pReal), dimension(:,:,:), allocatable :: CoSy
   real(pReal), dimension(:,:),   allocatable :: system
@@ -2320,6 +2324,6 @@ subroutine unitTest
   if(dNeq(lambda*0.5_pReal/(lambda+equivalent_mu(C,'reuss')),equivalent_nu(C,'reuss'),1.0e-12_pReal)) &
     call IO_error(0,ext_msg='equivalent_nu/reuss')
 
-end subroutine unitTest
+end subroutine selfTest
 
 end module lattice

src/math.f90

@@ -79,7 +79,7 @@ module math
 
 !---------------------------------------------------------------------------------------------------
  private :: &
-   unitTest
+   selfTest
 
 contains
 
@@ -113,7 +113,7 @@ subroutine math_init
 
   call random_seed(put = randInit)
 
-  call unitTest
+  call selfTest
 
 end subroutine math_init
 
@@ -1192,7 +1192,7 @@ end function math_clip
 !--------------------------------------------------------------------------------------------------
 !> @brief check correctness of some math functions
 !--------------------------------------------------------------------------------------------------
-subroutine unitTest
+subroutine selfTest
 
   integer, dimension(2,4) :: &
     sort_in_   = reshape([+1,+5,  +5,+6,  -1,-1,  +3,-2],[2,4])
@@ -1330,6 +1330,6 @@ subroutine unitTest
   if(dNeq0(math_LeviCivita(ijk(1),ijk(2),ijk(3))))&
     call IO_error(0,ext_msg='math_LeviCivita')
 
-end subroutine unitTest
+end subroutine selfTest
 
 end module math

src/mesh/mesh_mech_FEM.f90

@@ -17,11 +17,9 @@ module mesh_mech_FEM
   use prec
   use FEM_utilities
   use discretization_mesh
-  use IO
   use DAMASK_interface
   use numerics
   use FEM_quadrature
-  use FEsolving
   use homogenization
   use math
 

src/prec.f90

@@ -75,7 +75,7 @@ module prec
     emptyStringArray = [character(len=pStringLen)::]
 
   private :: &
-    unitTest
+    selfTest
 
 contains
 
@@ -94,7 +94,7 @@ subroutine prec_init
   write(6,'(a,e10.3)') '   Minimum value:        ',tiny(0.0_pReal)
   write(6,'(a,i3)')    '   Decimal precision:    ',precision(0.0_pReal)
 
-  call unitTest
+  call selfTest
 
 end subroutine prec_init
 
@@ -233,7 +233,7 @@ end function cNeq
 !--------------------------------------------------------------------------------------------------
 !> @brief check correctness of some prec functions
 !--------------------------------------------------------------------------------------------------
-subroutine unitTest
+subroutine selfTest
 
   integer, allocatable, dimension(:) :: realloc_lhs_test
   real(pReal), dimension(2) :: r
@@ -249,6 +249,6 @@ subroutine unitTest
   realloc_lhs_test = [1,2]
   if (any(realloc_lhs_test/=[1,2])) call quit(9000)
 
-end subroutine unitTest
+end subroutine selfTest
 
 end module prec

src/quaternions.f90

@@ -112,7 +112,7 @@ contains
 subroutine quaternions_init
 
   write(6,'(/,a)') ' <<<+-  quaternions init  -+>>>'; flush(6)
-  call unitTest
+  call selfTest
 
 end subroutine quaternions_init
 
@@ -457,7 +457,7 @@ end function inverse
 !--------------------------------------------------------------------------------------------------
 !> @brief check correctness of some quaternions functions
 !--------------------------------------------------------------------------------------------------
-subroutine unitTest
+subroutine selfTest
 
   real(pReal), dimension(4) :: qu
   type(quaternion)          :: q, q_2
@@ -524,7 +524,7 @@ subroutine unitTest
   endif
 #endif
 
-end subroutine unitTest
+end subroutine selfTest
 
 
 end module quaternions

src/rotations.f90

@@ -105,7 +105,7 @@ subroutine rotations_init
 
   call quaternions_init
   write(6,'(/,a)') ' <<<+-  rotations init  -+>>>'; flush(6)
-  call unitTest
+  call selfTest
 
 end subroutine rotations_init
 
@@ -1340,7 +1340,7 @@ end function GetPyramidOrder
 !--------------------------------------------------------------------------------------------------
 !> @brief check correctness of some rotations functions
 !--------------------------------------------------------------------------------------------------
-subroutine unitTest
+subroutine selfTest
 
   type(rotation)                  :: R
   real(pReal), dimension(4)       :: qu, ax, ro
@@ -1443,7 +1443,7 @@ subroutine unitTest
 
   enddo
 
-end subroutine unitTest
+end subroutine selfTest
 
 
 end module rotations