more work on array ordering in new solver. (now working but not tested)

code/DAMASK_spectral_SolverAL.f90

@@ -55,7 +55,7 @@ module DAMASK_spectral_SolverAL
  
 !--------------------------------------------------------------------------------------------------
 ! common pointwise data
-  real(pReal), private, dimension(:,:,:,:,:), allocatable ::  F_lastInc, F_lambda_lastInc, P
+  real(pReal), private, dimension(:,:,:,:,:), allocatable ::  F_lastInc, F_lambda_lastInc
   real(pReal), private, dimension(:,:,:,:),   allocatable ::  coordinates
   real(pReal), private, dimension(:,:,:),     allocatable ::  temperature
  
@@ -113,11 +113,12 @@ module DAMASK_spectral_SolverAL
     implicit none
     
     integer(pInt) :: i,j,k
+    real(pReal), dimension(:,:,:,:,:), allocatable ::  P
     
     PetscErrorCode :: ierr_psc
     PetscObject :: dummy
     PetscMPIInt :: rank
-    PetscScalar, pointer :: xx_psc(:,:,:,:)
+    PetscScalar, pointer :: xx_psc(:,:,:,:), F(:,:,:,:), F_lambda(:,:,:,:)
     
     call Utilities_init()
     
@@ -153,27 +154,23 @@ module DAMASK_spectral_SolverAL
  !--------------------------------------------------------------------------------------------------
  ! init fields                 
     call DMDAVecGetArrayF90(da,solution_vec,xx_psc,ierr_psc)
-
+    F => xx_psc(0:8,:,:,:)
+    F_lambda => xx_psc(9:17,:,:,:)
     if (restartInc == 1_pInt) then                                                                     ! no deformation (no restart)
-      F_lastInc         = spread(spread(spread(math_I3,3,res(1)),4,res(2)),5,res(3))                           ! initialize to identity
+      F_lastInc         = spread(spread(spread(math_I3,3,res(1)),4,res(2)),5,res(3))                   ! initialize to identity
       F_lambda_lastInc = F_lastInc
-
+      F = reshape(F_lastInc,[9,res(1),res(2),res(3)])
-      xx_psc(0:8,:,:,:) = reshape(F_lastInc,[9,res(1),res(2),res(3)])
+      F_lambda = F
-
-      xx_psc(9:17,:,:,:) = xx_psc(0:8,:,:,:)
-
-    call flush(6)
       do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt, res(1)
         coordinates(i,j,k,1:3) = geomdim/real(res,pReal)*real([i,j,k],pReal) &
                                - geomdim/real(2_pInt*res,pReal)
       enddo; enddo; enddo
-      print*, 'init done12'
     elseif (restartInc > 1_pInt) then                                                                  ! using old values from file                                                      
       if (debugRestart) write(6,'(a,i6,a)') 'Reading values of increment ',&
                                                 restartInc - 1_pInt,' from file' 
       call IO_read_jobBinaryFile(777,'convergedSpectralDefgrad',&
                                                    trim(getSolverJobName()),size(F_lastInc))
-      read (777,rec=1) xx_psc(0:8,:,:,:)
+      read (777,rec=1) F
       close (777)
       call IO_read_jobBinaryFile(777,'convergedSpectralDefgrad_lastInc',&
                                                    trim(getSolverJobName()),size(F_lastInc))
@@ -181,7 +178,7 @@ module DAMASK_spectral_SolverAL
       close (777)
       call IO_read_jobBinaryFile(777,'convergedSpectralDefgradLambda',&
                                                    trim(getSolverJobName()),size(F_lambda_lastInc))
-      read (777,rec=1) xx_psc(9:17,:,:,:)
+      read (777,rec=1) F_lambda
       close (777)
       call IO_read_jobBinaryFile(777,'convergedSpectralDefgradLambda_lastInc',&
                                                    trim(getSolverJobName()),size(F_lambda_lastInc))
@@ -197,12 +194,9 @@ module DAMASK_spectral_SolverAL
       coordinates = 0.0 ! change it later!!!
     endif
    
-    call Utilities_constitutiveResponse(coordinates,reshape(xx_psc(0:8,:,:,:),shape(F_lastInc)),&
+    call Utilities_constitutiveResponse(coordinates,F,F,temperature,0.0_pReal,P,C,P_av,.false.,math_I3)
-                                                    reshape(xx_psc(0:8,:,:,:),shape(F_lastInc)),&
-                                                    temperature,0.0_pReal,P,C,P_av,.false.,math_I3)
-    print*, 'const response'
     call DMDAVecRestoreArrayF90(da,solution_vec,xx_psc,ierr_psc)
- print*, 'restored'
+
  !--------------------------------------------------------------------------------------------------
  ! reference stiffness
     if (restartInc == 1_pInt) then
