diff --git a/src/coreComponents/finiteElement/elementFormulations/Qk_Hexahedron_Lagrange_GaussLobatto.hpp b/src/coreComponents/finiteElement/elementFormulations/Qk_Hexahedron_Lagrange_GaussLobatto.hpp
index 8b1c73c1b85..0e9d7abb971 100644
--- a/src/coreComponents/finiteElement/elementFormulations/Qk_Hexahedron_Lagrange_GaussLobatto.hpp
+++ b/src/coreComponents/finiteElement/elementFormulations/Qk_Hexahedron_Lagrange_GaussLobatto.hpp
@@ -1234,7 +1234,7 @@ computeStiffnessxyTerm( localIndex const q,
   GL_BASIS::TensorProduct3D::multiIndex( q, qa, qb, qc );
   real64 B[6] = {0};
   real64 J[3][3] = {{0}};
-  computeBxyMatrix( qa, qb, qc, X, J, B ); // The only change!
+  computeBxyMatrix( qa, qb, qc, X, J, B );
   computeGradPhiBGradPhi( qa, qb, qc, B, func );
 }
 
@@ -1252,7 +1252,7 @@ computeStiffnesszTerm( localIndex const q,
   GL_BASIS::TensorProduct3D::multiIndex( q, qa, qb, qc );
   real64 B[6] = {0};
   real64 J[3][3] = {{0}};
-  computeBzMatrix( qa, qb, qc, X, J, B ); // The only change!
+  computeBzMatrix( qa, qb, qc, X, J, B );
   computeGradPhiBGradPhi( qa, qb, qc, B, func );
 }
 
diff --git a/src/coreComponents/physicsSolvers/wavePropagation/sem/acoustic/secondOrderEqn/anisotropic/AcousticVTIFletcherWaveEquationSEM.cpp b/src/coreComponents/physicsSolvers/wavePropagation/sem/acoustic/secondOrderEqn/anisotropic/AcousticVTIFletcherWaveEquationSEM.cpp
index 9621382fb97..deb236409e2 100644
--- a/src/coreComponents/physicsSolvers/wavePropagation/sem/acoustic/secondOrderEqn/anisotropic/AcousticVTIFletcherWaveEquationSEM.cpp
+++ b/src/coreComponents/physicsSolvers/wavePropagation/sem/acoustic/secondOrderEqn/anisotropic/AcousticVTIFletcherWaveEquationSEM.cpp
@@ -19,6 +19,7 @@
 
 #include "AcousticVTIFletcherWaveEquationSEM.hpp"
 #include "AcousticVTIFletcherWaveEquationSEMKernel.hpp"
+#include "AcousticVTIFletcherAdjointWaveEquationSEMKernel.hpp"
 
 #include "finiteElement/FiniteElementDiscretization.hpp"
 #include "fieldSpecification/FieldSpecificationManager.hpp"
diff --git a/src/coreComponents/physicsSolvers/wavePropagation/sem/acoustic/secondOrderEqn/anisotropic/AcousticVTIZhangWaveEquationSEM.cpp b/src/coreComponents/physicsSolvers/wavePropagation/sem/acoustic/secondOrderEqn/anisotropic/AcousticVTIZhangWaveEquationSEM.cpp
index 90acb707ba5..ab0232afd63 100644
--- a/src/coreComponents/physicsSolvers/wavePropagation/sem/acoustic/secondOrderEqn/anisotropic/AcousticVTIZhangWaveEquationSEM.cpp
+++ b/src/coreComponents/physicsSolvers/wavePropagation/sem/acoustic/secondOrderEqn/anisotropic/AcousticVTIZhangWaveEquationSEM.cpp
@@ -19,6 +19,7 @@
 
 #include "AcousticVTIZhangWaveEquationSEM.hpp"
 #include "AcousticVTIZhangWaveEquationSEMKernel.hpp"
+#include "AcousticVTIZhangAdjointWaveEquationSEMKernel.hpp"
 
 #include "finiteElement/FiniteElementDiscretization.hpp"
 #include "fieldSpecification/FieldSpecificationManager.hpp"
diff --git a/src/coreComponents/physicsSolvers/wavePropagation/sem/acoustic/shared/AcousticMatricesSEMKernel.hpp b/src/coreComponents/physicsSolvers/wavePropagation/sem/acoustic/shared/AcousticMatricesSEMKernel.hpp
index fceec37d5f8..d5e85ec18f0 100644
--- a/src/coreComponents/physicsSolvers/wavePropagation/sem/acoustic/shared/AcousticMatricesSEMKernel.hpp
+++ b/src/coreComponents/physicsSolvers/wavePropagation/sem/acoustic/shared/AcousticMatricesSEMKernel.hpp
@@ -321,13 +321,23 @@ struct AcousticMatricesSEM
               }
             }
             constexpr localIndex numNodesPerFace = FE_TYPE::numNodesPerFace;
+            real32 epsi = std::fabs( vti_epsilon[e] );
+            if( std::fabs( epsi ) < 1e-5 )
+              epsi = 0;
+            real32 delt = std::fabs( vti_delta[e] );
+            if( std::fabs( delt ) < 1e-5 )
+              delt = 0;
+            if( delt > epsi )
+              delt = epsi;
+            real32 sqrtEpsi = sqrt( 1 + 2 * epsi );
+            real32 sqrtDelta = sqrt( 1 + 2 * vti_delta[e] );
             if( lateralSurfaceFaceIndicator[f] == 1 )
             {
               // ABC coefficients updated to fit horizontal velocity
-              real32 alpha = 1.0 / (velocity[e] * density[e] * sqrt( 1+2*vti_epsilon[e] ));
+              real32 alpha = 1.0 / (velocity[e] * density[e] * sqrtEpsi);
               // VTI coefficients
-              real32 vti_p_xy  = (1+2*vti_epsilon[e]);
-              real32 vti_qp_xy = sqrt( 1+2*vti_delta[e] );
+              real32 vti_p_xy  = 1 + 2 * epsi;
+              real32 vti_qp_xy = sqrtDelta;
 
               for( localIndex q = 0; q < numNodesPerFace; ++q )
               {
@@ -344,7 +354,7 @@ struct AcousticMatricesSEM
               // ABC coefficients updated to fit horizontal velocity
               real32 alpha = 1.0 / (velocity[e] * density[e]);
               // VTI coefficients
-              real32 vti_pq_z = sqrt( 1+2*vti_delta[e] );
+              real32 vti_pq_z = sqrtDelta;
               real32 vti_q_z  = 1;
               for( localIndex q = 0; q < numNodesPerFace; ++q )
               {