test critical functions

Author: coding4Acause

tests/test.cpp

@@ -0,0 +1,30 @@
+#define CATCH_CONFIG_MAIN
+#include <catch2/catch.hpp>
+#include <Eigen/Dense>
+#include "geometry.h"
+
+TEST_CASE("nodal_coordinates_initial produces correct number of nodes", "[geometry]") {
+    int n = 10;
+    double c = 1.0, q = 0.2, p = 0.04, t_m = 0.12;
+    int trailing_edge_type = 0;
+    Eigen::VectorXd x0(n), y0(n);
+    
+    nodal_coordinates_initial(n, c, q, p, trailing_edge_type, t_m, x0, y0);
+    
+    REQUIRE(x0.size() == n);
+    REQUIRE(y0.size() == n);
+    REQUIRE(x0(0) == Approx(0.0).epsilon(0.01)); // Leading edge
+    REQUIRE(y0(0) == Approx(0.0).epsilon(0.01));
+}
+
+TEST_CASE("nodal_coordinates_initial handles invalid inputs", "[geometry]") {
+    int n = 0;
+    double c = 1.0, q = 0.2, p = 0.04, t_m = 0.12;
+    int trailing_edge_type = 0;
+    Eigen::VectorXd x0(n), y0(n);
+    
+    nodal_coordinates_initial(n, c, q, p, trailing_edge_type, t_m, x0, y0);
+    
+    REQUIRE(x0.size() == 0);
+    REQUIRE(y0.size() == 0);
+}

tests/test_geometry.cpp

@@ -0,0 +1,12 @@
+#include <catch2/catch.hpp>
+#include <Eigen/Dense>
+#include "InfluenceMatrix.h"
+
+TEST_CASE("influence_matrix returns correct size", "[influencematrix]") {
+    double x1 = 0.0, y1 = 0.0, x2 = 1.0, y2 = 0.0, x = 0.5, y = 0.1;
+    
+    Eigen::MatrixXd matrix = influence_matrix(x1, y1, x2, y2, x, y);
+    
+    REQUIRE(matrix.rows() == 2);
+    REQUIRE(matrix.cols() == 2);
+}

tests/test_influencematrix.cpp

@@ -0,0 +1,68 @@
+{
+  "__geometry_explain": {
+    "n": "Number of nodes along the airfoil surface (n-1 panels) [It can be even or odd]",
+    "ymc": "1st digit of NACA 4-digit series (camber percentage)",
+    "xmc": "2nd digit of NACA 4-digit series (location of max camber)",
+    "tmax": "Last two digits of NACA 4-digit series (thickness percentage)",
+    "trailing_edge_type": "1 = open trailing edge, else closed",
+    "c": "Chord length of the airfoil [m]"
+  },
+  "geometry": {
+    "n": 101,
+    "c": 0.1,
+    "ymc": 0.0,
+    "xmc": 0.0,
+    "tmax": 12.0,
+    "trailing_edge_type": 2
+  },
+  "__flow_explain": {
+    "rho": "Fluid density [kg/m^3]",
+    "mu": "Dynamic viscosity [Pa.s]",
+    "Re": "Reynolds number (dimensionless)",
+    "Uinf": "Freestream X-velocity [optional: null means auto-compute from Re]",
+    "Vinf": "Freestream Y-velocity (usually 0)"
+  },
+  "flow": {
+    "rho": 998.2,
+    "mu": 0.0010016,
+    "Re": 40000.0,
+    "Qinf": null,
+    "Vinf": 0.0
+  },
+  "__motion_explain": {
+    "k": "Reduced frequency (k = omega * c / (2 * Uinf))",
+    "h1": "Amplitude of plunge motion [m]",
+    "h0": "Mean plunge offset [m]",
+    "alpha1": "pitch amplitude angle in degrees",
+    "alpha0": "Mean pitch angle in degrees",
+    "phi_h": "Plunge phase angle in degrees",
+    "x_pitch": "pitching axis x coord",
+    "y_pitch": "pitching axis y coord"
+  },
+  "motion": {
+    "k": 1.2,
+    "h1": null,
+    "h0": 0.0,
+    "alpha1": null,
+    "alpha0": 0.0,
+    "phi_h": 0.0,
+    "x_pitch": null,
+    "y_pitch": null
+  },
+  "__simulation_explain": {
+    "wake": "0 = free wake, 1 = prescribed wake",
+    "tolerance": "Convergence tolerance for Newton iteration",
+    "epsilon": "Small number for perturbation / finite differences",
+    "ncycles": "Number of oscillation cycles to simulate",
+    "nsteps": "Number of time steps per cycle",
+    "z": "Number of points for stagnation streamline integration"
+  },
+  "simulation": {
+    "wake": 0,
+    "tolerance": 1e-6,
+    "epsilon": 1e-8,
+    "ncycles": 1,
+    "nsteps": 40,
+    "z": 200
+  }
+}