@@ -253,19 +247,18 @@ module DAMASK_spectral_SolverAL
    real(pReal), intent(in) :: timeinc, timeinc_old, temperature_bc, guessmode
    type(boundaryCondition),      intent(in) :: P_BC,F_BC
    real(pReal), dimension(3,3), intent(in) :: rotation_BC
-   SNESConvergedReason reason
    real(pReal), dimension(3,3)            :: deltaF_aim, &
                                              F_aim_lab
 !--------------------------------------------------------------------------------------------------
 ! loop variables, convergence etc.
    real(pReal), dimension(3,3)            :: temp33_Real 
-   integer(pInt) :: ctr, i, j, k
  
 !--------------------------------------------------------------------------------------------------
 ! 
-   PetscScalar, pointer :: xx_psc(:,:,:,:)
+   PetscScalar, pointer :: xx_psc(:,:,:,:), F(:,:,:,:), F_lambda(:,:,:,:)
-   PetscErrorCode ierr_psc
+   PetscErrorCode ierr_psc   
-   
+   SNESConvergedReason reason
+
 !--------------------------------------------------------------------------------------------------
 ! restart information for spectral solver
    if (restartWrite) then
@@ -296,10 +289,10 @@ module DAMASK_spectral_SolverAL
 ! update local deformation gradient and coordinates
    deltaF_aim = math_rotate_backward33(deltaF_aim,rotation_BC)
    call DMDAVecGetArrayF90(da,solution_vec,xx_psc,ierr_psc)
-   call Utilities_forwardField(deltaF_aim,timeinc,timeinc_old,guessmode,F_lastInc, &
+   F => xx_psc(0:8,:,:,:)
-                               xx_psc(1:9,1:res(1),1:res(2),1:res(3)))
+   F_lambda => xx_psc(9:17,:,:,:)
-   !call Utilities_forwardField(deltaF_aim,timeinc,timeinc_old,guessmode,F_lambda_lastInc,&
+   call Utilities_forwardField(deltaF_aim,timeinc,timeinc_old,guessmode,F_lastInc,F)
-    !                           reshape(xx_psc(10:18,1:res(1),1:res(2),1:res(3)),shape(F_lambda_lastInc)))
+   call Utilities_forwardField(deltaF_aim,timeinc,timeinc_old,guessmode,F_lambda_lastInc,F_lambda)
    call DMDAVecRestoreArrayF90(da,solution_vec,xx_psc,ierr_psc)
    call deformed_fft(res,geomdim,math_rotate_backward33(F_aim,rotation_BC),1.0_pReal,F_lastInc,coordinates)
   
@@ -345,15 +338,22 @@ module DAMASK_spectral_SolverAL
    implicit none
 
    integer(pInt) :: i,j,k,l,  ctr
-   real(pReal), dimension (3,3) ::  temp33_real
+   real(pReal), dimension(3,3)            :: temp33_Real 
    
    DMDALocalInfo :: in(DMDA_LOCAL_INFO_SIZE)
-   PetscScalar :: x_scal(in(DMDA_LOCAL_INFO_DOF),XG_RANGE,YG_RANGE,ZG_RANGE)  
+   PetscScalar, target :: x_scal(3,3,2,XG_RANGE,YG_RANGE,ZG_RANGE)  
-   PetscScalar :: f_scal(in(DMDA_LOCAL_INFO_DOF),X_RANGE,Y_RANGE,Z_RANGE) 
+   PetscScalar, target :: f_scal(3,3,2,X_RANGE,Y_RANGE,Z_RANGE) 
+   PetscScalar, pointer :: F(:,:,:,:,:), F_lambda(:,:,:,:,:) 
+   PetscScalar, pointer :: residual_F(:,:,:,:,:), residual_F_lambda(:,:,:,:,:)
    PetscInt :: iter, nfuncs
    PetscObject :: dummy
    PetscErrorCode :: ierr_psc
 
+   F => x_scal(:,:,1,:,:,:)
+   F_lambda => x_scal(:,:,2,:,:,:)
+   residual_F => f_scal(:,:,1,:,:,:)
+   residual_F_lambda => f_scal(:,:,2,:,:,:)
+   
    call SNESGetNumberFunctionEvals(snes,nfuncs,ierr_psc)
    call SNESGetIterationNumber(snes,iter,ierr_psc)
   
@@ -367,9 +367,8 @@ module DAMASK_spectral_SolverAL
 
  !--------------------------------------------------------------------------------------------------
  ! evaluate constitutive response
-   call Utilities_constitutiveResponse(coordinates,F_lastInc,reshape(x_scal(1:9,1:res(1),1:res(2),1:res(3)),shape(F_lastInc)),&
+   call Utilities_constitutiveResponse(coordinates,F_lastInc,F,temperature,params%timeinc, &
-                                temperature,params%timeinc,&
+                                       residual_F,C,P_av,ForwardData,params%rotation_BC)
-                                 P,C,P_av,ForwardData,params%rotation_BC)
    ForwardData = .False.
    
 !--------------------------------------------------------------------------------------------------
