Update unsteady pressure computation

sdpm/src/sdpmSolver.cpp

@@ -522,20 +522,23 @@ void Solver::postUnsteady(size_t imd, size_t ifq, sdpmVectorXcd const &etau, sdp
         std::vector<std::vector<sdpmComplex>> cp012(elems.size(), std::vector<sdpmComplex>(3));
         for (size_t i = 0; i < elems.size(); ++i)
         {
-            // si = phi_t + 1/2 conv(u_k,u_k)_k - 1/2 u_inf^2
-            sdpmVectorXcd si(5);
-            si << 0., std::conj(-sdpmComplex(0, 1) * omega * emu(_rows[elems[i]])), -0.5 * ui.dot(ui), -sdpmComplex(0, 1) * omega * emu(_rows[elems[i]]), 0.;
+            // phi_t, phi_x
+            sdpmVector3cd phit(std::conj(-sdpmComplex(0, 1) * omega * emu(_rows[elems[i]])), 0., -sdpmComplex(0, 1) * omega * emu(_rows[elems[i]]));
+            sdpmVector3cd phix(_u1[imd][ifq][elems[i]->getId() - 1].conjugate()(0), _u[elems[i]->getId() - 1](0) - ui(0), _u1[imd][ifq][elems[i]->getId() - 1](0));
+            // 1 - phi_t
+            sdpmVectorXcd ecp(5);
+            ecp << 0., -2. * phit(0) / ui.dot(ui), 1., -2. * phit(2) / ui.dot(ui), 0.;
+            // phi_xx, phi_tt, phi_tx
+            ecp += mi * mi / ui.dot(ui) * (convolve(phix, phix) + convolve(phit, phit) + convolve(phix, phit) / ui.norm());
+            // conv(u_k,u_k)
             for (size_t j = 0; j < 3; ++j)
             {
                 sdpmVector3cd utrm(_u1[imd][ifq][elems[i]->getId() - 1].conjugate()(j), _u[elems[i]->getId() - 1](j), _u1[imd][ifq][elems[i]->getId() - 1](j));
-                si += 0.5 * convolve(utrm, utrm);
+                ecp -= convolve(utrm, utrm) / ui.dot(ui);
             }
-            // ecp = -2 / u_inf^2 * si + 1/u_inf^2/a_inf^2 * conv(si, si)
-            sdpmVectorXcd ecp(9);
-            ecp << 0., 0., -2 * si / ui.dot(ui), 0., 0.;
-            ecp += convolve(si, si) * mi * mi / (ui.dot(ui) * ui.dot(ui));
+            // extract 0th, 1st and 2nd harmonics
             for (size_t u = 0; u < 3; ++u)
-                cp012[i][u] = ecp(4 + u);
+                cp012[i][u] = ecp(2 + u);
             _cp1[imd][ifq][elems[i]->getId() - 1] = cp012[i][1];
         }
 

sdpm/tests/lann.py

@@ -69,10 +69,10 @@ def main():
     tests.add(CTest('CD0', sol.getDragCoeff(),  0.0020, 5e-2)) # MATLAB: 0.0031
     tests.add(CTest('CS0', sol.getSideforceCoeff(), 0.008, 5e-2))
     tests.add(CTest('CM0', sol.getMomentCoeff(), -0.079, 5e-2))
-    tests.add(CTest('Abs(CL1)', abs(sol.getUnsteadyLiftCoeff(0, 0)) / amp / math.pi, 1.54, 5e-2)) # MATLAB: 1.62
-    tests.add(CTest('Ang(CL1)', cmath.phase(sol.getUnsteadyLiftCoeff(0, 0)) * 180 / math.pi, 0.66, 5e-2)) # MATLAB: 2.6
-    tests.add(CTest('Abs(CD1)', abs(sol.getUnsteadyDragCoeff(0, 0)) / amp / math.pi, 0.038, 5e-2)) # MATLAB: 0.041
-    tests.add(CTest('Ang(CD1)', cmath.phase(sol.getUnsteadyDragCoeff(0, 0)) * 180 / math.pi, 14.6, 5e-2)) # MATLAB: 16.8
+    tests.add(CTest('Abs(CL1)', abs(sol.getUnsteadyLiftCoeff(0, 0)) / amp / math.pi, 1.61, 5e-2)) # MATLAB: 1.69
+    tests.add(CTest('Ang(CL1)', cmath.phase(sol.getUnsteadyLiftCoeff(0, 0)) * 180 / math.pi, 0.46, 5e-2)) # MATLAB: 0.64
+    tests.add(CTest('Abs(CD1)', abs(sol.getUnsteadyDragCoeff(0, 0)) / amp / math.pi, 0.036, 5e-2)) # MATLAB: 0.037
+    tests.add(CTest('Ang(CD1)', cmath.phase(sol.getUnsteadyDragCoeff(0, 0)) * 180 / math.pi, 16.4, 5e-2)) # MATLAB: 19.9
     tests.add(CTest('Abs(CS1)', abs(sol.getUnsteadySideforceCoeff(0, 0)) / amp / math.pi, 0.037, 5e-2))
     tests.add(CTest('Ang(CS1)', cmath.phase(sol.getUnsteadySideforceCoeff(0, 0)) * 180 / math.pi,  -10.0, 5e-2))
     tests.add(CTest('Abs(CM1)', abs(sol.getUnsteadyMomentCoeff(0, 0)) / amp / math.pi, 0.33, 5e-2))