tests/test_input.json

@@ -0,0 +1,12 @@
+#include <catch2/catch.hpp>
+#include "kinematics.h"
+#include "constants.h"
+
+TEST_CASE("alpha_instantaneous computes correct angle", "[kinematics]") {
+    double alpha0 = 0.0, alpha1 = 15.0 * DEG2RAD, phi_alpha = 90.0 * DEG2RAD;
+    double t = 0.0, omega = 1.0;
+    
+    double alpha = alpha_instantaneous(alpha0, alpha1, phi_alpha, t, omega);
+    
+    REQUIRE(alpha == Approx(alpha0).epsilon(0.01)); // At t=0, expect alpha0
+}

tests/test_kinematics.cpp

@@ -0,0 +1,23 @@
+#include <catch2/catch.hpp>
+#include <Eigen/Dense>
+#include "NewtonRaphsonNonLinear.h"
+#include "geometry.h"
+#include "velocity.h"
+
+TEST_CASE("newton_raphson returns residuals of correct size", "[newtonraphson]") {
+    int n = 10;
+    double dt = 0.1, t = 0.0, lwp = 0.1, theta_wp = 0.0;
+    Eigen::VectorXd freestream(2), vtotal_wp_cp(2), x_pp(n), y_pp(n), x_cp(n-1), y_cp(n-1), l(n-1);
+    Eigen::VectorXd B_unsteady(n+1), gamma_unsteady(n+1), gamma_bound(n);
+    std::vector<double> gamma_wake_strength, gamma_wake_x_location, gamma_wake_y_location;
+    Eigen::VectorXd wake_panel_cp(2), wake_panel_normal(2);
+    Eigen::MatrixXd A_unsteady(n+1, n+1), unit_normal(n-1, 2), wake_panel_coordinates(2, 2);
+    double gamma_old = 0.0;
+    
+    freestream << 1.0, 0.0;
+    x_pp.setZero(); y_pp.setZero(); x_cp.setZero(); y_cp.setZero(); l.setOnes();
+    
+    Eigen::VectorXd residuals = newton_raphson(n, dt, t, lwp, theta_wp, freestream, vtotal_wp_cp, x_pp, y_pp, x_cp, y_cp, l, B_unsteady, gamma_unsteady, gamma_old, gamma_bound, gamma_wake_strength, gamma_wake_x_location, gamma_wake_y_location, wake_panel_cp, wake_panel_normal, A_unsteady, unit_normal, wake_panel_coordinates);
+    
+    REQUIRE(residuals.size() == 2);
+}

tests/test_newtonraphson.cpp

@@ -0,0 +1,15 @@
+#include <catch2/catch.hpp>
+#include <Eigen/Dense>
+#include "velocity.h"
+#include "constants.h"
+
+TEST_CASE("velocity_at_surface_of_the_body_inertial_frame computes correct velocity", "[velocity]") {
+    double Qinf = 1.0, x_pitch = 0.5, y_pitch = 0.0, h0 = 0.0, h1 = 0.25, phi_h = 0.0;
+    double alpha0 = 0.0, alpha1 = 0.0, phi_alpha = 0.0, t = 0.0, omega = 1.0;
+    double x = 0.5, y = 0.0;
+    
+    Eigen::Vector2d velocity = velocity_at_surface_of_the_body_inertial_frame(Qinf, x_pitch, y_pitch, h0, h1, phi_h, alpha0, alpha1, phi_alpha, t, omega, x, y);
+    
+    REQUIRE(velocity(0) == Approx(Qinf).epsilon(0.01)); // Freestream velocity
+    REQUIRE(velocity(1) == Approx(0.0).epsilon(0.01));
+}