@@ -381,37 +380,22 @@ module DAMASK_spectral_SolverAL
  !--------------------------------------------------------------------------------------------------
  ! doing Fourier transform
    field_real = 0.0_pReal
-   
-   do j = 1_pInt, 3_pInt; do i = 1_pInt, 3_pInt
-     ctr = 1_pInt
-     do k = 1_pInt, 3_pInt; do l = 1_pInt, 3_pInt
-       field_real(1:res(1),1:res(2),1:res(3),i,j) = field_real(1:res(1),1:res(2),1:res(3),i,j) + &
-                                                    C_scale(i,j,l,k)*(x_scal(9+ctr,1:res(1),1:res(2),1:res(3)) - &
-                                                                      x_scal(ctr,1:res(1),1:res(2),1:res(3)))
-       ! P_temp(i,j,1:res(1),1:res(2),1:res(3)) = P_temp(i,j,1:res(1),1:res(2),1:res(3)) + &
-                                                    ! S_scale(i,j,l,k)*P(l,k,1:res(1),1:res(2),1:res(3))
-       ctr = ctr + 1_pInt                                         
-     enddo; enddo
-   enddo; enddo
-   
    do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt, res(1)
-     P(1:3,1:3,i,j,k) = math_mul3333xx33(S_scale,P(1:3,1:3,i,j,k)) + math_I3
+     temp33_Real = math_mul3333xx33(S_scale,residual_F(1:3,1:3,i,j,k)) + math_I3
+     residual_F(1:3,1:3,i,j,k) = temp33_Real
+     field_real(i,j,k,1:3,1:3) = math_mul3333xx33(C_scale,F_lambda(1:3,1:3,i,j,k)-F(1:3,1:3,i,j,k))
    enddo; enddo; enddo
   
    call Utilities_forwardFFT()
    call Utilities_fourierConvolution(math_rotate_backward33(F_aim,params%rotation_BC)) 
    call Utilities_backwardFFT()
                             
-   f_scal(1:9,1:res(1),1:res(2),1:res(3)) = reshape(P,[9,res(1),res(2),res(3)]) - x_scal(10:18,1:res(1),1:res(2),1:res(3)) +& 
+   residual_F_lambda = F - reshape(field_real(1:res(1),1:res(2),1:res(3),1:3,1:3),&
-                       x_scal(1:9,1:res(1),1:res(2),1:res(3)) - &
+                                   [3,3,res(1),res(2),res(3)],order=[3,4,5,1,2])
-                       reshape(field_real(1:res(1),1:res(2),1:res(3),1:3,1:3),[9,res(1),res(2),res(3)],order=[2,3,4,1])
+   residual_F = residual_F - F_lambda + residual_F_lambda
-   f_scal(10:18,1:res(1),1:res(2),1:res(3)) = x_scal(1:9,1:res(1),1:res(2),1:res(3)) - &
-                                              reshape(field_real(1:res(1),1:res(2),1:res(3),1:3,1:3),&
-                                                      [9,res(1),res(2),res(3)],order=[2,3,4,1])
    
-   err_f = wgt*sqrt(sum(f_scal(10:18,1:res(1),1:res(2),1:res(3))**2.0_pReal))
+   err_f = wgt*sqrt(sum(residual_F_lambda**2.0_pReal))
-   err_p = wgt*sqrt(sum((f_scal( 1:9 ,1:res(1),1:res(2),1:res(3)) &
+   err_p = wgt*sqrt(sum((residual_F - residual_F_lambda)**2.0_pReal))
-                        -f_scal(10:18,1:res(1),1:res(2),1:res(3)))**2.0_pReal))
 
  end subroutine AL_formResidual
 
@@ -470,6 +454,7 @@ module DAMASK_spectral_SolverAL
   
    call VecDestroy(solution_vec,ierr_psc)
    call SNESDestroy(snes,ierr_psc)
+   call DMDestroy(da,ierr_psc)
    call PetscFinalize(ierr_psc)
    call Utilities_destroy()
 

code/DAMASK_spectral_Utilities.f90

@@ -357,6 +357,7 @@ subroutine Utilities_fourierConvolution(fieldAim)
 
  write(6,'(a)') ''
  write(6,'(a)') '... doing convolution .................'
+ write(6,'(a)') ''
   
 !--------------------------------------------------------------------------------------------------
 ! to the actual spectral method calculation (mechanical equilibrium)
@@ -536,6 +537,7 @@ subroutine Utilities_constitutiveResponse(coordinates,F_lastInc,F,temperature,ti
   
   write(6,'(a)') ''
   write(6,'(a)') '... evaluating constitutive response .................'
+  write(6,'(a)') ''
   
   if (ForwardData) then
     CPFEM_mode = 1_pInt