Strain calculation

The calculation of the strain tensor (pressure gradient) is concisely written in the subroutine del() . The notation follows the mathematical formula in equation (1). The only aberration from the mathematical notation becomes obvious where I use abbreviated indices for practical reasons. It is transparent which of the calculated components are evaluated in a staggered fashion. If conj==0 then subroutine del() will evaluate $\nabla^T$ in (1), meaning we perform an outer product. Displacement uu and strain ss are conjugates of each other. The actual derivative calculation is carried out in subroutine conv1dc() , which is a finite difference approximation to a partial derivative. Instead of calling a convolution we could have called a Fourier derivative, then we would perform pseudo-spectral modeling. Or we could alternate between them depending on which space axis is being calculated.

# apply symmetric Lagrangian derivative ( with 8 point convolution)
#
subroutine del (conj,uu,t,ss,bnd)
implicit none
#include "commons"
integer         conj,   t
WFIELD(uu)
WFIELDLAY(uu)
TENSOR(ss)
TENSORLAY(ss)
BOUND(bnd)
BOUNDLAY(bnd)
integer ivec,iten,direction
real scale

if(conj==0) {ss = 0.;        }
else        {uu(t,:,:,:,:)=0.}

do ivec=1,ncomp {
do iten=1,ntcomp{    direction=trudim(deriv(tcomp(iten),icomp(ivec)))
 if (direction>0){   scale = scaling(tcomp(iten),icomp(ivec))
  if (conj==0)  { stagger = stagc0(tcomp(iten),icomp(ivec)) }
  else          { stagger = stagc1(tcomp(iten),icomp(ivec)) }
    call conv1dc(conj,uu,t,ivec,ss,iten,direction,ddx(direction)*scale,bnd)
  }
}}

return
end