diff --git a/.github/workflows/build.yml b/.github/workflows/build.yml index e7ea76f7..f23a6a5f 100644 --- a/.github/workflows/build.yml +++ b/.github/workflows/build.yml @@ -206,4 +206,4 @@ jobs: shell: bash # Execute tests defined by the CMake configuration. # See https://cmake.org/cmake/help/latest/manual/ctest.1.html for more detail - run: ctest -C $BUILD_TYPE -V \ No newline at end of file + run: ctest -C $BUILD_TYPE --output-on-failure diff --git a/sandbox/CMakeLists.txt b/sandbox/CMakeLists.txt index 5f13f075..cfeddb25 100644 --- a/sandbox/CMakeLists.txt +++ b/sandbox/CMakeLists.txt @@ -14,24 +14,26 @@ function(create_sandbox source_list) endfunction(create_sandbox) set(SANDBOX_SRCS - advection2d - steady_vortex - steady_vortex_adjoint - unsteady_vortex - viscous_shock - airfoil_steady - airfoil_chaotic - magnetostatic_box + #advection2d + #steady_vortex + #steady_vortex_adjoint + #unsteady_vortex + #viscous_shock + #airfoil_steady + #airfoil_chaotic + #magnetostatic_box magnetostatic_motor - magnetostatic_wire - navier_stokes_mms + #magnetostatic_wire + #navier_stokes_mms thermal_square - # joule_wire - joule_box - mesh_move - mesh_move_2 - surface_distance - joule_wire + #joule_wire + #joule_box + #mesh_move + #mesh_move_2 + #surface_distance + #rans_freestream + rans_walltest + #rans_mms ) create_sandbox("${SANDBOX_SRCS}") diff --git a/sandbox/rans_freestream.cpp b/sandbox/rans_freestream.cpp new file mode 100644 index 00000000..522a2790 --- /dev/null +++ b/sandbox/rans_freestream.cpp @@ -0,0 +1,214 @@ +// set this const expression to true in order to use entropy variables for state (doesn't work for rans) +constexpr bool entvar = false; + +#include +#include "adept.h" + +#include "mfem.hpp" +#include "rans.hpp" +#include +#include + +using namespace std; +using namespace mfem; +using namespace mach; + +std::default_random_engine gen(std::random_device{}()); +std::uniform_real_distribution uniform_rand(0.0,1.0); + +static double pert_fs; +static double mu; +static double mach_fs; +static double aoa_fs; +static int iroll; +static int ipitch; +static double chi_fs; +static bool mms; + +/// \brief Defines the random function for the jacobian check +/// \param[in] x - coordinate of the point at which the state is needed +/// \param[out] u - conservative + SA variables stored as a 5-vector +void pert(const Vector &x, Vector& p); + +/// \brief Defines the exact solution for the rans freestream problem +/// \param[in] x - coordinate of the point at which the state is needed +/// \param[out] u - conservative + SA variables stored as a 5-vector +void uexact(const Vector &x, Vector& u); + +/// \brief Defines a perturbed solution for the rans freestream problem +/// \param[in] x - coordinate of the point at which the state is needed +/// \param[out] u - conservative + SA variables stored as a 5-vector +void uinit_pert(const Vector &x, Vector& u); + +void uinit_pert_mms(const Vector &x, Vector& u); + + +int main(int argc, char *argv[]) +{ + const char *options_file = "rans_freestream_options.json"; + + // Initialize MPI + int num_procs, rank; + MPI_Init(&argc, &argv); + MPI_Comm_size(MPI_COMM_WORLD, &num_procs); + MPI_Comm_rank(MPI_COMM_WORLD, &rank); + ostream *out = getOutStream(rank); + + // Parse command-line options + OptionsParser args(argc, argv); + int degree = 2.0; + int nx = 10; + int ny = 10; + args.AddOption(&nx, "-nx", "--numx", + "Number of elements in x direction"); + args.AddOption(&ny, "-ny", "--numy", + "Number of elements in y direction"); + args.Parse(); + if (!args.Good()) + { + args.PrintUsage(*out); + return 1; + } + + try + { + // construct the mesh + string opt_file_name(options_file); + nlohmann::json file_options; + std::ifstream opts(opt_file_name); + opts >> file_options; + pert_fs = 1.0 + file_options["init-pert"].template get(); + mu = file_options["flow-param"]["mu"].template get(); + mach_fs = file_options["flow-param"]["mach"].template get(); + aoa_fs = file_options["flow-param"]["aoa"].template get()*M_PI/180; + iroll = file_options["flow-param"]["roll-axis"].template get(); + ipitch = file_options["flow-param"]["pitch-axis"].template get(); + chi_fs = file_options["flow-param"]["chi"].template get(); + mms = file_options["flow-param"]["rans-mms"].template get(); + + // generate a simple tri mesh and a distance function + std::unique_ptr smesh(new Mesh(nx, ny, + Element::TRIANGLE, true /* gen. edges */, 1.0, + 1.0, true)); + *out << "Number of elements " << smesh->GetNE() <<'\n'; + + // construct the solver and set initial conditions + auto solver = createSolver>(opt_file_name, + move(smesh)); + if (mms) + solver->setInitialCondition(uinit_pert_mms); + else + solver->setInitialCondition(uinit_pert); + solver->printSolution("rans_init", 0); + solver->printResidual("rans_res_init", 0); + + + // get the initial density error + double l2_error_init = (static_cast&>(*solver) + .calcConservativeVarsL2Error(uexact, 0)); + *out << "\n|| rho_h - rho ||_{L^2} init = " << l2_error_init << endl; + + double res_error = solver->calcResidualNorm(); + *out << "\ninitial residual norm = " << res_error << endl; + solver->checkJacobian(pert); + solver->solveForState(); + solver->printSolution("rans_final",0); + // get the final density error + double l2_error_final = (static_cast&>(*solver) + .calcConservativeVarsL2Error(uexact, 0)); + res_error = solver->calcResidualNorm(); + + *out << "\nfinal residual norm = " << res_error; + *out << "\n|| rho_h - rho ||_{L^2} final = " << l2_error_final << endl; + } + catch (MachException &exception) + { + exception.print_message(); + } + catch (std::exception &exception) + { + cerr << exception.what() << endl; + } + +#ifdef MFEM_USE_PETSC + MFEMFinalizePetsc(); +#endif + + MPI_Finalize(); +} + +// perturbation function used to check the jacobian in each iteration +void pert(const Vector &x, Vector& p) +{ + p.SetSize(5); + for (int i = 0; i < 5; i++) + { + p(i) = 2.0 * uniform_rand(gen) - 1.0; + } +} + +// Exact solution; same as freestream bc +void uexact(const Vector &x, Vector& q) +{ + // q.SetSize(4); + // Vector u(4); + q.SetSize(5); + Vector u(5); + + u = 0.0; + u(0) = 1.0; + u(1) = u(0)*mach_fs*cos(aoa_fs); + u(2) = u(0)*mach_fs*sin(aoa_fs); + u(3) = 1/(euler::gamma*euler::gami) + 0.5*mach_fs*mach_fs; + u(4) = chi_fs*mu; + + if (entvar == false) + { + q = u; + } + else + { + throw MachException("No entvar for this"); + } +} + +// initial guess perturbed from exact +void uinit_pert(const Vector &x, Vector& q) +{ + q.SetSize(5); + Vector u(5); + + u = 0.0; + u(0) = pert_fs*1.0; + u(1) = u(0)*mach_fs*cos(aoa_fs); + u(2) = u(0)*mach_fs*sin(aoa_fs); + u(3) = pert_fs*1/(euler::gamma*euler::gami) + 0.5*mach_fs*mach_fs; + u(4) = pert_fs*chi_fs*mu; + + q = u; +} + +// initial guess perturbed from exact +void uinit_pert_mms(const Vector &x, Vector& q) +{ + const double rho0 = 1.0; + const double rhop = 0.05; + const double U0 = 0.5; + const double Up = 0.05; + const double T0 = 1.0; + const double Tp = 0.05; + const double chif = 100.0; + const double mu = 1.0; + q.SetSize(5); + q(0) = rho0 + rhop*pow(sin(M_PI*x(0)),2)*sin(M_PI*x(1)); + q(1) = 4.0*U0*x(1)*(1.0 - x(1)) + Up*sin(2 * M_PI * x(1)) * pow(sin(M_PI * x(0)),2); + q(2) = -Up*pow(sin(2 * M_PI * x(0)),2) * sin(M_PI * x(1)); + double T = T0 + Tp*(pow(x(0), 4) - 2 * pow(x(0), 3) + pow(x(0), 2) + + pow(x(1), 4) - 2 * pow(x(1), 3) + pow(x(1), 2)); + double p = q(0)*T; // T is nondimensionalized by 1/(R*a_infty^2) + q(3) = p/euler::gami + 0.5*q(0)*(q(1)*q(1) + q(2)*q(2)); + q(1) *= q(0); + q(2) *= q(0); + + q(4) = q(0)*chif*mu*x(1); +} diff --git a/sandbox/rans_freestream_options.json b/sandbox/rans_freestream_options.json new file mode 100644 index 00000000..7201eb7f --- /dev/null +++ b/sandbox/rans_freestream_options.json @@ -0,0 +1,52 @@ +{ + "flow-param": { + "mach": 1.0, + "aoa": 0.3, + "roll-axis": 0, + "pitch-axis": 1, + "chi": 3.0, + "Re": 1000.0, + "Pr": 0.75, + "mu": 1.0, + "rans-mms": false, + "sa-consts": [0.1355, 0.622, 0.666666666666667, 0.41, 0.3, 2, 7.1, 1.2, 0.5, 10, 16, 0.7, 0.9], + "sa-srcs": [1.0, 1.0] + }, + "space-dis": { + "degree": 1, + "lps-coeff": 1.0, + "basis-type": "csbp" + }, + "time-dis": { + "steady": true, + "steady-abstol": 1e-12, + "steady-restol": 1e-10, + "const-cfl": true, + "ode-solver": "PTC", + "t-final": 100, + "dt": 1e6, + "cfl": 1.0 + }, + "nonlin-solver": { + "abstol": 1e-12, + "maxiter": 100, + "printlevel": 1, + "reltol": 1e-2, + "type": "newton" + }, + "lin-solver": { + "reltol": 1e-2, + "abstol": 1e-12, + "printlevel": 0, + "maxiter": 100 + }, + "bcs": { + "far-field": [1, 1, 1, 1] + }, + "wall-func": { + "type": "const", + "val": 1.0 + }, + "file-names": "rans_free", + "init-pert": 0.0 +} diff --git a/sandbox/rans_mms.cpp b/sandbox/rans_mms.cpp new file mode 100644 index 00000000..c7e41a38 --- /dev/null +++ b/sandbox/rans_mms.cpp @@ -0,0 +1,214 @@ +// set this const expression to true in order to use entropy variables for state (doesn't work for rans) +constexpr bool entvar = false; + +#include +#include "adept.h" + +#include "mfem.hpp" +#include "rans.hpp" +#include +#include + +using namespace std; +using namespace mfem; +using namespace mach; + +std::default_random_engine gen(std::random_device{}()); +std::uniform_real_distribution uniform_rand(0.0,1.0); + +static double pert_fs; +static double mu; +static double mach_fs; +static double aoa_fs; +static int iroll; +static int ipitch; +static double chi_fs; +static bool mms; + +/// \brief Defines the random function for the jacobian check +/// \param[in] x - coordinate of the point at which the state is needed +/// \param[out] u - conservative + SA variables stored as a 5-vector +void pert(const Vector &x, Vector& p); + +/// \brief Defines the exact solution for the rans freestream problem +/// \param[in] x - coordinate of the point at which the state is needed +/// \param[out] u - conservative + SA variables stored as a 5-vector +void uexact(const Vector &x, Vector& u); + +/// \brief Defines a perturbed solution for the rans freestream problem +/// \param[in] x - coordinate of the point at which the state is needed +/// \param[out] u - conservative + SA variables stored as a 5-vector +void uinit_pert(const Vector &x, Vector& u); + +void uinit_pert_mms(const Vector &x, Vector& u); + + +int main(int argc, char *argv[]) +{ + const char *options_file = "rans_mms_options.json"; + + // Initialize MPI + int num_procs, rank; + MPI_Init(&argc, &argv); + MPI_Comm_size(MPI_COMM_WORLD, &num_procs); + MPI_Comm_rank(MPI_COMM_WORLD, &rank); + ostream *out = getOutStream(rank); + + // Parse command-line options + OptionsParser args(argc, argv); + int degree = 2.0; + int nx = 10; + int ny = 10; + args.AddOption(&nx, "-nx", "--numx", + "Number of elements in x direction"); + args.AddOption(&ny, "-ny", "--numy", + "Number of elements in y direction"); + args.Parse(); + if (!args.Good()) + { + args.PrintUsage(*out); + return 1; + } + + try + { + // construct the mesh + string opt_file_name(options_file); + nlohmann::json file_options; + std::ifstream opts(opt_file_name); + opts >> file_options; + pert_fs = 1.0 + file_options["init-pert"].template get(); + mu = file_options["flow-param"]["mu"].template get(); + mach_fs = file_options["flow-param"]["mach"].template get(); + aoa_fs = file_options["flow-param"]["aoa"].template get()*M_PI/180; + iroll = file_options["flow-param"]["roll-axis"].template get(); + ipitch = file_options["flow-param"]["pitch-axis"].template get(); + chi_fs = file_options["flow-param"]["chi"].template get(); + mms = file_options["flow-param"]["sa-mms"].template get(); + + // generate a simple tri mesh and a distance function + std::unique_ptr smesh(new Mesh(nx, ny, + Element::TRIANGLE, true /* gen. edges */, 1.0, + 1.0, true)); + *out << "Number of elements " << smesh->GetNE() <<'\n'; + + // construct the solver and set initial conditions + auto solver = createSolver>(opt_file_name, + move(smesh)); + if (mms) + solver->setInitialCondition(uinit_pert_mms); + else + solver->setInitialCondition(uinit_pert); + solver->printSolution("rans_mms_init", 0); + solver->printResidual("rans_res_mms_init", 0); + + + // get the initial density error + double l2_error_init = (static_cast&>(*solver) + .calcConservativeVarsL2Error(uexact, 0)); + *out << "\n|| rho_h - rho ||_{L^2} init = " << l2_error_init << endl; + + double res_error = solver->calcResidualNorm(); + *out << "\ninitial residual norm = " << res_error << endl; + solver->checkJacobian(pert); + solver->solveForState(); + solver->printSolution("rans_mms_final",0); + // get the final density error + double l2_error_final = (static_cast&>(*solver) + .calcConservativeVarsL2Error(uexact, 0)); + res_error = solver->calcResidualNorm(); + + *out << "\nfinal residual norm = " << res_error; + *out << "\n|| rho_h - rho ||_{L^2} final = " << l2_error_final << endl; + } + catch (MachException &exception) + { + exception.print_message(); + } + catch (std::exception &exception) + { + cerr << exception.what() << endl; + } + +#ifdef MFEM_USE_PETSC + MFEMFinalizePetsc(); +#endif + + MPI_Finalize(); +} + +// perturbation function used to check the jacobian in each iteration +void pert(const Vector &x, Vector& p) +{ + p.SetSize(5); + for (int i = 0; i < 5; i++) + { + p(i) = 2.0 * uniform_rand(gen) - 1.0; + } +} + +// Exact solution; same as freestream bc +void uexact(const Vector &x, Vector& q) +{ + // q.SetSize(4); + // Vector u(4); + q.SetSize(5); + Vector u(5); + + u = 0.0; + u(0) = 1.0; + u(1) = u(0)*mach_fs*cos(aoa_fs); + u(2) = u(0)*mach_fs*sin(aoa_fs); + u(3) = 1/(euler::gamma*euler::gami) + 0.5*mach_fs*mach_fs; + u(4) = chi_fs*mu; + + if (entvar == false) + { + q = u; + } + else + { + throw MachException("No entvar for this"); + } +} + +// initial guess perturbed from exact +void uinit_pert(const Vector &x, Vector& q) +{ + q.SetSize(5); + Vector u(5); + + u = 0.0; + u(0) = pert_fs*1.0; + u(1) = u(0)*mach_fs*cos(aoa_fs); + u(2) = u(0)*mach_fs*sin(aoa_fs); + u(3) = pert_fs*1/(euler::gamma*euler::gami) + 0.5*mach_fs*mach_fs; + u(4) = pert_fs*chi_fs*mu; + + q = u; +} + +// initial guess perturbed from exact +void uinit_pert_mms(const Vector &x, Vector& q) +{ + const double rho0 = 1.0; + const double rhop = 0.05; + const double U0 = 0.5; + const double Up = 0.05; + const double T0 = 1.0; + const double Tp = 0.05; + const double chif = 100.0; + const double mu = 1.0; + q.SetSize(5); + q(0) = rho0 + rhop*pow(sin(M_PI*x(0)),2)*sin(M_PI*x(1)); + q(1) = 4.0*U0*x(1)*(1.0 - x(1)) + Up*sin(2 * M_PI * x(1)) * pow(sin(M_PI * x(0)),2); + q(2) = -Up*pow(sin(2 * M_PI * x(0)),2) * sin(M_PI * x(1)); + double T = T0 + Tp*(pow(x(0), 4) - 2 * pow(x(0), 3) + pow(x(0), 2) + + pow(x(1), 4) - 2 * pow(x(1), 3) + pow(x(1), 2)); + double p = q(0)*T; // T is nondimensionalized by 1/(R*a_infty^2) + q(3) = p/euler::gami + 0.5*q(0)*(q(1)*q(1) + q(2)*q(2)); + q(1) *= q(0); + q(2) *= q(0); + + q(4) = q(0)*chif*mu*x(1); +} diff --git a/sandbox/rans_mms_options.json b/sandbox/rans_mms_options.json new file mode 100644 index 00000000..05c94606 --- /dev/null +++ b/sandbox/rans_mms_options.json @@ -0,0 +1,52 @@ +{ + "flow-param": { + "mach": 1.0, + "aoa": 0.3, + "roll-axis": 0, + "pitch-axis": 1, + "chi": 3.0, + "Re": 1.0, + "Pr": 0.75, + "mu": 1.0, + "sa-mms": true, + "sa-consts": [0.1355, 0.622, 0.666666666666667, 0.41, 0.3, 2, 7.1, 1.2, 0.5, 10, 16, 0.7, 0.9], + "sa-srcs": [1.0, 1.0] + }, + "space-dis": { + "degree": 1, + "lps-coeff": 1.0, + "basis-type": "csbp" + }, + "time-dis": { + "steady": true, + "steady-abstol": 1e-12, + "steady-restol": 1e-10, + "const-cfl": true, + "ode-solver": "PTC", + "t-final": 100, + "dt": 1e6, + "cfl": 1.0 + }, + "nonlin-solver": { + "abstol": 1e-12, + "maxiter": 100, + "printlevel": 1, + "reltol": 1e-2, + "type": "newton" + }, + "lin-solver": { + "reltol": 1e-2, + "abstol": 1e-12, + "printlevel": 0, + "maxiter": 100 + }, + "bcs": { + "far-field": [1, 1, 1, 1] + }, + "wall-func": { + "type": "const", + "val": 1.0 + }, + "file-names": "rans_mms", + "init-pert": 0.0 +} diff --git a/sandbox/rans_walltest.cpp b/sandbox/rans_walltest.cpp new file mode 100644 index 00000000..6f0d0401 --- /dev/null +++ b/sandbox/rans_walltest.cpp @@ -0,0 +1,377 @@ +// set this const expression to true in order to use entropy variables for state (doesn't work for rans) +constexpr bool entvar = false; + +#include +#include +#include +#include "adept.h" + +#include "mfem.hpp" +#include "euler.hpp" +#include "rans.hpp" +#include +#include + +using namespace std; +using namespace mfem; +using namespace mach; + +std::default_random_engine gen(std::random_device{}()); +std::uniform_real_distribution uniform_rand(0.0,1.0); + +static double pert_fs; +static double mu; +static double mach_fs; +static double aoa_fs; +static int iroll; +static int ipitch; +static double chi_fs; + +static double m_offset; +static double m_coeff; +static int m_x; +static int m_y; + +/// \brief Defines the random function for the jacobian check +/// \param[in] x - coordinate of the point at which the state is needed +/// \param[out] u - conservative + SA variables stored as a 5-vector +void pert(const Vector &x, Vector& p); + +/// \brief Defines the random function for the jacobian check +/// \param[in] x - coordinate of the point at which the state is needed +/// \param[out] u - conservative variables stored as a 4-vector +void pert_ns(const Vector &x, Vector& p); + +/// \brief Defines the exact solution for the rans freestream problem +/// \param[in] x - coordinate of the point at which the state is needed +/// \param[out] u - conservative + SA variables stored as a 5-vector +void uexact(const Vector &x, Vector& u); + +/// \brief Defines a perturbed solution for the rans freestream problem +/// \param[in] x - coordinate of the point at which the state is needed +/// \param[out] u - conservative + SA variables stored as a 5-vector +void uinit_pert(const Vector &x, Vector& u); + +/// \brief Defines a perturbed solution for the ns freestream problem +/// \param[in] x - coordinate of the point at which the state is needed +/// \param[out] u - conservative variables stored as a 4-vector +void uinit_pert_ns(const Vector &x, Vector& u); + +/// \brief Defines a pseudo wall-distance function +/// \param[in] x - coordinate of the point at which the distance is needed +double walldist(const Vector &x); + +/// Generate a quadrilateral wall mesh, more dense near the wall +/// \param[in] num_x - number of nodes in x +/// \param[in] num_y - number of nodes in y +std::unique_ptr buildWalledMesh(int num_x, + int num_y); + +int main(int argc, char *argv[]) +{ + const char *options_file = "rans_walltest_options.json"; + + // Initialize MPI + int num_procs, rank; + MPI_Init(&argc, &argv); + MPI_Comm_size(MPI_COMM_WORLD, &num_procs); + MPI_Comm_rank(MPI_COMM_WORLD, &rank); + ostream *out = getOutStream(rank); + + // Parse command-line options + OptionsParser args(argc, argv); + int degree = 2.0; + int nx = 14; + int ny = 10; + args.AddOption(&nx, "-nx", "--numx", + "Number of elements in x direction"); + args.AddOption(&ny, "-ny", "--numy", + "Number of elements in y direction"); + args.Parse(); + if (!args.Good()) + { + args.PrintUsage(*out); + return 1; + } + + try + { + // construct the mesh + string opt_file_name(options_file); + nlohmann::json file_options; + std::ifstream opts(opt_file_name); + opts >> file_options; + pert_fs = 1.0 + file_options["init-pert"].template get(); + mu = file_options["flow-param"]["mu"].template get(); + mach_fs = file_options["flow-param"]["mach"].template get(); + aoa_fs = file_options["flow-param"]["aoa"].template get()*M_PI/180; + iroll = file_options["flow-param"]["roll-axis"].template get(); + ipitch = file_options["flow-param"]["pitch-axis"].template get(); + chi_fs = file_options["flow-param"]["chi"].template get(); + + m_offset = file_options["mesh"]["offset"].template get(); //thickness of elements on wall + m_coeff = file_options["mesh"]["coeff"].template get(); + m_x = file_options["mesh"]["num-x"].template get(); + m_y = file_options["mesh"]["num-y"].template get(); + + string file = file_options["file-names"].template get(); + stringstream fileinit; stringstream filefinal; + fileinit << file << "_init"; + filefinal << file << "_final"; + + std::vector u_ns(1); + std::vector u_rans(1); + + // generate the walled mesh, and its distance function + std::unique_ptr smesh = buildWalledMesh(m_x, m_y); + + if (file_options["try-rans"].template get()) + { + // construct the solver and set initial conditions + auto solver = createSolver>(opt_file_name, + move(smesh)); + solver->setInitialCondition(uinit_pert); + solver->printSolution(fileinit.str(), 0); + + double res_error = solver->calcResidualNorm(); + *out << "\ninitial rans residual norm = " << res_error << endl; + //solver->checkJacobian(pert); + solver->solveForState(); + solver->printSolution(filefinal.str(),0); + // get the final density error + res_error = solver->calcResidualNorm(); + + out->precision(15); + *out << "\nfinal rans residual norm = " << res_error; + double drag = abs(solver->calcOutput("drag")); + *out << "\nDrag = " << drag << endl; + } + + } + catch (MachException &exception) + { + exception.print_message(); + } + catch (std::exception &exception) + { + cerr << exception.what() << endl; + } + +#ifdef MFEM_USE_PETSC + MFEMFinalizePetsc(); +#endif + + MPI_Finalize(); +} + +// perturbation function used to check the jacobian in each iteration +void pert(const Vector &x, Vector& p) +{ + p.SetSize(5); + for (int i = 0; i < 5; i++) + { + p(i) = 2.0 * uniform_rand(gen) - 1.0; + } +} + +// perturbation function used to check the jacobian in each iteration +void pert_ns(const Vector &x, Vector& p) +{ + p.SetSize(4); + for (int i = 0; i < 4; i++) + { + p(i) = 2.0 * uniform_rand(gen) - 1.0; + } +} + +// Exact solution; same as freestream bc +void uexact(const Vector &x, Vector& q) +{ + // q.SetSize(4); + // Vector u(4); + q.SetSize(5); + Vector u(5); + + u = 0.0; + u(0) = 1.0; + u(1) = u(0)*mach_fs*cos(aoa_fs); + u(2) = u(0)*mach_fs*sin(aoa_fs); + u(3) = 1/(euler::gamma*euler::gami) + 0.5*mach_fs*mach_fs; + u(4) = chi_fs*mu; + + if (entvar == false) + { + q = u; + } + else + { + throw MachException("No entvar for this"); + } +} + +// initial guess perturbed from exact +void uinit_pert(const Vector &x, Vector& q) +{ + // q.SetSize(4); + // Vector u(4); + q.SetSize(5); + Vector u(5); + + u = 0.0; + u(0) = pert_fs*1.0; + u(1) = u(0)*mach_fs*cos(aoa_fs); + u(2) = u(0)*mach_fs*sin(aoa_fs); + u(3) = pert_fs*1/(euler::gamma*euler::gami) + 0.5*mach_fs*mach_fs; + u(4) = pert_fs*chi_fs*abs(mu); + + // if(x(1) == 0.0) + // { + // u(1) = 1e-10; + // } + + q = u; +} + +// initial guess perturbed from exact +void uinit_pert_ns(const Vector &x, Vector& q) +{ + // q.SetSize(4); + // Vector u(4); + q.SetSize(4); + Vector u(4); + + u = 0.0; + u(0) = pert_fs*1.0; + u(1) = u(0)*mach_fs*cos(aoa_fs); + u(2) = u(0)*mach_fs*sin(aoa_fs); + u(3) = pert_fs*1/(euler::gamma*euler::gami) + 0.5*mach_fs*mach_fs; + // if(x(1) == 0.0) + // { + // u(1) = 1e-10; + // } + + q = u; +} + +std::unique_ptr buildWalledMesh(int num_x, int num_y) +{ + auto mesh_ptr = unique_ptr(new Mesh(num_x, num_y, + Element::TRIANGLE, true /* gen. edges */, + 2.33333, 1.0, true)); + + // Randomly reorder elements + int numelem = mesh_ptr->GetNE(); + unsigned seed = std::chrono::system_clock::now().time_since_epoch().count(); + Array randelem(numelem); + for (int n = 0; n < numelem; n++) + randelem[n] = n; + //randelem[2*num_x - 1] = 0; + //randelem[0] = 2*num_x - 1; + //std::shuffle(randelem.begin(), randelem.end(), std::default_random_engine(seed)); + mesh_ptr->GetGeckoElementOrdering(randelem); + mesh_ptr->ReorderElements(randelem); + + // strategy: + // 1) generate a fes for Lagrange elements of desired degree + // 2) create a Grid Function using a VectorFunctionCoefficient + // 4) use mesh_ptr->NewNodes(nodes, true) to set the mesh nodes + + // Problem: fes does not own fec, which is generated in this function's scope + // Solution: the grid function can own both the fec and fes + H1_FECollection *fec = new H1_FECollection(1, 2 /* = dim */); + FiniteElementSpace *fes = new FiniteElementSpace(mesh_ptr.get(), fec, 2, + Ordering::byVDIM); + + // Lambda function increases element density towards wall + double offset = m_offset; + double coeff = m_coeff; + // double s1 = m_s1; + // double s2 = m_s2; + + // double A = sqrt(s2)/sqrt(s1); + // double B = 1.0/sqrt(s1*s2); + + // // need to solve transcendental equation sinh(x)/x = B for x, use Newton's method + // double del = 5.0; + // double res = sinh(del) - B*del; + // double dres; + // while(res > 1e-15) + // { + // dres = cosh(del) - B; + // del = del - res/dres; + // res = sinh(del) - B*del; + // } + // double delta = del; + + auto xy_fun = [coeff, offset, num_y/*A, delta*/](const Vector& rt, Vector &xy) + { + xy(0) = rt(0) - 0.33333; + xy(1) = rt(1); + + double c = 1.0/num_y; + double b = log(offset)/(c - 1.0); + double a = 1.0/exp(1.0*b); + + // double u = 0.5*(1 + tanh(delta*(rt(1) - 0.5))/tanh(delta))); + + // double s = u/(A + (1-A)*u); + + + //Condense mesh near wall + xy(1) = (1 + tanh(coeff*(rt(1) - 1))/tanh(coeff)); + + // if(rt(1) > 0.0 && rt(1) < 1.0) + // { + // xy(1) = coeff*a*exp(b*rt(1)); + // //std::cout << xy(1) << std::endl; + // } + + ///TODO: condense mesh near wall transition as well + //double cx1 = 2.0/ + // if(rt(0) > 0.0 && rt(0) < 2.0) + // { + // xy(1) = coeff*a*exp(b*rt(0)); + // } + }; + VectorFunctionCoefficient xy_coeff(2, xy_fun); + GridFunction *xy = new GridFunction(fes); + xy->MakeOwner(fec); + xy->ProjectCoefficient(xy_coeff); + + mesh_ptr->NewNodes(*xy, true); + + // Assign extra boundary attribute + for (int i = 0; i < mesh_ptr->GetNBE(); ++i) + { + Element *elem = mesh_ptr->GetBdrElement(i); + + Array verts; + elem->GetVertices(verts); + + bool before = true; //before wall + for (int j = 0; j < 2; ++j) + { + auto vtx = mesh_ptr->GetVertex(verts[j]); + if (vtx[1] == 0.0 && vtx[0] <= 0.3333334) + { + before = before & true; + } + else + { + before = before & false; + } + } + if (before) + { + elem->SetAttribute(5); + } + } + mesh_ptr->bdr_attributes.Append(5); + // mesh_ptr->SetAttributes(); + // mesh_ptr->Finalize(); + + cout << "Number of bdr attr " << mesh_ptr->bdr_attributes.Size() <<'\n'; + cout << "Number of elements " << mesh_ptr->GetNE() <<'\n'; + + return mesh_ptr; +} diff --git a/sandbox/rans_walltest_options.json b/sandbox/rans_walltest_options.json new file mode 100644 index 00000000..8b1cbc8f --- /dev/null +++ b/sandbox/rans_walltest_options.json @@ -0,0 +1,77 @@ +{ + "flow-param": { + "mach": 0.2, + "aoa": 0.0, + "roll-axis": 0, + "pitch-axis": 1, + "chi": 3.0, + "Re": 5000000.0, + "Pr": 0.75, + "mu": 1.0, + "viscous-mms": false, + "rans-mms": false, + "sa-consts": [0.1355, 0.622, 0.666666666666667, 0.41, 0.3, 2, 7.1, 1.2, 0.5, 10, 16, 0.7, 0.9], + "sa-srcs": [1.0, 1.0] + }, + "space-dis": { + "degree": 1, + "lps-coeff": 1.0, + "basis-type": "csbp" + }, + "time-dis": { + "steady": true, + "steady-abstol": 1e-12, + "steady-restol": 1e-10, + "const-cfl": true, + "ode-solver": "PTC", + "t-final": 100.0, + "dt": 10.0, + "cfl": 0.5, + "res-exp": 2.0, + "start-up": true + }, + "nonlin-solver": { + "abstol": 1e-12, + "maxiter": 100, + "printlevel": 1, + "reltol": 1e-2, + "type": "inexact-newton" + }, + "lin-solver": { + "reltol": 1e-5, + "abstol": 1e-12, + "printlevel": 1, + "maxiter": 300 + }, + "lin-prec" : { + "type": "hypreilu", + "lev-fill": 3, + "ilu-type": 10, + "ilu-reorder": 1, + "printlevel": 0 + }, + "bcs": { + "far-field": [0, 1, 1, 1, 0], + "slip-wall": [0, 0, 0, 0, 1], + "no-slip-adiabatic": [1, 0, 0, 0, 0] + }, + "wall-func": { + "type": "true" + }, + "mesh": { + "num-x": 21, + "num-y": 12, + "offset": 1e-3, + "coeff": 1.0 + }, + "outputs": { + "drag": [1, 0, 0, 0, 0] + }, + "init-pert": 0.0, + "file-names": "rans_wturb_1e-3", + "try-rans": true, + "try-ns": false, + "compare": false + +} + diff --git a/src/common/default_options.cpp b/src/common/default_options.cpp index b87117a4..e1f730de 100644 --- a/src/common/default_options.cpp +++ b/src/common/default_options.cpp @@ -18,7 +18,10 @@ const nlohmann::json default_options {"pitch-axis", 1}, // axis in the "vertical" direction {"Re", 0.0}, // far-field Reynolds number {"Pr", 0.72}, // the Prandtl number - {"mu", -1.0} // nondimensional viscosity (if negative, use Sutherland's) + {"mu", -1.0}, // nondimensional viscosity (if negative, use Sutherland's) + {"sa-consts", {0.1355, 0.622, 0.666666666666667, 0.41, 0.3, 2, 7.1, 1.2, 0.5, 10, 16, 0.7, 0.9}}, + //Spalart-Allmaras turbulence model constants, defined in rans_fluxes.hpp + {"sa-srcs", {1.0, 1.0}} //for debug purposes, disable SA source terms }}, {"space-dis", // options related to spatial discretization diff --git a/src/common/inexact_newton.cpp b/src/common/inexact_newton.cpp index 8098973b..e5d4c0cc 100755 --- a/src/common/inexact_newton.cpp +++ b/src/common/inexact_newton.cpp @@ -1,5 +1,6 @@ #include "inexact_newton.hpp" - +#include +#include #include "utils.hpp" using namespace std; @@ -99,6 +100,9 @@ void InexactNewton::Mult(const Vector &b, Vector &x) const int it; double norm0, norm, norm_goal; const bool have_b = (b.Size() == Height()); + + Vector scaled_r(r.Size()); + Vector scaling_v(r.Size()); if (!iterative_mode) { x = 0.0; @@ -145,8 +149,36 @@ void InexactNewton::Mult(const Vector &b, Vector &x) const } jac = &oper->GetGradient(x); + + /// scaling stuff + // for(int k = 0; k < r.Size(); k++) + // { + // scaling_v(k) = 1/r(k); + // } + + //scaling_m = new SparseMatrix(scaling_v); + //jac_scaled = new ProductOperator(scaling_m, jac, false, false); + //HypreParMatrix* jac_scaled(jac); + //jac_scaled->ScaleRows(scaling_v); + + // if (it == 12) + // { + // stringstream nssolname; stringstream nssolname2; + // nssolname << "jac_" <PrintMatlab(matlab); + // //jac_scaled->PrintMatlab(matlab2); + // } + // stringstream nssolname; + // nssolname << "jac_" <PrintMatlab(matlab); // std::cout << "Get the jacobian matrix.\n"; + prec->SetOperator(*jac); + //scaling_m->Mult(r, scaled_r); //std::cout << "jac is set as one operator.\n"; prec->Mult(r, c); // c = [DF(x_i)]^{-1} [F(x_i)-b] // std::cout << "Solve for the newton step.\n"; diff --git a/src/common/inexact_newton.hpp b/src/common/inexact_newton.hpp index 2cb0027e..46a2d95c 100755 --- a/src/common/inexact_newton.hpp +++ b/src/common/inexact_newton.hpp @@ -45,6 +45,10 @@ class InexactNewton : public mfem::NewtonSolver protected: /// Jacobian of the nonlinear operator; needed by ComputeStepSize(); mutable Operator *jac; + /// Diagonal row scaling matrix + mutable SparseMatrix *scaling_m; + /// Jacobian scaled by an automatic matrix to improve condition number + mutable ProductOperator *jac_scaled; /// member vector saves the new x position. mutable mfem::Vector x_new; /// Parameters for inexact newton method. diff --git a/src/common/sbp_fe.cpp b/src/common/sbp_fe.cpp index 9d840e86..33a0f3c5 100644 --- a/src/common/sbp_fe.cpp +++ b/src/common/sbp_fe.cpp @@ -195,6 +195,16 @@ double SBPFiniteElement::getSymEntry(int di, int i, return Eij; } +void SBPFiniteElement::getSymEntryRevDiff(int di, int i, double Sij_bar, + mfem::DenseMatrix &adjJ_i_bar) const +{ + for (int k = 0; k < GetDim(); ++k) + { + //Eij += adjJ_i(k,di)*Q[k](i,i); + adjJ_i_bar(k,di) += Sij_bar*Q[k](i,i); + } +} + double SBPFiniteElement::getQEntry(int di, int i, int j, const mfem::DenseMatrix &adjJ_i, const mfem::DenseMatrix &adjJ_j) const @@ -209,6 +219,23 @@ double SBPFiniteElement::getQEntry(int di, int i, int j, } } +void SBPFiniteElement::getQEntryRevDiff(int di, int i, int j, double Sij_bar, + mfem::DenseMatrix &adjJ_i_bar, + mfem::DenseMatrix &adjJ_j_bar) const +{ + if (i == j) + { + getSymEntryRevDiff(di, i, Sij_bar, adjJ_i_bar); + //adjJ_j_bar *= 0.0; + } + else + { + getSkewEntryRevDiff(di, i, j, 0.5*Sij_bar, adjJ_i_bar, adjJ_j_bar); + // adjJ_i_bar *= 0.5; + // adjJ_j_bar *= 0.5; + } +} + void SBPFiniteElement::getProjOperator(DenseMatrix &P) const { MFEM_ASSERT( P.Size() == dof, ""); diff --git a/src/common/sbp_fe.hpp b/src/common/sbp_fe.hpp index a107cfbf..cb03bfe4 100644 --- a/src/common/sbp_fe.hpp +++ b/src/common/sbp_fe.hpp @@ -137,6 +137,9 @@ class SBPFiniteElement : public NodalFiniteElement /// \returns \f$ (E_{di})_{i,j} \f$ in physical space double getSymEntry(int di, int i, const mfem::DenseMatrix &adjJ_i) const; + void getSymEntryRevDiff(int di, int i, double Sij_bar, + mfem::DenseMatrix &adjJ_i_bar) const; + /// `(i,j)`th entry of skew-symmetric matrix \f$ S_{di} \f$ in physical space /// \param[in] di - desired physical space coordinate direction /// \param[in] i - row index for \f$ S_{di} \f$ @@ -163,6 +166,10 @@ class SBPFiniteElement : public NodalFiniteElement double getQEntry(int di, int i, int j, const mfem::DenseMatrix &adjJ_i, const mfem::DenseMatrix &adjJ_j) const; + void getQEntryRevDiff(int di, int i, int j, double Sij_bar, + mfem::DenseMatrix &adjJ_i_bar, + mfem::DenseMatrix &adjJ_j_bar) const; + /// Attempts to find the index corresponding to a given IntegrationPoint /// \param[in] ip - try to match the coordinates of this point /// \returns index - the index of the node corresponding to `ip` diff --git a/src/physics/fluidflow/CMakeLists.txt b/src/physics/fluidflow/CMakeLists.txt index 83e4d356..bd72ef0a 100644 --- a/src/physics/fluidflow/CMakeLists.txt +++ b/src/physics/fluidflow/CMakeLists.txt @@ -17,10 +17,16 @@ target_sources(mach navier_stokes_fluxes.hpp navier_stokes_integ_def.hpp navier_stokes_integ.hpp + navier_stokes_sens_integ.hpp navier_stokes.cpp navier_stokes.hpp viscous_integ_def.hpp viscous_integ.hpp + rans_fluxes.hpp + rans_integ.hpp + rans_integ_def.hpp + rans.cpp + rans.hpp ) target_include_directories(mach diff --git a/src/physics/fluidflow/euler.cpp b/src/physics/fluidflow/euler.cpp index 85e039a5..b90496de 100644 --- a/src/physics/fluidflow/euler.cpp +++ b/src/physics/fluidflow/euler.cpp @@ -285,12 +285,21 @@ double EulerSolver::calcStepSize(int iter, double t, { if (options["time-dis"]["steady"].template get()) { + double dt = 0.0; + // run a start-up period at constant dt until a threshold is reached + // if (start_up) + // { + // //dt = options["time-dis"]["dt"].template get(); + + // return dt_old; + // } + // ramp up time step for pseudo-transient continuation // TODO: the l2 norm of the weak residual is probably not ideal here // A better choice might be the l1 norm double res_norm = calcResidualNorm(state); double exponent = options["time-dis"]["res-exp"]; - double dt = options["time-dis"]["dt"].template get() * + dt = options["time-dis"]["dt"].template get() * pow(res_norm0 / res_norm, exponent); return max(dt, dt_old); } @@ -372,7 +381,7 @@ double EulerSolver::calcConservativeVarsL2Error( // This lambda function computes the error at a node // Beware: this is not particularly efficient, given the conditionals // Also **NOT thread safe!** - Vector qdiscrete(dim+2), qexact(dim+2); // define here to avoid reallocation + Vector qdiscrete(num_state), qexact(num_state); // define here to avoid reallocation auto node_error = [&](const Vector &discrete, const Vector &exact) -> double { if (entvar) @@ -389,7 +398,7 @@ double EulerSolver::calcConservativeVarsL2Error( double err = 0.0; if (entry < 0) { - for (int i = 0; i < dim+2; ++i) + for (int i = 0; i < num_state; ++i) { double dq = qdiscrete(i) - qexact(i); err += dq*dq; diff --git a/src/physics/fluidflow/euler.hpp b/src/physics/fluidflow/euler.hpp index 1de7841c..e3afe6b1 100644 --- a/src/physics/fluidflow/euler.hpp +++ b/src/physics/fluidflow/euler.hpp @@ -75,6 +75,8 @@ class EulerSolver : public AbstractSolver /// used to store the initial residual norm for PTC and convergence checks double res_norm0 = -1.0; + bool start_up = false; + /// Class constructor (protected to prevent misuse) /// \param[in] json_options - json object containing the options /// \param[in] smesh - if provided, defines the mesh for the problem diff --git a/src/physics/fluidflow/mesh_sens_integ.hpp b/src/physics/fluidflow/mesh_sens_integ.hpp index c6826871..b4dbf64c 100644 --- a/src/physics/fluidflow/mesh_sens_integ.hpp +++ b/src/physics/fluidflow/mesh_sens_integ.hpp @@ -67,6 +67,80 @@ class DyadicMeshSensIntegrator : public mfem::LinearFormIntegrator } }; +/// Integrator for mesh sensitivity of symmetrics viscous integrators +/// \tparam Derived - a class Derived from this one (needed for CRTP) +template +class SymmetricViscousMeshSensIntegrator : public mfem::LinearFormIntegrator +{ +public: + /// Constructs an integrator for dyadic-integrator mesh sensitivities + /// \param[in] state_vec - the state at which to evaluate the senstivity + /// \param[in] adjoint_vec - the adjoint that weights the residual + /// \param[in] num_state_vars - the number of state variables per node + /// \param[in] a - used to move residual to lhs (1.0) or rhs(-1.0) + SymmetricViscousMeshSensIntegrator(const mfem::GridFunction &state_vec, + const mfem::GridFunction &adjoint_vec, + int num_state_vars = 1, double a = 1.0) + : state(state_vec), adjoint(adjoint_vec), num_states(num_state_vars), + alpha(a) {} + + /// Construct the element local contribution to dF/dx + /// \param[in] el - the finite element whose dF/dx contribution we want + /// \param[in] trans - defines the reference to physical element mapping + /// \param[out] elvect - element local dF/dx + virtual void AssembleRHSElementVect(const mfem::FiniteElement &el, + mfem::ElementTransformation &trans, + mfem::Vector &elvect); + +protected: + /// The state vector used to evaluate fluxes + const mfem::GridFunction &state; + /// The adjoint vector that weights the residuals + const mfem::GridFunction &adjoint; + /// number of states (could extract from state or adjoint) + int num_states; + /// scales the terms; can be used to move to rhs/lhs + double alpha; +#ifndef MFEM_THREAD_SAFE + /// stores the result of calling the flux function + mfem::Vector fluxij; + /// used to store the adjugate of the mapping Jacobian at node i + mfem::DenseMatrix adjJ_i; + /// used to store the adjugate of the mapping Jacobian at node j + mfem::DenseMatrix adjJ_j; + /// used to store the adjugate of the mapping Jacobian at node j + mfem::DenseMatrix adjJ_k; + /// stores derivatives w.r.t. adjugate Jacobian at node i + mfem::DenseMatrix adjJ_i_bar; + /// stores derivatives w.r.t. adjugate Jacobian at node k + mfem::DenseMatrix adjJ_k_bar; + /// stores derivatives w.r.t. mesh nodes + mfem::DenseMatrix PointMat_bar; +#endif + + /// converts working variables to another set (e.g. conservative to entropy) + /// \param[in] u - working states that are to be converted + /// \param[out] w - transformed variables + /// \note This uses the CRTP, so it wraps a call to `convertVars` in Derived. + void convert(const mfem::Vector &u, mfem::Vector &w) + { + static_cast(this)->convertVars(u, w); + } + + /// applies symmetric matrices \f$ C_{d,:}(u) \f$ to input `Dw` + /// \param[in] d - index `d` in \f$ C_{d,:} \f$ matrices + /// \param[in] x - coordinate location at which scaling is evaluated + /// \param[in] u - state at which the symmetric matrices `C` are evaluated + /// \param[in] Dw - `Dw[:,d2]` stores derivative of `w` in direction `d2`. + /// \param[out] CDw - product of the multiplication between the `C` and `Dw`. + /// \note This uses the CRTP, so it wraps call to `applyScaling` in Derived. + void scale(int d, const mfem::Vector &x, const mfem::Vector &u, + const mfem::DenseMatrix &Dw, mfem::Vector &CDw) + { + static_cast(this)->applyScaling(d, x, u, Dw, CDw); + } +}; + /// Integrator for mesh sensitivity associated with boundary integrators /// \tparam Derived - a class Derived from this one (needed for CRTP) template diff --git a/src/physics/fluidflow/mesh_sens_integ_def.hpp b/src/physics/fluidflow/mesh_sens_integ_def.hpp index 78cbc7aa..4e0ce29b 100644 --- a/src/physics/fluidflow/mesh_sens_integ_def.hpp +++ b/src/physics/fluidflow/mesh_sens_integ_def.hpp @@ -78,6 +78,168 @@ void DyadicMeshSensIntegrator::AssembleRHSElementVect( } } + +template +void SymmetricViscousMeshSensIntegrator::AssembleRHSElementVect( + const mfem::FiniteElement &el, + mfem::ElementTransformation &trans, + mfem::Vector &elvect) +{ + using namespace mfem; + + // reverse-diff functions we need are only defined for IsoparametricTrans + IsoparametricTransformation &isotrans = + dynamic_cast(trans); + // extract the relevant sbp operator for this element + const FiniteElementSpace *fes = state.FESpace(); // Should check that fes match with adjoint + const FiniteElement *fe = fes->GetFE(trans.ElementNo); + const SBPFiniteElement &sbp = dynamic_cast(*fe); + // extract the state and adjoint values for this element + const IntegrationRule& ir = sbp.GetNodes(); + DenseMatrix u, psi, Dwi, Dpsi; + state.GetVectorValues(isotrans, ir, u); + adjoint.GetVectorValues(isotrans, ir, psi); + u.Transpose(); psi.Transpose(); + + int num_nodes = sbp.GetDof(); // number of state dofs + int ndof = el.GetDof(); // number of coord node dofs != num_nodes + int dim = el.GetDim(); + Vector xi, ui, uj, wi, wj, psii, psj, CDwi, CDpsi; +#ifdef MFEM_THREAD_SAFE + Vector fluxij; + DenseMatrix adjJ_i_bar, adjJ_k_bar, PointMat_bar; +#endif + fluxij.SetSize(num_states); + xi.SetSize(dim); + wj.SetSize(num_states); + ui.SetSize(num_states); + uj.SetSize(num_states); + psii.SetSize(num_states); + psj.SetSize(num_states); + Dwi.SetSize(num_states,dim); + Dpsi.SetSize(num_states,dim); + CDwi.SetSize(num_states); + CDpsi.SetSize(num_states); + adjJ_i.SetSize(dim); + adjJ_j.SetSize(dim); + adjJ_k.SetSize(dim); + adjJ_i_bar.SetSize(dim); + adjJ_k_bar.SetSize(dim); + elvect.SetSize(dim*ndof); + PointMat_bar.SetSize(dim, ndof); // PointMat_bar = dfdx + PointMat_bar = 0.0; + for (int i = 0; i < num_nodes; ++i) + { + // get the Jacobian (trans.Weight) and cubature weight (node.weight) + const IntegrationPoint &node = sbp.GetNodes().IntPoint(i); + isotrans.SetIntPoint(&node); + trans.SetIntPoint(&node); + double Hinv = 1.0 / (sbp.getDiagNormEntry(i) * trans.Weight()); + CalcAdjugate(trans.Jacobian(), adjJ_i); + u.GetRow(i, ui); + psi.GetRow(i, psii); + trans.Transform(node, xi); + + // compute Dwi and Dpsi derivatives at node i + Dwi = 0.0; Dpsi = 0.0; + for (int j = 0; j < num_nodes; ++j) + { + // Get mapping Jacobian adjugate and transform state to entropy vars + trans.SetIntPoint(&sbp.GetNodes().IntPoint(j)); + CalcAdjugate(trans.Jacobian(), adjJ_j); + u.GetRow(j, uj); psi.GetRow(j, psj); + convert(uj, wj); + for (int d = 0; d < dim; ++d) + { + double Qij = sbp.getQEntry(d, i, j, adjJ_i, adjJ_j); + for (int s = 0; s < num_states; ++s) + { + Dwi(s,d) += Qij * wj(s); + Dpsi(s,d) += Qij * psj(s); + } + } // loop over space dimensions d + } // loop over element nodes j + Dwi *= Hinv; + Dpsi *= Hinv; + double dHinv = -1.0 / (trans.Weight()); // multiply Dwi by this afterward + + // next, scale the derivatives (e.g. using \hat{C} matrices), and then + // apply adjugate jacobian sensitivities + adjJ_i_bar = 0.0; // prepare adjJ_i_bar for accumulation + for (int d = 0; d < dim; ++d) + { + scale(d, xi, ui, Dwi, CDwi); + scale(d, xi, ui, Dpsi, CDpsi); + + for (int k = 0; k < num_nodes; ++k) + { + adjJ_k_bar = 0.0; + trans.SetIntPoint(&sbp.GetNodes().IntPoint(k)); + CalcAdjugate(trans.Jacobian(), adjJ_k); + //double Qik = sbp.getQEntry(d, i, k, adjJ_i, adjJ_k); + // apply transposed derivative in `d` direction + // this evaluates Q_d'*(C_d,1 D_1 + C_d,2 D_2 + ...) w + // begin reverse sweep... + double Sij_u_bar = 0.0; + double Sij_psi_bar = 0.0; + for (int s = 0; s < num_states; ++s) + { + // res(k, s) += alpha * Qik * CDwi(s); + Sij_u_bar += alpha*CDwi(s)*psii(s); + Sij_psi_bar += alpha*CDpsi(s)*ui(s); + } + sbp.getQEntryRevDiff(d, i, k, Sij_u_bar, adjJ_i_bar, adjJ_k_bar); + sbp.getQEntryRevDiff(d, i, k, Sij_psi_bar, adjJ_i_bar, adjJ_k_bar); + isotrans.SetIntPoint(&ir.IntPoint(k)); + // adjJ_i = isotrans.AdjugateJacobian(); + isotrans.AdjugateJacobianRevDiff(adjJ_k_bar, PointMat_bar); + } + } // loop over space dimensions d + isotrans.SetIntPoint(&ir.IntPoint(i)); + // adjJ_i = isotrans.AdjugateJacobian(); + isotrans.AdjugateJacobianRevDiff(adjJ_i_bar, PointMat_bar); + + // Hinv sensitivity: compute dHinv, and accumulate energy in H_bar + Dwi *= dHinv; + double H_bar = 0.0; + for (int d = 0; d < dim; ++d) { + scale(d, xi, ui, Dwi, CDwi); + for (int k = 0; k < num_nodes; ++k) + { + psi.GetRow(k, psj); + // Get mapping Jacobian adjugate + trans.SetIntPoint(&sbp.GetNodes().IntPoint(k)); + CalcAdjugate(trans.Jacobian(), adjJ_k); + double Qik = sbp.getQEntry(d, i, k, adjJ_i, adjJ_k); + // apply transposed derivative in `d` direction + for (int s = 0; s < num_states; ++s) + { + H_bar += alpha * Qik * CDwi(s)*psj(s); + } + } // loop over element nodes k + } // loop over space dimensions d + DenseMatrix PointMat_bar_w(dim, ndof); + PointMat_bar_w = 0.0; + isotrans.WeightRevDiff(PointMat_bar_w); + + // add to total sensitivity + PointMat_bar.Add(H_bar, PointMat_bar_w); + } // i node loop + + + + // Insert PointMat_bar = dfdx into elvect + for (int i = 0; i < ndof; ++i) + { + for (int d = 0; d < dim; ++d) + { + elvect(d*ndof + i) = PointMat_bar(d, i); + } + } +} + + + template void BoundaryMeshSensIntegrator::AssembleRHSElementVect( const mfem::FiniteElement &el, mfem::ElementTransformation &trans, diff --git a/src/physics/fluidflow/navier_stokes_sens_integ.hpp b/src/physics/fluidflow/navier_stokes_sens_integ.hpp new file mode 100644 index 00000000..a6b8fcbc --- /dev/null +++ b/src/physics/fluidflow/navier_stokes_sens_integ.hpp @@ -0,0 +1,68 @@ +#ifndef MACH_NAVIER_STOKES_SENS_INTEG +#define MACH_NAVIER_STOKES_SENS_INTEG + +#include "adept.h" +#include "mfem.hpp" + +#include "mesh_sens_integ.hpp" +#include "navier_stokes_fluxes.hpp" + +using adept::adouble; + +namespace mach +{ + +/// Integrator for the mesh sensitivity of the Ismail-Roe domain integrator +/// \tparam dim - number of spatial dimensions (1, 2, or 3) +/// \tparam entvar - if true, the state variables are the entropy variables +/// \note This derived class uses the CRTP +template +class ESViscousMeshSensIntegrator : public SymmetricViscousMeshSensIntegrator< + ESViscousMeshSensIntegrator> +{ +public: + /// Construct an integrator for the Ismail-Roe flux over domains + /// \param[in] a - factor, usually used to move terms to rhs + ESViscousMeshSensIntegrator(const mfem::GridFunction &state_vec, + const mfem::GridFunction &adjoint_vec, + double Re_num, double Pr_num, + double vis = -1.0, double a = 1.0) + : SymmetricViscousMeshSensIntegrator>( + state_vec, adjoint_vec, dim+2, a), + Re(Re_num), Pr(Pr_num), mu(vis) {} + + void convertVars(const mfem::Vector &q, mfem::Vector &w) + { + calcEntropyVars(q.GetData(), w.GetData()); + } + + void applyScaling(int d, const mfem::Vector &x, const mfem::Vector &q, + const mfem::DenseMatrix &Dw, mfem::Vector &CDw) + { + double mu_Re = mu; + if (mu < 0.0) + { + mu_Re = calcSutherlandViscosity(q.GetData()); + } + mu_Re /= Re; + // applyViscousScaling(d, mu_Re, Pr, q.GetData(), Dw.GetData(), + // CDw.GetData()); + for(int i = 0; i < CDw.Size(); i++) + { + CDw(i) = mu_Re*Dw(i,d); + } + } +private: + /// Reynolds number + double Re; + /// Prandtl number + double Pr; + /// nondimensional dynamic viscosity + double mu; +}; + +#include "navier_stokes_sens_integ_def.hpp" + +} // namespace mach + +#endif \ No newline at end of file diff --git a/src/physics/fluidflow/navier_stokes_sens_integ_def.hpp b/src/physics/fluidflow/navier_stokes_sens_integ_def.hpp new file mode 100644 index 00000000..e69de29b diff --git a/src/physics/fluidflow/rans.cpp b/src/physics/fluidflow/rans.cpp new file mode 100644 index 00000000..2a31d90a --- /dev/null +++ b/src/physics/fluidflow/rans.cpp @@ -0,0 +1,413 @@ +#include +#include + +#include "rans.hpp" +#include "rans_integ.hpp" + +using namespace mfem; +using namespace std; + +namespace mach +{ + +template +RANavierStokesSolver::RANavierStokesSolver(const nlohmann::json &json_options, + unique_ptr smesh, + MPI_Comm comm) + : NavierStokesSolver(json_options, move(smesh), comm) +{ + if (entvar) + { + throw MachException("Entropy variables not implemented for RANS!"); + } + + // define free-stream parameters; may or may not be used, depending on case + chi_fs = this->options["flow-param"]["chi"].template get(); + mu = this->options["flow-param"]["mu"].template get(); + this->start_up = this->options["time-dis"]["start-up"].template get(); + vector sa = this->options["flow-param"]["sa-consts"].template get>(); + sacs.SetSize(13); + sacs = sa.data(); + getDistanceFunction(); +} + +template +void RANavierStokesSolver::addResVolumeIntegrators(double alpha) +{ + vector srcs = this->options["flow-param"]["sa-srcs"].template get>(); + Vector q_ref(dim+3); + this->getFreeStreamState(q_ref); + double d0 = getZeroDistance(); + + // MMS option ///NOTE: Doesn't work + bool mms = false;//this->options["flow-param"]["sa-mms"].template get(); + if (mms) + { + if (dim != 2) + { + throw MachException("SA MMS problem only available for 2D!"); + } + this->res->AddDomainIntegrator( + new SAMMSIntegrator(this->re_fs, this->pr_fs)); + this->res->AddDomainIntegrator(new SASourceIntegrator( + this->diff_stack, *dist, this->re_fs, sacs, mu, alpha, srcs[0], srcs[1], d0)); + } + else + { + + this->res->AddDomainIntegrator(new SAInviscidIntegrator( //Inviscid term + this->diff_stack, alpha)); + this->res->AddDomainIntegrator(new SAViscousIntegrator( //SAViscousIntegrator + this->diff_stack, this->re_fs, this->pr_fs, sacs, mu, alpha)); + // now add RANS integrators + this->res->AddDomainIntegrator(new SASourceIntegrator( + this->diff_stack, *dist, this->re_fs, sacs, mu, alpha, srcs[0], srcs[1], d0)); + // add LPS stabilization + double lps_coeff = this->options["space-dis"]["lps-coeff"].template get(); + this->res->AddDomainIntegrator(new SALPSIntegrator( + this->diff_stack, alpha, lps_coeff)); + } + +} + +template +void RANavierStokesSolver::addResBoundaryIntegrators(double alpha) +{ + auto &bcs = this->options["bcs"]; + double mu = this->options["flow-param"]["mu"].template get(); + + int idx = 0; + + if (bcs.find("slip-wall") != bcs.end()) + { // slip-wall boundary condition + vector tmp = bcs["slip-wall"].template get>(); + this->bndry_marker[idx].SetSize(tmp.size(), 0); + this->bndry_marker[idx].Assign(tmp.data()); + + Vector q_ref(dim+3); + this->getFreeStreamState(q_ref); + this->res->AddBdrFaceIntegrator( + new SAViscousSlipWallBC(this->diff_stack, this->fec.get(), + this->re_fs, this->pr_fs, sacs, mu, alpha), + this->bndry_marker[idx]); + idx++; + } + + if (bcs.find("no-slip-adiabatic") != bcs.end()) + { // no-slip-wall boundary condition + vector tmp = bcs["no-slip-adiabatic"].template get>(); + this->bndry_marker[idx].SetSize(tmp.size(), 0); + this->bndry_marker[idx].Assign(tmp.data()); + // reference state needed by penalty flux + Vector q_ref(dim+3); + this->getFreeStreamState(q_ref); + // Add the adiabatic flux BC + this->res->AddBdrFaceIntegrator( + new SANoSlipAdiabaticWallBC( + this->diff_stack, this->fec.get(), this->re_fs, + this->pr_fs, sacs, q_ref, mu, alpha), + this->bndry_marker[idx]); + idx++; + } +#if 0 + if (bcs.find("rans-inflow") != bcs.end()) + { + vector tmp = bcs["rans-inflow"].template get>(); + this->bndry_marker[idx].SetSize(tmp.size(), 0); + this->bndry_marker[idx].Assign(tmp.data()); + // get the in-flow state needed by the integrator + Vector q_in(dim+2); + getViscousInflowState(q_in); + // Add the rans-inflow integrator + this->res->AddBdrFaceIntegrator( + new RANSInflowBC(this->diff_stack, this->fec.get(), + re_fs, pr_fs, q_in, mu, alpha), + this->bndry_marker[idx]); + idx++; + } + if (bcs.find("rans-outflow") != bcs.end()) + { + vector tmp = bcs["rans-outflow"].template get>(); + this->bndry_marker[idx].SetSize(tmp.size(), 0); + this->bndry_marker[idx].Assign(tmp.data()); + // get the out-flow state needed by the integrator + Vector q_out(dim+2); + getRANSOutflowState(q_out); + // Add the rans-inflow integrator + this->res->AddBdrFaceIntegrator( + new RANSOutflowBC(this->diff_stack, this->fec.get(), + re_fs, pr_fs, q_out, mu, alpha), + this->bndry_marker[idx]); + idx++; + } +#endif + if (bcs.find("far-field") != bcs.end()) + { + // far-field boundary conditions + vector tmp = bcs["far-field"].template get>(); + mfem::Vector qfar(dim+3); + this->getFreeStreamState(qfar); + this->bndry_marker[idx].SetSize(tmp.size(), 0); + this->bndry_marker[idx].Assign(tmp.data()); + this->res->AddBdrFaceIntegrator( + new SAFarFieldBC(this->diff_stack, this->fec.get(), + this->re_fs, this->pr_fs, qfar, alpha), + this->bndry_marker[idx]); + idx++; + } +#if 0 + if (bcs.find("viscous-shock") != bcs.end()) + { + throw MachException("Viscous shock bc not implemented!"); + } +#endif + +} + +template +void RANavierStokesSolver::getRANSInflowState(Vector &q_in) +{ + vector tmp = this->options["flow-param"] + ["inflow-state"] + .template get>(); + if (tmp.size() != dim+3) + { + throw MachException("inflow-state option has wrong number of entries" + " for problem dimension!"); + } + for (int i = 0; i < dim+3; ++i) + { + q_in(i) = tmp[i]; + } +} + +template +void RANavierStokesSolver::getRANSOutflowState(Vector &q_out) +{ + vector tmp = this->options["flow-param"] + ["outflow-state"] + .template get>(); + if (tmp.size() != dim+3) + { + throw MachException("outflow-state option has wrong number of entries" + " for problem dimension!"); + } + for (int i = 0; i < dim+3; ++i) + { + q_out(i) = tmp[i]; + } +} + +template +void RANavierStokesSolver::getFreeStreamState(mfem::Vector &q_ref) +{ + q_ref = 0.0; + q_ref(0) = 1.0; + if (dim == 1) + { + q_ref(1) = q_ref(0)*this->mach_fs; // ignore angle of attack + } + else + { + q_ref(this->iroll+1) = q_ref(0)*this->mach_fs*cos(this->aoa_fs); + q_ref(this->ipitch+1) = q_ref(0)*this->mach_fs*sin(this->aoa_fs); + } + q_ref(dim+1) = 1/(euler::gamma*euler::gami) + 0.5*this->mach_fs*this->mach_fs; + q_ref(dim+2) = this->chi_fs*mu; +} + +static void pert(const Vector &x, Vector& p); + +template +void RANavierStokesSolver::iterationHook(int iter, + double t, double dt, + const ParGridFunction &state) +{ + string file = this->options["file-names"].template get(); + stringstream filename; + filename << file <<"_last"; + int num_rans = this->res->Height(); + Array disable_ns(num_rans); + //disable_ns = NULL; + //this->res->SetEssentialTrueDofs(disable_ns); + + this->checkJacobian(pert); + this->printSolution(filename.str(), 0); + + stringstream solname; + solname << file <<"_last.gf"; + ofstream ssol(solname.str()); ssol.precision(16); + this->u->Save(ssol); + + stringstream rsolname; + rsolname << file <<"_last_res.gf"; + ofstream rsol(rsolname.str()); rsol.precision(16); + GridFunction r(this->u->FESpace()); + this->res->Mult(*this->u, r); + r.Save(rsol); + double res_norm = this->calcResidualNorm(); + + cout << "Iter "<start_up) + // cout << "Start-Up Phase"< 1 && res_norm/this->res_norm0 > 1e-9 && sa_conv == false) + // { + // cout << "Navier Stokes Disabled"<start_up = true; + // for(int s = 0; s < num_rans; s++) + // { + // if((s+1)%(dim+3) == 0) + // disable_ns[s] = s-1; + // else + // disable_ns[s] = s; + // } + // } + // else if (iter <= 1) + // disable_ns = NULL; + // else + // { + // disable_ns = NULL; + // sa_conv = true; + // this->start_up = false; + // } + // this->res->SetEssentialTrueDofs(disable_ns); +} + + +template +void RANavierStokesSolver::terminalHook(int iter, double t_final, + const ParGridFunction &state) +{ + // double entropy = ent->GetEnergy(*u); + // entropylog << t_final << ' ' << entropy << endl; + // entropylog.close(); + + string file = this->options["file-names"].template get(); + stringstream filename; + filename << file <<"_final"; + stringstream solname; + solname << file <<"_final.gf"; + ofstream ssol(solname.str()); ssol.precision(30); + this->u->Save(ssol); +} + + +template +void RANavierStokesSolver::getDistanceFunction() +{ + std::string wall_type = + this->options["wall-func"]["type"].template get(); + int fe_order = this->options["space-dis"]["degree"].template get(); + H1_FECollection *dfec = new H1_FECollection(fe_order, dim); + FiniteElementSpace *dfes = new FiniteElementSpace(this->mesh.get(), dfec); + dist.reset(new GridFunction(dfes)); + + if (wall_type == "const") // uniform d value + { + double val = this->options["wall-func"]["val"].template get(); + ConstantCoefficient wall_coeff(val); + dist->ProjectCoefficient(wall_coeff); + } + if (wall_type == "y-dist") // y distance from the origin + { + auto walldist = [](const Vector &x) + { + // if(x(1) == 0.0) + // return 0.25*offset; //// not going to do it this way + // else + return x(1); + }; + FunctionCoefficient wall_coeff(walldist); + dist->ProjectCoefficient(wall_coeff); + } + if (wall_type == "true") // true computed distance function + { + // assemble wall attribute vector + auto &bcs = this->options["bcs"]; + vector tmp = bcs["no-slip-adiabatic"].template get>(); + vector tmp2 = bcs["slip-wall"].template get>(); + for(int i = 0; i < tmp.size(); i++) + { + tmp[i] = tmp[i] + tmp2[i]; + } + mfem::Array wall_marker; + wall_marker.SetSize(tmp.size(), 0); + wall_marker.Assign(tmp.data()); + + // create surface object + Surface surf(*this->mesh, wall_marker); + + // project surface function + auto walldist = [&surf](const Vector &x) + { + return surf.calcDistance(x); + }; + FunctionCoefficient wall_coeff(walldist); + dist->ProjectCoefficient(wall_coeff); + } + ///TODO: Differentiate true wall function +} + +template +double RANavierStokesSolver::getZeroDistance() +{ + std::string wall_type = + this->options["wall-func"]["type"].template get(); + // H1_FECollection *dfec = new H1_FECollection(1, dim); + // FiniteElementSpace *dfes = new FiniteElementSpace(this->mesh.get(), dfec); + // unique_ptr distcheck; + // distcheck.reset(new GridFunction(dfes)); + + double d0 = 0.0; + + if (wall_type == "const") // uniform d value + { + // doesn't make any difference + double val = this->options["wall-func"]["val"].template get(); + d0 = val; + } + else + { + // this should probably work for any distance function + double work = 1e15; + // exhaustive search of the mesh for the smallest d value + for(int i = 0; i < dist->Size(); i++) + { + if(dist->Elem(i) < work && dist->Elem(i) > 1e-15) + { + work = dist->Elem(i); + } + } + + // half the smallest computed distance from the wall + d0 = work/2; + } + ///TODO: Add option for proper wall distance function + + cout << "At-wall distance: "< normal_rand_ns(-1.0,1.0); + +// perturbation function used to check the jacobian in each iteration +void pert(const Vector &x, Vector& p) +{ + p.SetSize(5); + for (int i = 0; i < 5; i++) + { + p(i) = normal_rand_ns(gen_ns); + } +} + +// explicit instantiation +template class RANavierStokesSolver<1, false>; +template class RANavierStokesSolver<2, false>; +template class RANavierStokesSolver<3, false>; + +}//namespace mach + diff --git a/src/physics/fluidflow/rans.hpp b/src/physics/fluidflow/rans.hpp new file mode 100644 index 00000000..4dbe55f7 --- /dev/null +++ b/src/physics/fluidflow/rans.hpp @@ -0,0 +1,87 @@ +#ifndef MACH_RANS +#define MACH_RANS + +#include "mfem.hpp" + +#include "navier_stokes.hpp" +#include "surface.hpp" + +namespace mach +{ + +/// Solver for Reynolds-Averaged Navier-Stokes flows +/// dim - number of spatial dimensions (1, 2, or 3) +/// entvar - if true, the entropy variables are used in the integrators +template + +class RANavierStokesSolver : public NavierStokesSolver +{ +public: + /// Class constructor. + /// \param[in] opt_file_name - file where options are stored + /// \param[in] smesh - if provided, defines the mesh for the problem + /// \param[in] dim - number of dimensions + RANavierStokesSolver(const nlohmann::json &json_options, + std::unique_ptr smesh, + MPI_Comm comm); + + /// Sets `q_in` to the inflow conservative + SA variables + void getRANSInflowState(mfem::Vector &q_in); + + /// Sets `q_out` to the outflow conservative + SA variables + void getRANSOutflowState(mfem::Vector &q_out); + + /// Sets `q_ref` to the free-stream conservative + SA variables + void getFreeStreamState(mfem::Vector &q_ref); + + /// convert conservative variables to entropy variables + /// \param[in/out] state - the conservative/entropy variables + virtual void convertToEntvar(mfem::Vector &state) override + { + throw MachException("Entropy variables not implemented!!!"); + } + +protected: + + + /// Add volume/domain integrators to `res` based on `options` + /// \param[in] alpha - scales the data; used to move terms to rhs or lhs + /// \note This function calls NavierStokes::addVolumeIntegrators() first + virtual void addResVolumeIntegrators(double alpha); + + /// Add boundary-face integrators to `res` based on `options` + /// \param[in] alpha - scales the data; used to move terms to rhs or lhs + /// \note This function calls EulerSolver::addBoundaryIntegrators() first + virtual void addResBoundaryIntegrators(double alpha); + + /// Return the number of state variables + virtual int getNumState() override {return dim+3; } + + virtual void iterationHook(int iter, double t, double dt, const mfem::ParGridFunction &state) override; + + virtual void terminalHook(int iter, double t_final, const mfem::ParGridFunction &state) override; + + friend SolverPtr createSolver>( + const nlohmann::json &json_options, + std::unique_ptr smesh, + MPI_Comm comm); + + /// create mesh distance function from options + void getDistanceFunction(); + /// get the distance value to calculate at the wall itself + double getZeroDistance(); + /// SA constants vector + mfem::Vector sacs; + /// Gridfunction for wall distance function (unused for now) + std::unique_ptr dist; + /// free-stream SA viscosity ratio (nu_tilde/nu_material) + double chi_fs; + /// material dynamic viscosity + double mu; + + bool sa_conv = false; +}; + +} //namespace mach + +#endif \ No newline at end of file diff --git a/src/physics/fluidflow/rans_fluxes.hpp b/src/physics/fluidflow/rans_fluxes.hpp new file mode 100644 index 00000000..cabed2dd --- /dev/null +++ b/src/physics/fluidflow/rans_fluxes.hpp @@ -0,0 +1,751 @@ +/// Functions related to RANS/SA equations +#ifndef MACH_RANS_FLUXES +#define MACH_RANS_FLUXES + +#include // std::max + +#include "utils.hpp" +#include "navier_stokes_fluxes.hpp" + +/// Reference +/// sacs[0] = cb1 +/// sacs[1] = cb2 +/// sacs[2] = sigma +/// sacs[3] = kappa +/// sacs[4] = cw2 +/// sacs[5] = cw3 +/// sacs[6] = cv1 + +/// Not including trip terms for now + +/// sacs[7] = ct3 +/// sacs[8] = ct4 +/// sacs[9] = rlim +/// sacs[10] = cn1 + +///NOTE: These are related only to preventing negative St +/// sacs[11] = cv2 +/// sacs[12] = cv3 + +namespace mach +{ + +/// For constants related to the RANS/SA equations +namespace rans +{ + +} // namespace rans + +/// Returns the laminar suppression term in SA +/// \param[in] q - state used to define the laminar suppression +/// \param[in] mu - **nondimensionalized** viscosity +/// \param[in] sacs - Spalart-Allmaras constants +/// \returns ft2 - laminar suppression term +/// \tparam xdouble - typically `double` or `adept::adouble` +/// \tparam dim - number of spatial dimensions (1, 2, or 3) +template +xdouble calcSALaminarSuppression(const xdouble *q, const xdouble mu, + const xdouble *sacs) +{ + xdouble ct3 = sacs[7]; + xdouble ct4 = sacs[8]; + xdouble nu_tilde = q[dim+2];///q[0]; + xdouble nu_mat = mu/q[0]; + xdouble chi = nu_tilde/nu_mat; + xdouble ft2 = ct3*exp(-ct4*chi*chi); + return ft2; +} + +/// Returns the turbulent viscosity coefficient in SA +/// \param[in] q - state used to define the destruction +/// \param[in] mu - **nondimensionalized** viscosity +/// \param[in] sacs - Spalart-Allmaras constants +/// \returns fv1 - turbulent viscosity coefficient +/// \tparam xdouble - typically `double` or `adept::adouble` +/// \tparam dim - number of spatial dimensions (1, 2, or 3) +template +xdouble calcSACoefficient(const xdouble *q, const xdouble mu, + const xdouble *sacs) +{ + xdouble cv1 = sacs[6]; + xdouble nu_tilde = q[dim+2];///q[0]; + xdouble nu_mat = mu/q[0]; + xdouble chi = nu_tilde/nu_mat; + xdouble fv1 = chi*chi*chi/(cv1*cv1*cv1 + chi*chi*chi); + return fv1; +} + +/// Returns the production coefficient in SA +/// \param[in] q - state used to define the destruction +/// \param[in] mu - **nondimensionalized** viscosity +/// \param[in] sacs - Spalart-Allmaras constants +/// \returns fv2 - production coefficient +/// \tparam xdouble - typically `double` or `adept::adouble` +/// \tparam dim - number of spatial dimensions (1, 2, or 3) +template +xdouble calcSAProductionCoefficient(const xdouble *q, const xdouble mu, + const xdouble *sacs) +{ + xdouble nu_tilde = q[dim+2];///q[0]; + xdouble nu_mat = mu/q[0]; + xdouble chi = nu_tilde/nu_mat; + xdouble fv1 = calcSACoefficient(q, mu, sacs); + xdouble fv2 = 1.0 - chi/(1.0 + chi*fv1); + return fv2; +} + +/// Returns the modified vorticity in SA +/// \param[in] q - state used to define the destruction +/// \param[in] S - vorticity magnitude +/// \param[in] mu - **nondimensionalized** viscosity +/// \param[in] d - wall distance value +/// \param[in] Re - Reynold's number, if needed +/// \param[in] sacs - Spalart-Allmaras constants +/// \returns St - modified vorticity +/// \tparam xdouble - typically `double` or `adept::adouble` +/// \tparam dim - number of spatial dimensions (1, 2, or 3) +template +xdouble calcSAModifiedVorticity(const xdouble *q, const xdouble S, + const xdouble mu, const xdouble d, + const xdouble Re, const xdouble *sacs) +{ + xdouble nu_tilde = q[dim+2];///q[0]; + xdouble kappa = sacs[3]; + xdouble cv2 = sacs[11]; + xdouble cv3 = sacs[12]; + xdouble fv2 = calcSAProductionCoefficient(q, mu, sacs); + xdouble work = nu_tilde*fv2/(kappa*kappa*d*d); + xdouble St; + if (work/Re < -cv2*S) + St = S + (S*(cv2*cv2*S*Re + cv3*work))/((cv3-2*cv2)*S*Re - work); + else + St = S + work/Re; + return Re*St; +} + +/// Returns the destruction coefficient in SA +/// \param[in] q - state used to define the destruction +/// \param[in] mu - **nondimensionalized** viscosity +/// \param[in] d - wall distance value +/// \param[in] S - vorticity magnitude +/// \param[in] Re - Reynold's number, if needed +/// \param[in] sacs - Spalart-Allmaras constants +/// \returns fw - destruction coefficient +/// \tparam xdouble - typically `double` or `adept::adouble` +/// \tparam dim - number of spatial dimensions (1, 2, or 3) +template +xdouble calcSADestructionCoefficient(const xdouble *q, + const xdouble mu, const xdouble d, + const xdouble S, const xdouble Re, + const xdouble *sacs) +{ + using std::min; + using adept::min; + using std::pow; + using adept::pow; + xdouble nu_tilde = q[dim+2];///q[0]; + xdouble kappa = sacs[3]; + xdouble cw2 = sacs[4]; + xdouble cw3 = sacs[5]; + xdouble rlim = sacs[9]; + + xdouble St = calcSAModifiedVorticity(q, S, mu, d, Re, sacs); + //xdouble fv2 = calcSAProductionCoefficient(q, mu, sacs); + + xdouble work = nu_tilde/(St*kappa*kappa*d*d); + + xdouble r = min(work, rlim); + xdouble r6 = r*r*r*r*r*r; + xdouble g = r + cw2*(/*pow(r, 6.0)*/ r6 - r); + xdouble g6 = g*g*g*g*g*g; + xdouble cw36 = cw3*cw3*cw3*cw3*cw3*cw3; + xdouble work2 = (1.0 + /*pow(cw3, 6.0)*/ cw36)/ + (/*pow(g, 6.0)*/ g6 + /*pow(cw3, 6.0)*/ cw36); + xdouble fw = g*pow(work2, (1.0/6.0)); + return fw; +} + +/// Returns the production term in SA +/// \param[in] q - state used to define the production +/// \param[in] d - wall distance value +/// \param[in] S - vorticity magnitude (how to compute in h_grad space?) +/// \param[in] Re - Reynold's number, if needed +/// \param[in] sacs - Spalart-Allmaras constants +/// \returns P - production term +/// \tparam xdouble - typically `double` or `adept::adouble` +/// \tparam dim - number of spatial dimensions (1, 2, or 3) +template +xdouble calcSAProduction(const xdouble *q, + const xdouble mu, const xdouble d, + const xdouble S, const xdouble Re, + const xdouble *sacs) +{ + xdouble cb1 = sacs[0]; + xdouble nu_tilde = q[dim+2];///q[0]; + xdouble St = calcSAModifiedVorticity(q, S, mu, d, Re, sacs); + xdouble ft2 = calcSALaminarSuppression(q, mu, sacs); + xdouble P = cb1*(1.0-ft2)*St*nu_tilde; + return /*q[0]*/P/Re; +} + +/// Returns the destruction term in SA +/// \param[in] q - state used to define the destruction +/// \param[in] d - wall distance (precomputed from gridfunction) +/// \param[in] sacs - Spalart-Allmaras constants +/// \returns D - destruction term +/// \tparam xdouble - typically `double` or `adept::adouble` +/// \tparam dim - number of spatial dimensions (1, 2, or 3) +template +xdouble calcSADestruction(const xdouble *q, + const xdouble mu, const xdouble d, + const xdouble S, const xdouble Re, + const xdouble *sacs) +{ + xdouble cb1 = sacs[0]; + xdouble cb2 = sacs[1]; + xdouble sigma = sacs[2]; + xdouble kappa = sacs[3]; + xdouble chi_d = q[dim+2]/(d/*q[0]*/); + xdouble cw1 = cb1/(kappa*kappa) + (1+cb2)/sigma; + xdouble fw = calcSADestructionCoefficient(q, mu, d, S, Re, sacs); + xdouble ft2 = calcSALaminarSuppression(q, mu, sacs); + xdouble D = (cw1*fw - (cb1/(kappa*kappa))*ft2)*chi_d*chi_d; + return /*q[0]*/-D/Re; +} + +/// Returns the production term in negative SA +/// \param[in] q - state used to define the production +/// \param[in] S - vorticity magnitude (how to compute in h_grad space?) +/// \param[in] sacs - Spalart-Allmaras constants +/// \returns Pn - production term +/// \tparam xdouble - typically `double` or `adept::adouble` +/// \tparam dim - number of spatial dimensions (1, 2, or 3) +template +xdouble calcSANegativeProduction(const xdouble *q, const xdouble S, + const xdouble *sacs) +{ + xdouble cb1 = sacs[0]; + xdouble ct3 = sacs[7]; + xdouble nu_tilde = q[dim+2];///q[0]; + xdouble Pn = cb1*(1.0-ct3)*S*nu_tilde; + return /*q[0]*/Pn; +} + +/// Returns the destruction term in negative SA +/// \param[in] q - state used to define the destruction +/// \param[in] d - wall distance (precomputed from gridfunction) +/// \param[in] sacs - Spalart-Allmaras constants +/// \returns Dn - destruction term +/// \tparam xdouble - typically `double` or `adept::adouble` +/// \tparam dim - number of spatial dimensions (1, 2, or 3) +template +xdouble calcSANegativeDestruction(const xdouble *q, const xdouble d, + const xdouble Re, const xdouble *sacs) +{ + xdouble cb1 = sacs[0]; + xdouble cb2 = sacs[1]; + xdouble sigma = sacs[2]; + xdouble kappa = sacs[3]; + xdouble chi_d = q[dim+2]/(d/*q[0]*/); + xdouble cw1 = cb1/(kappa*kappa) + (1+cb2)/sigma; + xdouble Dn = -cw1*chi_d*chi_d; + return /*q[0]*/-Dn/Re; +} + +/// Returns the turbulent viscosity coefficient in negative SA +/// \param[in] q - state used to define the destruction +/// \param[in] mu - **nondimensionalized** viscosity +/// \param[in] sacs - Spalart-Allmaras constants +/// \returns fn - negative turbulent viscosity coefficient +/// \tparam xdouble - typically `double` or `adept::adouble` +/// \tparam dim - number of spatial dimensions (1, 2, or 3) +template +xdouble calcSANegativeCoefficient(const xdouble *q, const xdouble mu, + const xdouble *sacs) +{ + xdouble cn1 = sacs[10]; + xdouble nu_tilde = q[dim+2];///q[0]; + xdouble nu_mat = mu/q[0]; + xdouble chi = nu_tilde/nu_mat; + xdouble fn = (cn1 + chi*chi*chi)/(cn1 - chi*chi*chi); + if(q[dim+2] >= 0) + fn = 1.0; + return fn; +} + +/// Returns the "source" term in SA +/// \param[in] q - state used to define the destruction +/// \param[in] dir - turbulent viscosity gradient +/// \param[in] sacs - Spalart-Allmaras constants +/// \returns Sr - "source" term +/// \tparam xdouble - typically `double` or `adept::adouble` +/// \tparam dim - number of spatial dimensions (1, 2, or 3) +template +xdouble calcSASource(const xdouble *q, const xdouble *dir, + const xdouble *sacs) +{ + xdouble cb2 = sacs[1]; + xdouble sigma = sacs[2]; + xdouble Sr = (cb2/sigma)*dot(dir, dir); + return /*q[0]*/Sr; +} + +/// Returns the viscous "source" term in SA +/// \param[in] q - state used to define the destruction +/// \param[in] dir - turbulent viscosity gradient +/// \param[in] sacs - Spalart-Allmaras constants +/// \returns Sr - "source" term +/// \tparam xdouble - typically `double` or `adept::adouble` +/// \tparam dim - number of spatial dimensions (1, 2, or 3) +template +xdouble calcSASource2(const xdouble *q, const xdouble mu, const xdouble *dir, + const xdouble *dir2, const xdouble *sacs) +{ + xdouble nu_tilde = q[dim+2]/q[0]; + xdouble sigma = sacs[2]; + xdouble fn = calcSANegativeCoefficient(q, mu, sacs); + xdouble Sr = (1.0/sigma)*(mu/q[0] + fn*nu_tilde)*dot(dir, dir2); + return 0.0;//Sr; +} + + + +/// Returns the full sum of source terms +/// \param[in] q - state used to define the destruction +/// \param[in] mu - **nondimensionalized** viscosity +/// \param[in] d - wall distance value +/// \param[in] S - vorticity magnitude +/// \param[in] Re - Reynold's number, if needed +/// \param[in] d0 - at-wall distance value +/// \param[in] dir - turbulent viscosity gradient +/// \param[in] dir2 - density gradient +/// \param[in] sacs - Spalart-Allmaras constants +/// \returns src - full source term sum +/// \tparam xdouble - typically `double` or `adept::adouble` +/// \tparam dim - number of spatial dimensions (1, 2, or 3) +template +xdouble calcSAFullSource(const xdouble *q, + const xdouble mu, const xdouble d, + const xdouble S, const xdouble Re, + const xdouble d0, const xdouble *dir, + const xdouble *dir2, const xdouble *sacs, + const xdouble prod, const xdouble dest) +{ + xdouble src = 0.0; + src += calcSASource( + q, dir, sacs)/Re; + //src -= calcSASource2( + // q, mu, dir, dir2, sacs)/Re; + if (fabs(d) > 1e-12) + { + if (q[dim+2] < 0) + { + src += prod*calcSANegativeProduction( + q, S, sacs); + src += dest*calcSANegativeDestruction( + q, d, Re, sacs); + } + else + { + src += prod*calcSAProduction( + q, mu, d, S, Re, sacs); + src += dest*calcSADestruction( + q, mu, d, S, Re, sacs); + } + } + else + { + if (q[dim+2] < 0) + { + src += dest*calcSANegativeDestruction( + q, d0, Re, sacs); + } + else + { + src += dest*calcSADestruction( + q, mu, d0, S, Re, sacs); + } + } + return src; +} + + + +#if 0 +// Compute vorticity on an SBP element, needed for SA model terms +/// \param[in] q - the state over the element +/// \param[in] sbp - the sbp element whose shape functions we want +/// \param[in] Trans - defines the reference to physical element mapping +/// \param[out] curl - the curl of the velocity field at each node/int point +/// \tparam xdouble - typically `double` or `adept::adouble` +/// \tparam dim - number of spatial dimensions (1, 2, or 3) +template +void calcVorticitySBP(const mfem::DenseMatrix &q, const mfem::FiniteElement &fe, + mfem::ElementTransformation &Trans, mfem::DenseMatrix curl) +{ + using namespace mfem; + const SBPFiniteElement &sbp = dynamic_cast(fe); + int num_nodes = sbp.GetDof(); + DenseMatrix dq(q.Height(), q.Width()); //contains state derivatives in reference space + DenseMatrix dxi(dim*q.Height(), dim); //contains velocity derivatives in reference space + dq = 0.0; + for(int di = 0; di < dim; di++) + { + sbp.multWeakOperator(di, q, dq); + //need to scale with 1/H, probably when looping over nodes + dxi.CopyMN(dq, q.Height(), dim, 1, dim+1, di*q.Height(), 0); + } + + DenseMatrix dx(dim, dim); //contains velocity derivatives in absolute space + DenseMatrix dxin(dim, dim); //contains velocity derivatives in reference space for a node + Vector dxrow(dim); Vector curln(dim); + for (int i = 0; i < num_nodes; ++i) + { + // get the Jacobian + IntegrationPoint &node = sbp.GetNodes().IntPoint(i); + Trans.SetIntPoint(&node); + + // store nodal derivatives in a single matrix + for(int di = 0; di < dim; di++) + { + dxi.GetRow(di*q.Height() + i, dxrow); + dxin.SetRow(di, dxrow); + } + + // compute absolute derivatives + MultAtB(Trans.InverseJacobian(), dxin, dx); + + // compute curl at node and append + if(dim == 2) + { + curl(i, 0) = 0; + curl(i, 1) = 0; + curl(i, 2) = dx(1,0) - dx(0,1); + } + if(dim == 3) + { + curl(i, 0) = dx(2,1) - dx(1,2); + curl(i, 1) = dx(0,2) - dx(2,0); + curl(i, 2) = dx(1,0) - dx(0,1); + } + } +} + +// Compute gradient for the turbulence variable on an SBP element, +// needed for SA model terms +/// \param[in] q - the state over the element +/// \param[in] sbp - the sbp element whose shape functions we want +/// \param[in] Trans - defines the reference to physical element mapping +/// \param[out] grad - the gradient of the turbulence variable at each node +/// \tparam xdouble - typically `double` or `adept::adouble` +/// \tparam dim - number of spatial dimensions (1, 2, or 3) +template +void calcGradSBP(const mfem::DenseMatrix &q, const mfem::FiniteElement &fe, + mfem::ElementTransformation &Trans, mfem::DenseMatrix grad) +{ + using namespace mfem; + const SBPFiniteElement &sbp = dynamic_cast(fe); + int num_nodes = sbp.GetDof(); + DenseMatrix dq(q.Height(), q.Width()); //contains state derivatives in reference space + DenseMatrix dnu(q.Height(), dim); //contains turb variable derivatives in reference space + dq = 0.0; + for(int di = 0; di < dim; di++) + { + sbp.multWeakOperator(di, q, dq); + //need to scale with 1/H + dnu.SetCol(di, dq.GetColumn(dim+2)); + } + + Vector dnurow(dim); Vector gradn(dim); + for (int i = 0; i < num_nodes; ++i) + { + // get the Jacobian + IntegrationPoint &node = sbp.GetNodes().IntPoint(i); + Trans.SetIntPoint(&node); + + // store nodal grad in a vector + dnu.GetRow(i, dnurow); + + // compute absolute derivatives + Trans.InverseJacobian().MultTranspose(dnurow, gradn); + + // append result + grad.SetRow(i, gradn); + } +} +#endif + +// Compute vorticity at a point on an element, needed for SA model terms +/// \param[in] Dq - the state gradient Dq +/// \param[in] jac_inv - inverse jacobian at node +/// \param[out] curl - the curl of the velocity field at the node +/// \tparam xdouble - typically `double` or `adept::adouble` +/// \tparam dim - number of spatial dimensions (1, 2, or 3) +/// NOTE: Need inverse jacobian! Check EvalJacobian +template +void calcVorticity(const xdouble *Dq, const xdouble *jac_inv, + xdouble *curl) +{ + xdouble DqJ[dim*(dim+3)]; + + // compute absolute derivatives + for(int di = 0; di < dim+3; di++) + { + for(int d = 0; d < dim; d++) + { + xdouble work = 0; + for(int k = 0; k < dim; k++) + { + work += (Dq[di+k*(dim+3)])*jac_inv[k+d*dim]; + } + DqJ[di+d*(dim+3)] = work; + } + } + + // compute curl at node and append + if(dim == 2) + { + curl[0] = 0; + curl[1] = 0; + curl[2] = (DqJ[2] - DqJ[1 + dim+3]); + } + ////!!!!!!! + //curl[2] = 1.0;// min(abs(curl[2]), 0.0); + ////!!!!!!! + if(dim == 3) + { + curl[0] = (DqJ[3 + dim+3] - DqJ[2 + 2*(dim+3)]); + curl[1] = (DqJ[1 + 2*(dim+3)] - DqJ[3]); + curl[2] = (DqJ[2] - DqJ[1 + dim+3]); + } +} + +// Compute vorticity jacobian w.r.t. state on an SBP element, needed for SA model terms +/// \param[in] stack - adept stack +/// \param[in] Dq - the state gradient Dq +/// \param[in] jac_inv - inverse jacobian at node +/// \param[out] jac_curl - the curl jacobian +/// \tparam xdouble - typically `double` or `adept::adouble` +/// \tparam dim - number of spatial dimensions (1, 2, or 3) +template +void calcVorticityJacDw(adept::Stack &stack, const double *Dq, const double *jac_inv, + std::vector &jac_curl) +{ + // vector of active input variables + std::vector Dq_a(dim*(dim+3)); + std::vector jac_a(dim*(dim+3)); + // initialize adouble inputs + adept::set_values(Dq_a.data(), dim*(dim+3), Dq); + adept::set_values(jac_a.data(), dim*dim, jac_inv); + // start recording + stack.new_recording(); + // create vector of active output variables + std::vector curl_a(3); + // run algorithm + calcVorticity(Dq_a.data(), jac_a.data(), curl_a.data()); + // identify independent and dependent variables + stack.independent(Dq_a.data(), Dq_a.size()); + stack.dependent(curl_a.data(), curl_a.size()); + // compute and store jacobian in jac_inv + mfem::Vector work(dim*(dim+3)*3); + stack.jacobian_reverse(work.GetData()); + for (int i = 0; i < dim; ++i) + { + jac_curl[i] = (work.GetData() + i*(dim+3)*3); + } +} + +// Compute gradient for the turbulence variable on an element node, +// needed for SA model terms +/// \param[in] i - the state variable to take +/// \param[in] Dq - the state gradient Dq +/// \param[in] jac_inv - transformation jacobian at node +/// \param[out] grad - the gradient of the turbulence variable at each node +/// \tparam xdouble - typically `double` or `adept::adouble` +/// \tparam dim - number of spatial dimensions (1, 2, or 3) +template +void calcGrad(const int i, const xdouble *Dq, const xdouble *jac_inv, + xdouble *grad) +{ + xdouble DqJ[dim*(dim+3)]; + + // compute absolute derivatives + for(int di = 0; di < dim+3; di++) + { + for(int d = 0; d < dim; d++) + { + xdouble work = 0; + for(int k = 0; k < dim; k++) + { + work += Dq[di+k*(dim+3)]*jac_inv[k+d*dim]; + } + DqJ[di+d*(dim+3)] = work; + } + } + + for (int j = 0; j < dim; ++j) + { + grad[j] = DqJ[i + j*(dim+3)]; + } +} + +// Compute grad jacobian w.r.t. state on an SBP element, needed for SA model terms +/// \param[in] stack - adept stack +/// \param[in] i - the state variable to take +/// \param[in] Dq - the state gradient Dq +/// \param[in] jac_inv - transformation jacobian at node +/// \param[out] jac_grad - the gradient jacobian +/// \tparam xdouble - typically `double` or `adept::adouble` +/// \tparam dim - number of spatial dimensions (1, 2, or 3) +template +void calcGradJacDw(adept::Stack &stack, const int i, const double *Dq, const double *jac_inv, + std::vector &jac_grad) +{ + // vector of active input variables + std::vector Dq_a(dim*(dim+3)); + std::vector jac_a(dim*(dim+3)); + // initialize adouble inputs + adept::set_values(Dq_a.data(), dim*(dim+3), Dq); + adept::set_values(jac_a.data(), dim*dim, jac_inv); + // start recording + stack.new_recording(); + // create vector of active output variables + std::vector grad_a(dim); + // run algorithm + calcGrad(i, Dq_a.data(), jac_a.data(), grad_a.data()); + // identify independent and dependent variables + stack.independent(Dq_a.data(), Dq_a.size()); + stack.dependent(grad_a.data(), grad_a.size()); + // compute and store jacobian in jac_grad + mfem::Vector work(dim*(dim+3)*dim); + stack.jacobian_reverse(work.GetData()); + for (int j = 0; j < dim; ++j) + { + jac_grad[j] = (work.GetData() + j*(dim+3)*dim); + } +} + +/// Computes the dual-consistent term for the no-slip wall penalty (temporary SA version) +/// \param[in] dir - desired (scaled) normal vector to the wall +/// \param[in] mu - nondimensionalized dynamic viscosity +/// \param[in] mu2 - SA transport viscosity +/// \param[in] Pr - Prandtl number +/// \param[in] q - state at the wall location +/// \param[out] fluxes - fluxes to be scaled by Dw (column major) +/// \tparam xdouble - typically `double` or `adept::adouble` +/// \tparam dim - number of spatial dimensions (1, 2, or 3) +template +void calcNoSlipDualFluxSA(const xdouble *dir, xdouble mu, xdouble mu2, + double Pr, const xdouble *q, xdouble *fluxes) +{ + int num_state = dim + 3; + // zero out the fluxes, since applyCijMatrix accummulates + for (int i = 0; i < num_state*dim; ++i) + fluxes[i] = 0.0; + // evaluate the difference w - w_bc, where w_bc = [w[0], 0, 0, ...,w[dim+1], 0] + xdouble dw[num_state]; + dw[0] = 0.0; + dw[dim+1] = 0.0; + xdouble p = pressure(q); + for (int d = 0; d < dim; ++d) + { + dw[d+1] = q[d+1]/p; + } + dw[dim+2] = q[dim+2]; + // loop over the normal components + for (int d = 0; d < dim; ++d) + { + // loop over the derivative directions + for (int d2 = 0; d2 < dim; ++d2) + { + // we "sneak" dir[d] into the computation via mu + // !!!! we also sneak the sign into the computation here + applyCijMatrix(d2, d, -mu * dir[d], Pr, q, dw, + fluxes + (d2 * num_state)); + fluxes[dim+2 + d2*num_state] = -mu2 * dir[d] * dw[dim+2]; + } + } +#if 1 + // The following appears to be redundant (as far as the error goes) + // Project out the normal component from the last variable + xdouble nrm[dim]; + xdouble fac = 1.0 / sqrt(dot(dir, dir)); + xdouble flux_nrm = 0.0; + for (int d = 0; d < dim; ++d) + { + nrm[d] = dir[d] * fac; + flux_nrm += fluxes[d*num_state + dim+1]*nrm[d]; + } + for (int d = 0; d < dim; ++d) + fluxes[d*num_state + dim+1] -= flux_nrm*nrm[d]; +#endif +} + +/// Compute the product \f$ \sum_{d'=0}^{dim} C_{d,d'} D_{d'} w \f$ +/// Applies over conservative vars only +/// \param[in] d - desired space component of the flux +/// \param[in] mu - nondimensionalized dynamic viscosity +/// \param[in] Pr - Prandtl number +/// \param[in] q - state used to evaluate the \f$ C_d,d'} \f$ matrices +/// \param[in] Dw - derivatives of entropy varaibles, stored column-major +/// \param[out] mat_vec - stores the resulting matrix vector product +/// \tparam xdouble - typically `double` or `adept::adouble` +/// \tparam dim - number of spatial dimensions (1, 2, or 3) +template +void applyViscousScalingSA(int d, xdouble mu, double Pr, const xdouble *q, + const xdouble *Dw, xdouble *mat_vec) +{ + for (int k = 0; k < dim+3; ++k) + { + mat_vec[k] = 0.0; + } + for (int d2 = 0; d2 < dim; ++d2) { + applyCijMatrix(d, d2, mu, Pr, q, Dw+(d2*(dim+3)), mat_vec); + } +} + + +/// MMS source term for a particular Navier-Stokes verification +/// \param[in] mu - nondimensionalized dynamic viscosity +/// \param[in] Pr - Prandtl number +/// \param[in] x - location at which to evaluate the source +/// \param[out] src - the source value +/// \tparam xdouble - typically `double` or `adept::adouble` +template +void calcSAMMS(double mu, double Pr, const xdouble *x, xdouble *src) +{ + double gamma = euler::gamma; + const double rho0 = 1.0; + const double rhop = 0.05; + const double U0 = 0.5; + const double Up = 0.05; + const double T0 = 1.0; + const double Tp = 0.05; + double kappa = mu * gamma / (Pr * euler::gami); + const double chif = 100.0; + // const double nut = chif*mu*x[1]; + // const double rho = rho0 + rhop*sin(M_PI*x[0])*sin(M_PI*x[0])*sin(M_PI*x[1]); + // const double chi = nut*rho/mu; + // const double fv1 = (chi*chi*chi)/(chi*chi*chi + 7.1*7.1*7.1); + // const double fv2 = 1 - chi/(1 + chi*fv1); + // const double mu_t = mu + rho*nut*fv1; + // const double dummy = sqrt(pow(-4*U0*x[1] + 4*U0*(1 - x[1]) + 2*M_PI*Up*pow(sin(M_PI*x[0]), 2)*cos(2*M_PI*x[1]) + 4*M_PI*Up*sin(2*M_PI*x[0])*sin(M_PI*x[1])*cos(2*M_PI*x[0]), 2)); + + src[0] = 0 ; + + src[1] = 0 ; + + src[2] = 0 ; + + src[3] = 0 ; + + //-P + //src[4] = chif*mu*x[1]*(0.13550000000000001 - 0.16259999999999999*exp(-0.5*pow(chif, 2)*pow(x[1], 2)*pow(rho0 + rhop*pow(sin(M_PI*x[0]), 2)*sin(M_PI*x[1]), 2)))*(-3.5693039857227848*chif*mu*(chif*x[1]*(rho0 + rhop*pow(sin(M_PI*x[0]), 2)*sin(M_PI*x[1]))/(pow(chif, 4)*pow(x[1], 4)*pow(rho0 + rhop*pow(sin(M_PI*x[0]), 2)*sin(M_PI*x[1]), 4)/(pow(chif, 3)*pow(x[1], 3)*pow(rho0 + rhop*pow(sin(M_PI*x[0]), 2)*sin(M_PI*x[1]), 3) + 357.91099999999994) + 1) - 1)/x[1] - 0.67333333333333334*sqrt(pow(-4*U0*x[1] + 4*U0*(1 - x[1]) + 2*M_PI*Up*pow(sin(M_PI*x[0]), 2)*cos(2*M_PI*x[1]) + 4*M_PI*Up*sin(2*M_PI*x[0])*sin(M_PI*x[1])*cos(2*M_PI*x[0]), 2))) ; + + //-P-src + src[4] = -1.0718113612004292*pow(chif, 2)*pow(mu, 2) - chif*mu*x[1]*(0.13550000000000001 - 0.16259999999999999*exp(-0.5*pow(chif, 2)*pow(x[1], 2)*pow(rho0 + rhop*pow(sin(M_PI*x[0]), 2)*sin(M_PI*x[1]), 2)))*(3.5693039857227848*chif*mu*(chif*x[1]*(rho0 + rhop*pow(sin(M_PI*x[0]), 2)*sin(M_PI*x[1]))/(pow(chif, 4)*pow(x[1], 4)*pow(rho0 + rhop*pow(sin(M_PI*x[0]), 2)*sin(M_PI*x[1]), 4)/(pow(chif, 3)*pow(x[1], 3)*pow(rho0 + rhop*pow(sin(M_PI*x[0]), 2)*sin(M_PI*x[1]), 3) + 357.91099999999994) + 1) - 1)/x[1] + 0.67333333333333334*sqrt(pow(-4*U0*x[1] + 4*U0*(1 - x[1]) + 2*M_PI*Up*pow(sin(M_PI*x[0]), 2)*cos(2*M_PI*x[1]) + 4*M_PI*Up*sin(2*M_PI*x[0])*sin(M_PI*x[1])*cos(2*M_PI*x[0]), 2))) ;} + + +} // namespace mach + +#endif diff --git a/src/physics/fluidflow/rans_integ.hpp b/src/physics/fluidflow/rans_integ.hpp new file mode 100644 index 00000000..deae219d --- /dev/null +++ b/src/physics/fluidflow/rans_integ.hpp @@ -0,0 +1,977 @@ +#ifndef MACH_RANS_INTEG +#define MACH_RANS_INTEG + +#include "adept.h" +#include "mfem.hpp" + +#include "navier_stokes_integ.hpp" +#include "euler_fluxes.hpp" +#include "navier_stokes_fluxes.hpp" +#include "rans_fluxes.hpp" + +using adept::adouble; +using namespace std; + +namespace mach +{ + +/// Volume integrator for inviscid terms, including SA variable +/// \tparam dim - number of spatial dimensions (1, 2, or 3) +/// \note This derived class uses the CRTP +template +class SAInviscidIntegrator : public DyadicFluxIntegrator< + SAInviscidIntegrator> +{ +public: + /// Construct an inviscid integrator with SA terms + /// \param[in] diff_stack - for algorithmic differentiation + /// \param[in] a - used to move residual to lhs (1.0) or rhs(-1.0) + SAInviscidIntegrator(adept::Stack &diff_stack, double a = 1.0) + : DyadicFluxIntegrator>( + diff_stack, dim+3, a) {} + + /// Ismail-Roe two-point (dyadic) flux function with additional variable + /// \param[in] di - physical coordinate direction in which flux is wanted + /// \param[in] qL - state variables at "left" state + /// \param[in] qR - state variables at "right" state + /// \param[out] flux - fluxes in the direction `di` + /// \note This is simply a wrapper for the function in `euler_fluxes.hpp` + void calcFlux(int di, const mfem::Vector &qL, + const mfem::Vector &qR, mfem::Vector &flux); + + /// Compute the Jacobians of `flux` with respect to `u_left` and `u_right` + /// \param[in] di - desired coordinate direction for flux + /// \param[in] qL - the "left" state + /// \param[in] qR - the "right" state + /// \param[out] jacL - Jacobian of `flux` w.r.t. `qL` + /// \param[out] jacR - Jacobian of `flux` w.r.t. `qR` + void calcFluxJacStates(int di, const mfem::Vector &qL, + const mfem::Vector &qR, + mfem::DenseMatrix &jacL, + mfem::DenseMatrix &jacR); +private: + mfem::Vector qfs; +}; + +/// Integrator for RANS SA Production, Destruction, Viscous-Like terms +/// \tparam dim - number of spatial dimensions (1, 2, or 3) +/// \note This derived class uses the CRTP +template +class SASourceIntegrator : public mfem::NonlinearFormIntegrator +{ +public: + /// Construct an integrator for RANS SA Production term + /// \param[in] diff_stack - for algorithmic differentiation + /// \param[in] distance - wall distance function projected onto a space + /// \param[in] re_fs - freestream reynolds number + /// \param[in] sa_params - Spalart-Allmaras model constants + /// \param[in] vis - nondimensional dynamic viscosity (use Sutherland if neg) + /// \param[in] a - used to move residual to lhs (1.0) or rhs(-1.0) + /// \param[in] P, D - control production and destruction terms for debugging + SASourceIntegrator(adept::Stack &diff_stack, mfem::GridFunction distance, double re_fs, + mfem::Vector sa_params, double vis = -1.0, double a = 1.0, double P = 1.0, double D = 1.0, double dmin = 1e-4) + : alpha(a), d0(dmin), prod(P), dest(D), mu(vis), stack(diff_stack), num_states(dim+3), Re(re_fs), sacs(sa_params), + dist(distance) {} + + /// Construct the element local residual + /// \param[in] fe - the finite element whose residual we want + /// \param[in] Trans - defines the reference to physical element mapping + /// \param[in] elfun - element local state function + /// \param[out] elvect - element local residual + virtual void AssembleElementVector(const mfem::FiniteElement &fe, + mfem::ElementTransformation &Trans, + const mfem::Vector &elfun, + mfem::Vector &elvect); + + /// Construct the element local Jacobian + /// \param[in] fe - the finite element whose Jacobian we want + /// \param[in] Trans - defines the reference to physical element mapping + /// \param[in] elfun - element local state function + /// \param[out] elmat - element local Jacobian + virtual void AssembleElementGrad(const mfem::FiniteElement &fe, + mfem::ElementTransformation &Trans, + const mfem::Vector &elfun, + mfem::DenseMatrix &elmat); + + +private: + +protected: + /// nondimensional dynamic viscosity + double mu; + /// number of states + int num_states; + /// scales the terms; can be used to move to rhs/lhs + double alpha; + /// Freestream Reynolds number + double Re; + /// 0-wall distance to evaluate + double d0; + /// vector of SA model parameters + mfem::Vector sacs; + /// stack used for algorithmic differentiation + adept::Stack &stack; + /// activate production and destruction terms + double prod; double dest; +#ifndef MFEM_THREAD_SAFE + /// the coordinates of node i + mfem::Vector xi; + /// used to reference the states at node i + mfem::Vector ui; + /// used to reference the conservative states at node i + mfem::Vector uci; + /// used to reference the gradient of nu at node i + mfem::Vector grad_nu_i; + /// used to reference the gradient of rho at node i + mfem::Vector grad_rho_i; + /// used to reference the curl at node i + mfem::Vector curl_i; + /// stores the result of calling the flux function + mfem::Vector fluxi; + /// used to store the gradient on the element + mfem::DenseMatrix grad; + /// used to store the gradient on the element + mfem::DenseMatrix curl; + /// used to store the flux Jacobian at node i + mfem::DenseMatrix flux_jaci; + /// used to store the flux at each node + mfem::DenseMatrix elflux; + /// used to store the residual in (num_states, Dof) format + mfem::DenseMatrix elres; + /// wall distance function + mfem::GridFunction dist; +#endif +}; + +/// Integrator for SA no-slip adiabatic-wall boundary condition +/// \tparam dim - number of spatial dimensions (1, 2, or 3) +/// \note This derived class uses the CRTP +template +class SANoSlipAdiabaticWallBC : public ViscousBoundaryIntegrator> +{ +public: + /// Constructs an integrator for a no-slip, adiabatic boundary flux + /// \param[in] diff_stack - for algorithmic differentiation + /// \param[in] fe_coll - used to determine the face elements + /// \param[in] Re_num - Reynolds number + /// \param[in] Pr_num - Prandtl number + /// \param[in] q_ref - a reference state (needed by penalty) + /// \param[in] vis - viscosity (if negative use Sutherland's law) + /// \param[in] a - used to move residual to lhs (1.0) or rhs(-1.0) + SANoSlipAdiabaticWallBC(adept::Stack &diff_stack, + const mfem::FiniteElementCollection *fe_coll, + double Re_num, double Pr_num, mfem::Vector sa_params, + const mfem::Vector &q_ref, double vis = -1.0, + double a = 1.0) + : ViscousBoundaryIntegrator>( + diff_stack, fe_coll, dim + 3, a), + Re(Re_num), Pr(Pr_num), sacs(sa_params), + qfs(q_ref), mu(vis), work_vec(dim + 3) {} + + /// converts conservative variables to entropy variables + /// \param[in] q - conservative variables that are to be converted + /// \param[out] w - entropy variables corresponding to `q` + /// \note a wrapper for the relevant function in `euler_fluxes.hpp` + void convertVars(const mfem::Vector &q, mfem::Vector &w) + { + calcEntropyVars(q.GetData(), w.GetData()); + w(dim+2) = q(dim+2); + } + + /// Compute the Jacobian of the mapping `convert` w.r.t. `u' + /// \param[in] q - conservative variables that are to be converted + /// \param[out] dwdu - Jacobian of entropy variables w.r.t. `u` + void convertVarsJacState(const mfem::Vector &q, mfem::DenseMatrix &dwdu) + { + dwdu = 0.0; + convertVarsJac(q, this->stack, dwdu); + dwdu(dim+2,dim+2) = 1.0; + } + + /// Compute entropy-stable, no-slip, adiabatic-wall boundary flux + /// \param[in] x - coordinate location at which flux is evaluated (not used) + /// \param[in] dir - vector normal to the boundary at `x` + /// \param[in] jac - mapping Jacobian (needed by no-slip penalty) + /// \param[in] q - conservative variables at which to evaluate the flux + /// \param[in] Dw - space derivatives of the entropy variables + /// \param[out] flux_vec - value of the flux + void calcFlux(const mfem::Vector &x, const mfem::Vector &dir, double jac, + const mfem::Vector &q, const mfem::DenseMatrix &Dw, + mfem::Vector &flux_vec); + + /// Compute Jacobian of entropy-stable, no-slip, adiabatic-wall boundary flux + /// w.r.t states + /// \param[in] x - coordinate location at which flux is evaluated (not used) + /// \param[in] dir - vector normal to the boundary at `x` + /// \param[in] jac - mapping Jacobian (needed by no-slip penalty) + /// \param[in] q - conservative variables at which to evaluate the flux + /// \param[in] Dw - space derivatives of the entropy variables + /// \param[out] flux_jac - jacobian of the flux w.r.t. states + void calcFluxJacState(const mfem::Vector &x, const mfem::Vector &dir, double jac, + const mfem::Vector &q, const mfem::DenseMatrix &Dw, + mfem::DenseMatrix &flux_jac); + + /// Compute Jacobian of entropy-stable, no-slip, adiabatic-wall boundary flux + /// w.r.t vector of entropy-variables' derivatives + /// \param[in] x - coordinate location at which flux is evaluated (not used) + /// \param[in] dir - vector normal to the boundary at `x` + /// \param[in] jac - mapping Jacobian (needed by no-slip penalty) + /// \param[in] q - conservative variables at which to evaluate the flux + /// \param[in] Dw - space derivatives of the entropy variables + /// \param[out] flux_jac - jacobian of the flux w.r.t. entropy-variables' derivatives + void calcFluxJacDw(const mfem::Vector &x, const mfem::Vector &dir, double jac, + const mfem::Vector &q, const mfem::DenseMatrix &Dw, + std::vector &flux_jac); + + /// Compute flux terms that are multiplied by test-function derivative + /// \param[in] x - coordinate location at which flux is evaluated (not used) + /// \param[in] dir - vector normal to the boundary at `x` + /// \param[in] q - conservative variables at which to evaluate the fluxes + /// \param[out] flux_mat[:,di] - to be scaled by derivative `D_[di] v` + void calcFluxDv(const mfem::Vector &x, const mfem::Vector &dir, + const mfem::Vector &q, mfem::DenseMatrix &flux_mat); + // { + // flux_mat = 0.0; + // } + + /// Compute the Jacobian of calcFluxDv w.r.t. state + /// \param[in] x - coordinate location at which fluxes are evaluated + /// \param[in] dir - vector normal to the boundary at `x` + /// \param[in] q - conservative variables at which to evaluate the flux + /// \param[in] flux_jac[di] - Jacobian of calcFluxDv[di] with respect to `q` + void calcFluxDvJacState(const mfem::Vector &x, const mfem::Vector dir, + const mfem::Vector &u, + std::vector &flux_jac); + // { + // for (int d = 0; d < dim; ++d) + // flux_jac[d] = 0.0; + // } + + /// Computes boundary node contribution to the surface force + /// \param[in] x - coordinate location at which function is evaluated + /// \param[in] dir - vector normal to the boundary at `x` + /// \param[in] jac - mapping Jacobian determinant (needed by some fluxes) + /// \param[in] u - state at which to evaluate the function + /// \param[in] Dw - `Dw[:,di]` is the derivative of `w` in direction `di` + /// \returns fun - stress at given point + double calcBndryFun(const mfem::Vector &x, const mfem::Vector &dir, + double jac, const mfem::Vector &u, + const mfem::DenseMatrix &Dw) {return 0.0;} + +private: + /// Reynolds number + double Re; + /// Prandtl number + double Pr; + /// nondimensionalized dynamic viscosity + double mu; + /// Fixed state used to compute no-slip penalty matrix + mfem::Vector qfs; + /// work space for flux computations + mfem::Vector work_vec; + /// vector of SA model parameters + mfem::Vector sacs; +}; + +/// Integrator for SA viscous slip-wall boundary condition +/// \tparam dim - number of spatial dimensions (1, 2, or 3) +/// \note This derived class uses the CRTP +/// \note This is the same as the inviscid slip wall, but it provides the +/// necessary entropy-variable gradient flux. +template +class SAViscousSlipWallBC : public ViscousBoundaryIntegrator> +{ +public: + /// Constructs an integrator for a viscous slip-wall boundary + /// \param[in] diff_stack - for algorithmic differentiation + /// \param[in] fe_coll - used to determine the face elements + /// \param[in] Re_num - Reynolds number + /// \param[in] Pr_num - Prandtl number + /// \param[in] vis - viscosity (if negative use Sutherland's law) + /// \param[in] a - used to move residual to lhs (1.0) or rhs(-1.0) + SAViscousSlipWallBC(adept::Stack &diff_stack, + const mfem::FiniteElementCollection *fe_coll, + double Re_num, double Pr_num, mfem::Vector sa_params, + double vis = -1.0, double a = 1.0) + : ViscousBoundaryIntegrator>( + diff_stack, fe_coll, dim + 3, a), + Re(Re_num), Pr(Pr_num), mu(vis), + sacs(sa_params), work_vec(dim + 3) {} + + /// converts conservative variables to entropy variables + /// \param[in] q - conservative variables that are to be converted + /// \param[out] w - entropy variables corresponding to `q` + /// \note a wrapper for the relevant function in `euler_fluxes.hpp` + void convertVars(const mfem::Vector &q, mfem::Vector &w) + { + calcEntropyVars(q.GetData(), w.GetData()); + w(dim+2) = q(dim+2); + } + + /// Compute the Jacobian of the mapping `convert` w.r.t. `u` + /// \param[in] q - conservative variables that are to be converted + /// \param[out] dwdu - Jacobian of entropy variables w.r.t. `u` + void convertVarsJacState(const mfem::Vector &q, mfem::DenseMatrix &dwdu) + { + dwdu = 0.0; + convertVarsJac(q, this->stack, dwdu); + dwdu(dim+2,dim+2) = 1.0; + } + /// Compute flux corresponding to a viscous slip-wall boundary + /// \param[in] x - coordinate location at which flux is evaluated (not used) + /// \param[in] dir - vector normal to the boundary at `x` + /// \param[in] jac - mapping Jacobian (needed by no-slip penalty) + /// \param[in] q - conservative variables at which to evaluate the flux + /// \param[in] Dw - space derivatives of the entropy variables + /// \param[out] flux_vec - value of the flux + void calcFlux(const mfem::Vector &x, const mfem::Vector &dir, double jac, + const mfem::Vector &q, const mfem::DenseMatrix &Dw, + mfem::Vector &flux_vec); + + /// Compute jacobian of flux corresponding to a viscous slip-wall boundary + /// w.r.t `states` + /// \param[in] x - coordinate location at which flux is evaluated + /// \param[in] dir - vector normal to the boundary at `x` + /// \param[in] jac - mapping Jacobian (needed by no-slip penalty) + /// \param[in] q - conservative variables at which to evaluate the flux + /// \param[in] Dw - space derivatives of the entropy variables (not used yet) + /// \param[out] flux_jac - jacobian of the flux + void calcFluxJacState(const mfem::Vector &x, const mfem::Vector &dir, + double jac, const mfem::Vector &q, + const mfem::DenseMatrix &Dw, + mfem::DenseMatrix &flux_jac); + + /// Compute jacobian of flux corresponding to a viscous slip-wall boundary + /// w.r.t `entrpy-variables' derivatives` + /// \param[in] x - coordinate location at which flux is evaluated + /// \param[in] dir - vector normal to the boundary at `x` + /// \param[in] jac - mapping Jacobian (needed by no-slip penalty) + /// \param[in] q - conservative variables at which to evaluate the flux + /// \param[in] Dw - space derivatives of the entropy variables + /// \param[out] flux_jac - jacobian of the flux + void calcFluxJacDw(const mfem::Vector &x, const mfem::Vector &dir, + double jac, const mfem::Vector &q, + const mfem::DenseMatrix &Dw, + std::vector &flux_jac); + + /// Compute flux terms that are multiplied by test-function derivative + /// \param[in] x - coordinate location at which flux is evaluated (not used) + /// \param[in] dir - vector normal to the boundary at `x` + /// \param[in] q - conservative variables at which to evaluate the fluxes + /// \param[out] flux_mat[:,di] - to be scaled by derivative `D_[di] v` + void calcFluxDv(const mfem::Vector &x, const mfem::Vector &dir, + const mfem::Vector &q, mfem::DenseMatrix &flux_mat) + { + flux_mat = 0.0; + } + + /// Compute the Jacobian of calcFluxDv w.r.t. state + /// \param[in] x - coordinate location at which fluxes are evaluated + /// \param[in] dir - vector normal to the boundary at `x` + /// \param[in] q - conservative variables at which to evaluate the flux + /// \param[in] flux_jac[di] - Jacobian of calcFluxDv[di] with respect to `q` + void calcFluxDvJacState(const mfem::Vector &x, const mfem::Vector dir, + const mfem::Vector &u, + std::vector &flux_jac) + { + for (int d = 0; d < dim; ++d) + flux_jac[d] = 0.0; + } + + /// Computes boundary node contribution to the surface force + /// \param[in] x - coordinate location at which function is evaluated + /// \param[in] dir - vector normal to the boundary at `x` + /// \param[in] jac - mapping Jacobian determinant (needed by some fluxes) + /// \param[in] u - state at which to evaluate the function + /// \param[in] Dw - `Dw[:,di]` is the derivative of `w` in direction `di` + /// \returns fun - stress at given point + double calcBndryFun(const mfem::Vector &x, const mfem::Vector &dir, + double jac, const mfem::Vector &u, + const mfem::DenseMatrix &Dw) {return 0.0;} + +private: + /// Reynolds number + double Re; + /// Prandtl number + double Pr; + /// nondimensionalized dynamic viscosity + double mu; + /// work space for flux computations + mfem::Vector work_vec; + /// vector of SA model parameters + mfem::Vector sacs; + + mfem::Vector qfs; +}; + +/// Integrator for SA far-field boundary conditions +/// \tparam dim - number of spatial dimensions (1, 2, or 3) +/// \note This derived class uses the CRTP +template +class SAFarFieldBC : public ViscousBoundaryIntegrator> +{ +public: + /// Constructs an integrator for a SA far-field boundary + /// \param[in] diff_stack - for algorithmic differentiation + /// \param[in] fe_coll - used to determine the face elements + /// \param[in] Re_num - Reynolds number + /// \param[in] Pr_num - Prandtl number + /// \param[in] q_far - state at the far-field + /// \param[in] vis - viscosity (if negative use Sutherland's law) + /// \param[in] a - used to move residual to lhs (1.0) or rhs(-1.0) + SAFarFieldBC(adept::Stack &diff_stack, + const mfem::FiniteElementCollection *fe_coll, + double Re_num, double Pr_num, + const mfem::Vector &q_far, double vis = -1.0, + double a = 1.0) + : ViscousBoundaryIntegrator>( + diff_stack, fe_coll, dim + 3, a), + Re(Re_num), Pr(Pr_num), qfs(q_far), mu(vis), work_vec(dim + 2) {} + + /// converts conservative variables to entropy variables + /// \param[in] q - conservative variables that are to be converted + /// \param[out] w - entropy variables corresponding to `q` + /// \note a wrapper for the relevant function in `euler_fluxes.hpp` + void convertVars(const mfem::Vector &q, mfem::Vector &w) + { + calcEntropyVars(q.GetData(), w.GetData()); + w(dim+2) = q(dim+2); + } + + /// Compute the Jacobian of the mapping `convert` w.r.t. `u` + /// \param[in] q - conservative variables that are to be converted + /// \param[out] dwdu - Jacobian of entropy variables w.r.t. `u` + void convertVarsJacState(const mfem::Vector &q, mfem::DenseMatrix &dwdu) + { + dwdu = 0.0; + convertVarsJac(q, this->stack, dwdu); + dwdu(dim+2,dim+2) = 1.0; + } + + /// Compute flux corresponding to a SA inflow boundary + /// \param[in] x - coordinate location at which flux is evaluated (not used) + /// \param[in] dir - vector normal to the boundary at `x` + /// \param[in] jac - mapping Jacobian (needed by no-slip penalty) + /// \param[in] q - conservative variables at which to evaluate the flux + /// \param[in] Dw - space derivatives of the entropy variables + /// \param[out] flux_vec - value of the flux + void calcFlux(const mfem::Vector &x, const mfem::Vector &dir, double jac, + const mfem::Vector &q, const mfem::DenseMatrix &Dw, + mfem::Vector &flux_vec); + + + /// Compute jacobian of flux corresponding to a SA far-field boundary + /// w.r.t `states` + /// \param[in] x - coordinate location at which flux is evaluated + /// \param[in] dir - vector normal to the boundary at `x` + /// \param[in] jac - mapping Jacobian (needed by no-slip penalty) + /// \param[in] q - conservative variables at which to evaluate the flux + /// \param[in] Dw - space derivatives of the entropy variables (not used yet) + /// \param[out] flux_jac - jacobian of the flux + void calcFluxJacState(const mfem::Vector &x, const mfem::Vector &dir, + double jac, const mfem::Vector &q, + const mfem::DenseMatrix &Dw, + mfem::DenseMatrix &flux_jac); + + /// Compute jacobian of flux corresponding to a SA far-field boundary + /// w.r.t `entrpy-variables' derivatives` + /// \param[in] x - coordinate location at which flux is evaluated (not used) + /// \param[in] dir - vector normal to the boundary at `x` + /// \param[in] jac - mapping Jacobian (needed by no-slip penalty) + /// \param[in] q - conservative variables at which to evaluate the flux + /// \param[in] Dw - space derivatives of the entropy variables + /// \param[out] flux_jac - jacobian of the flux + void calcFluxJacDw(const mfem::Vector &x, const mfem::Vector &dir, + double jac, const mfem::Vector &q, + const mfem::DenseMatrix &Dw, + std::vector &flux_jac); + + /// Compute flux terms that are multiplied by test-function derivative + /// \param[in] x - coordinate location at which flux is evaluated (not used) + /// \param[in] dir - vector normal to the boundary at `x` + /// \param[in] q - conservative variables at which to evaluate the fluxes + /// \param[out] flux_mat[:,di] - to be scaled by derivative `D_[di] v` + void calcFluxDv(const mfem::Vector &x, const mfem::Vector &dir, + const mfem::Vector &q, mfem::DenseMatrix &flux_mat) + { + flux_mat = 0.0; + } + + /// Compute the Jacobian of calcFluxDv w.r.t. state + /// \param[in] x - coordinate location at which fluxes are evaluated + /// \param[in] dir - vector normal to the boundary at `x` + /// \param[in] q - conservative variables at which to evaluate the flux + /// \param[in] flux_jac[di] - Jacobian of calcFluxDv[di] with respect to `q` + void calcFluxDvJacState(const mfem::Vector &x, const mfem::Vector dir, + const mfem::Vector &u, + std::vector &flux_jac) + { + for (int d = 0; d < dim; ++d) + flux_jac[d] = 0.0; + } + + /// Computes boundary node contribution to the surface force + /// \param[in] x - coordinate location at which function is evaluated + /// \param[in] dir - vector normal to the boundary at `x` + /// \param[in] jac - mapping Jacobian determinant (needed by some fluxes) + /// \param[in] u - state at which to evaluate the function + /// \param[in] Dw - `Dw[:,di]` is the derivative of `w` in direction `di` + /// \returns fun - stress at given point + double calcBndryFun(const mfem::Vector &x, const mfem::Vector &dir, + double jac, const mfem::Vector &u, + const mfem::DenseMatrix &Dw) {return 0.0;} + +private: + /// Reynolds number + double Re; + /// Prandtl number + double Pr; + /// nondimensionalized dynamic viscosity + double mu; + /// far-field boundary state + mfem::Vector qfs; + /// work space for flux computations + mfem::Vector work_vec; +}; + + +/// Integrator for SA inflow boundary condition +/// \tparam dim - number of spatial dimensions (1, 2, or 3) +/// \note This derived class uses the CRTP +template +class SAInflowBC : public ViscousBoundaryIntegrator> +{ +public: + /// Constructs an integrator for a SA inflow boundary + /// \param[in] diff_stack - for algorithmic differentiation + /// \param[in] fe_coll - used to determine the face elements + /// \param[in] Re_num - Reynolds number + /// \param[in] Pr_num - Prandtl number + /// \param[in] q_inflow - state at the inflow + /// \param[in] vis - viscosity (if negative use Sutherland's law) + /// \param[in] a - used to move residual to lhs (1.0) or rhs(-1.0) + SAInflowBC(adept::Stack &diff_stack, + const mfem::FiniteElementCollection *fe_coll, + double Re_num, double Pr_num, + const mfem::Vector &q_inflow, double vis = -1.0, + double a = 1.0) + : ViscousBoundaryIntegrator>( + diff_stack, fe_coll, dim + 2, a), + Re(Re_num), Pr(Pr_num), + q_in(q_inflow), mu(vis), work_vec(dim + 2) {} + + /// converts conservative variables to entropy variables + /// \param[in] q - conservative variables that are to be converted + /// \param[out] w - entropy variables corresponding to `q` + /// \note a wrapper for the relevant function in `euler_fluxes.hpp` + void convertVars(const mfem::Vector &q, mfem::Vector &w) + { + calcEntropyVars(q.GetData(), w.GetData()); + } + + /// Compute the Jacobian of the mapping `convert` w.r.t. `u` + /// \param[in] q - conservative variables that are to be converted + /// \param[out] dwdu - Jacobian of entropy variables w.r.t. `u` + void convertVarsJacState(const mfem::Vector &q, mfem::DenseMatrix &dwdu) + { + convertVarsJac(q, this->stack, dwdu); + } + + /// Compute flux corresponding to a SA inflow boundary + /// \param[in] x - coordinate location at which flux is evaluated (not used) + /// \param[in] dir - vector normal to the boundary at `x` + /// \param[in] jac - mapping Jacobian (needed by no-slip penalty) + /// \param[in] q - conservative variables at which to evaluate the flux + /// \param[in] Dw - space derivatives of the entropy variables + /// \param[out] flux_vec - value of the flux + void calcFlux(const mfem::Vector &x, const mfem::Vector &dir, double jac, + const mfem::Vector &q, const mfem::DenseMatrix &Dw, + mfem::Vector &flux_vec); + + /// Compute jacobian of flux corresponding to a SA inflow boundary + /// w.r.t `states` + /// \param[in] x - coordinate location at which flux is evaluated + /// \param[in] dir - vector normal to the boundary at `x` + /// \param[in] jac - mapping Jacobian (needed by no-slip penalty) + /// \param[in] q - conservative variables at which to evaluate the flux + /// \param[in] Dw - space derivatives of the entropy variables (not used yet) + /// \param[out] flux_jac - jacobian of the flux + void calcFluxJacState(const mfem::Vector &x, const mfem::Vector &dir, + double jac, const mfem::Vector &q, + const mfem::DenseMatrix &Dw, + mfem::DenseMatrix &flux_jac); + + /// Compute jacobian of flux corresponding to a SA inflow boundary + /// w.r.t `entrpy-variables' derivatives` + /// \param[in] x - coordinate location at which flux is evaluated + /// \param[in] dir - vector normal to the boundary at `x` + /// \param[in] jac - mapping Jacobian (needed by no-slip penalty) + /// \param[in] q - conservative variables at which to evaluate the flux + /// \param[in] Dw - space derivatives of the entropy variables + /// \param[out] flux_jac - jacobian of the flux + void calcFluxJacDw(const mfem::Vector &x, const mfem::Vector &dir, + double jac, const mfem::Vector &q, + const mfem::DenseMatrix &Dw, + vector &flux_jac); + + /// Compute flux terms that are multiplied by test-function derivative + /// \param[in] x - coordinate location at which flux is evaluated (not used) + /// \param[in] dir - vector normal to the boundary at `x` + /// \param[in] q - conservative variables at which to evaluate the fluxes + /// \param[out] flux_mat[:,di] - to be scaled by derivative `D_[di] v` + void calcFluxDv(const mfem::Vector &x, const mfem::Vector &dir, + const mfem::Vector &q, mfem::DenseMatrix &flux_mat) + { + flux_mat = 0.0; + } + + /// Compute the Jacobian of calcFluxDv w.r.t. state + /// \param[in] x - coordinate location at which fluxes are evaluated + /// \param[in] dir - vector normal to the boundary at `x` + /// \param[in] q - conservative variables at which to evaluate the flux + /// \param[in] flux_jac[di] - Jacobian of calcFluxDv[di] with respect to `q` + void calcFluxDvJacState(const mfem::Vector &x, const mfem::Vector dir, + const mfem::Vector &u, + std::vector &flux_jac) + { + for (int d = 0; d < dim; ++d) + flux_jac[d] = 0.0; + } + + /// Computes boundary node contribution to the surface force + /// \param[in] x - coordinate location at which function is evaluated + /// \param[in] dir - vector normal to the boundary at `x` + /// \param[in] jac - mapping Jacobian determinant (needed by some fluxes) + /// \param[in] u - state at which to evaluate the function + /// \param[in] Dw - `Dw[:,di]` is the derivative of `w` in direction `di` + /// \returns fun - stress at given point + double calcBndryFun(const mfem::Vector &x, const mfem::Vector &dir, + double jac, const mfem::Vector &u, + const mfem::DenseMatrix &Dw) {return 0.0;} + +private: + /// Reynolds number + double Re; + /// Prandtl number + double Pr; + /// nondimensionalized dynamic viscosity + double mu; + /// Inflow boundary state + mfem::Vector q_in; + /// work space for flux computations + mfem::Vector work_vec; +}; + +/// Integrator for SA outflow boundary condition +/// \tparam dim - number of spatial dimensions (1, 2, or 3) +/// \note This derived class uses the CRTP +template +class SAOutflowBC : public ViscousBoundaryIntegrator> +{ +public: + /// Constructs an integrator for a SA outflow boundary + /// \param[in] diff_stack - for algorithmic differentiation + /// \param[in] fe_coll - used to determine the face elements + /// \param[in] Re_num - Reynolds number + /// \param[in] Pr_num - Prandtl number + /// \param[in] q_outflow - state at the outflow + /// \param[in] vis - viscosity (if negative use Sutherland's law) + /// \param[in] a - used to move residual to lhs (1.0) or rhs(-1.0) + SAOutflowBC(adept::Stack &diff_stack, + const mfem::FiniteElementCollection *fe_coll, + double Re_num, double Pr_num, + const mfem::Vector &q_outflow, double vis = -1.0, + double a = 1.0) + : ViscousBoundaryIntegrator>( + diff_stack, fe_coll, dim + 2, a), + Re(Re_num), Pr(Pr_num), + q_out(q_outflow), mu(vis), work_vec(dim + 2) {} + + /// converts conservative variables to entropy variables + /// \param[in] q - conservative variables that are to be converted + /// \param[out] w - entropy variables corresponding to `q` + /// \note a wrapper for the relevant function in `euler_fluxes.hpp` + void convertVars(const mfem::Vector &q, mfem::Vector &w) + { + calcEntropyVars(q.GetData(), w.GetData()); + } + + /// Compute the Jacobian of the mapping `convert` w.r.t. `u` + /// \param[in] q - conservative variables that are to be converted + /// \param[out] dwdu - Jacobian of entropy variables w.r.t. `u` + void convertVarsJacState(const mfem::Vector &q, mfem::DenseMatrix &dwdu) + { + convertVarsJac(q, this->stack, dwdu); + } + + /// Compute flux corresponding to a SA inflow boundary + /// \param[in] x - coordinate location at which flux is evaluated (not used) + /// \param[in] dir - vector normal to the boundary at `x` + /// \param[in] jac - mapping Jacobian (needed by no-slip penalty) + /// \param[in] q - conservative variables at which to evaluate the flux + /// \param[in] Dw - space derivatives of the entropy variables + /// \param[out] flux_vec - value of the flux + void calcFlux(const mfem::Vector &x, const mfem::Vector &dir, double jac, + const mfem::Vector &q, const mfem::DenseMatrix &Dw, + mfem::Vector &flux_vec); + + /// Compute jacobian of flux corresponding to a SA inflow boundary + /// w.r.t `states` + /// \param[in] x - coordinate location at which flux is evaluated + /// \param[in] dir - vector normal to the boundary at `x` + /// \param[in] jac - mapping Jacobian (needed by no-slip penalty) + /// \param[in] q - conservative variables at which to evaluate the flux + /// \param[in] Dw - space derivatives of the entropy variables (not used yet) + /// \param[out] flux_jac - jacobian of the flux + void calcFluxJacState(const mfem::Vector &x, const mfem::Vector &dir, + double jac, const mfem::Vector &q, + const mfem::DenseMatrix &Dw, + mfem::DenseMatrix &flux_jac); + + /// Compute jacobian of flux corresponding to a SA inflow boundary + /// w.r.t `entrpy-variables' derivatives` + /// \param[in] x - coordinate location at which flux is evaluated (not used) + /// \param[in] dir - vector normal to the boundary at `x` + /// \param[in] q - conservative variables at which to evaluate the flux + /// \param[in] Dw - space derivatives of the entropy variables + /// \param[out] flux_jac - jacobian of the flux + void calcFluxJacDw(const mfem::Vector &x, const mfem::Vector &dir, + double jac, const mfem::Vector &q, + const mfem::DenseMatrix &Dw, + std::vector &flux_jac); + + /// Compute flux terms that are multiplied by test-function derivative + /// \param[in] x - coordinate location at which flux is evaluated (not used) + /// \param[in] dir - vector normal to the boundary at `x` + /// \param[in] q - conservative variables at which to evaluate the fluxes + /// \param[out] flux_mat[:,di] - to be scaled by derivative `D_[di] v` + void calcFluxDv(const mfem::Vector &x, const mfem::Vector &dir, + const mfem::Vector &q, mfem::DenseMatrix &flux_mat) + { + flux_mat = 0.0; + } + + /// Compute the Jacobian of calcFluxDv w.r.t. state + /// \param[in] x - coordinate location at which fluxes are evaluated + /// \param[in] dir - vector normal to the boundary at `x` + /// \param[in] q - conservative variables at which to evaluate the flux + /// \param[in] flux_jac[di] - Jacobian of calcFluxDv[di] with respect to `q` + void calcFluxDvJacState(const mfem::Vector &x, const mfem::Vector dir, + const mfem::Vector &u, + std::vector &flux_jac) + { + for (int d = 0; d < dim; ++d) + flux_jac[d] = 0.0; + } + + /// Computes boundary node contribution to the surface force + /// \param[in] x - coordinate location at which function is evaluated + /// \param[in] dir - vector normal to the boundary at `x` + /// \param[in] jac - mapping Jacobian determinant (needed by some fluxes) + /// \param[in] u - state at which to evaluate the function + /// \param[in] Dw - `Dw[:,di]` is the derivative of `w` in direction `di` + /// \returns fun - stress at given point + double calcBndryFun(const mfem::Vector &x, const mfem::Vector &dir, + double jac, const mfem::Vector &u, + const mfem::DenseMatrix &Dw) {return 0.0;} + +private: + /// Reynolds number + double Re; + /// Prandtl number + double Pr; + /// nondimensionalized dynamic viscosity + double mu; + /// Outflow boundary state + mfem::Vector q_out; + /// work space for flux computations + mfem::Vector work_vec; +}; + +/// Integrator for local-projection stabilization with SA variable +/// \tparam dim - number of spatial dimensions (1, 2, or 3) +/// \tparam entvar - if true, the state variables are the entropy variables +/// \note This derived class uses the CRTP +template +class SALPSIntegrator : public LPSIntegrator< + SALPSIntegrator> +{ +public: + /// Construct an LPS integrator with the SA variable + /// \param[in] diff_stack - for algorithmic differentiation + /// \param[in] a - used to move residual to lhs (1.0) or rhs(-1.0) + /// \param[in] coeff - the LPS coefficient + SALPSIntegrator(adept::Stack &diff_stack, double a = 1.0, + double coeff = 1.0) + : LPSIntegrator>( + diff_stack, dim + 3, a, coeff) {} + + /// converts state variables to entropy variables, if necessary + /// \param[in] q - state variables that are to be converted + /// \param[out] w - entropy variables corresponding to `q` + /// \note a wrapper for the relevant function in `euler_fluxes.hpp` + void convertVars(const mfem::Vector &q, mfem::Vector &w); + + /// Compute the Jacobian of the mapping `convert` w.r.t. `u` + /// \param[in] q - conservative variables that are to be converted + /// \param[out] dwdu - Jacobian of entropy variables w.r.t. `u` + void convertVarsJacState(const mfem::Vector &q, mfem::DenseMatrix &dwdu); + + /// Applies the matrix `dQ/dW` to `vec`, and scales by the avg. spectral radius + /// \param[in] adjJ - the adjugate of the mapping Jacobian + /// \param[in] q - the state at which `dQ/dW` and radius are to be evaluated + /// \param[in] vec - the vector being multiplied + /// \param[out] mat_vec - the result of the operation + /// \warning adjJ must be supplied transposed from its `mfem` storage format, + /// so we can use pointer arithmetic to access its rows. + /// \note a wrapper for the relevant function in `euler_fluxes.hpp` + void applyScaling(const mfem::DenseMatrix &adjJ, const mfem::Vector &q, + const mfem::Vector &vec, mfem::Vector &mat_vec); + + /// Computes the Jacobian of the product `A(adjJ,q)*v` w.r.t. `q` + /// \param[in] adjJ - adjugate of the mapping Jacobian + /// \param[in] q - state at which `dQ/dW` and radius are evaluated + /// \param[in] vec - vector that is being multiplied + /// \param[out] mat_vec_jac - Jacobian of product w.r.t. `q` + /// \warning adjJ must be supplied transposed from its `mfem` storage format, + /// so we can use pointer arithmetic to access its rows. + void applyScalingJacState(const mfem::DenseMatrix &adjJ, + const mfem::Vector &q, + const mfem::Vector &vec, + mfem::DenseMatrix &mat_vec_jac); + + /// Computes the Jacobian of the product `A(adjJ,u)*v` w.r.t. `vec` + /// \param[in] adjJ - adjugate of the mapping Jacobian + /// \param[in] q - state at which the symmetric matrix `A` is evaluated + /// \param[out] mat_vec_jac - Jacobian of product w.r.t. `vec` + /// \note `mat_vec_jac` stores derivatives treating `adjJ` is a 1d array. + /// \note The size of `mat_vec_jac` must be set before calling this function + void applyScalingJacV(const mfem::DenseMatrix &adjJ, + const mfem::Vector &q, + mfem::DenseMatrix &mat_vec_jac); + +}; + +/// Volume integrator for Navier-Stokes viscous terms with SA variable +/// \tparam dim - number of spatial dimensions (1, 2, or 3) +/// \note This derived class uses the CRTP +template +class SAViscousIntegrator : public SymmetricViscousIntegrator> +{ +public: + /// Construct an SA viscous integrator + /// \param[in] diff_stack - for algorithmic differentiation + /// \param[in] Re_num - Reynolds number + /// \param[in] Pr_num - Prandtl number + /// \param[in] vis - nondimensional dynamic viscosity (use Sutherland if neg) + /// \param[in] a - used to move residual to lhs (1.0) or rhs(-1.0) + SAViscousIntegrator(adept::Stack &diff_stack, double Re_num, double Pr_num, + mfem::Vector sa_params, double vis = -1.0, double a = 1.0) + : SymmetricViscousIntegrator>( + diff_stack, dim + 3, a), + Re(Re_num), Pr(Pr_num), mu(vis), sacs(sa_params) {} + + /// converts conservative variables to entropy variables + /// \param[in] q - conservative variables that are to be converted + /// \param[out] w - entropy variables corresponding to `q` + /// \note a wrapper for the relevant function in `euler_fluxes.hpp` + void convertVars(const mfem::Vector &q, mfem::Vector &w) + { + calcEntropyVars(q.GetData(), w.GetData()); + w(dim+2) = q(dim+2); + } + + /// Compute the Jacobian of the mapping `convert` w.r.t. `u` + /// \param[in] q - conservative variables that are to be converted + /// \param[out] dwdu - Jacobian of entropy variables w.r.t. `u` + void convertVarsJacState(const mfem::Vector &q, mfem::DenseMatrix &dwdu) + { + dwdu = 0.0; + convertVarsJac(q, this->stack, dwdu); + dwdu(dim+2, dim+2) = 1.0; + } + + /// applies symmetric matrices \f$ C_{d,:}(q) \f$ to input `Dw` + /// \param[in] d - index `d` in \f$ C_{d,:} \f$ matrices + /// \param[in] x - coordinate location at which scaling evaluated (not used) + /// \param[in] q - state at which the symmetric matrices `C` are evaluated + /// \param[in] Dw - `Du[:,d2]` stores derivative of `w` in direction `d2`. + /// \param[out] CDw - product of the multiplication between the `C` and `Dw`. + void applyScaling(int d, const mfem::Vector &x, const mfem::Vector &q, + const mfem::DenseMatrix &Dw, mfem::Vector &CDw); + + /// Computes the Jacobian of the product `C(q)*Dw` w.r.t. `q` + /// \param[in] d - index `d` in \f$ C_{d,:} \f$ matrices + /// \param[in] x - coordinate location at which scaling evaluated (not used) + /// \param[in] q - state at which the symmetric matrix `C` is evaluated + /// \param[in] Dw - vector that is being multiplied + /// \param[out] CDw_jac - Jacobian of product w.r.t. `q` + /// \note This uses the CRTP, so it wraps call to a func. in Derived. + void applyScalingJacState(int d, const mfem::Vector &x, + const mfem::Vector &q, + const mfem::DenseMatrix &Dw, + mfem::DenseMatrix &CDw_jac); + + /// Computes the Jacobian of the product `C(q)*Dw` w.r.t. `Dw` + /// \param[in] d - index `d` in \f$ C_{d,:} \f$ matrices + /// \param[in] x - coordinate location at which scaling evaluated (not used) + /// \param[in] q - state at which the symmetric matrix `C` is evaluated + /// \param[in] Dw - vector that is being multiplied + /// \param[out] CDw_jac - Jacobian of product w.r.t. `Dw` (i.e. `C`) + /// \note This uses the CRTP, so it wraps call to a func. in Derived. + void applyScalingJacDw(int d, const mfem::Vector &x, const mfem::Vector &q, + const mfem::DenseMatrix &Dw, + vector &CDw_jac); + + /// This allows the base class to access the number of dimensions + static const int ndim = dim; + +private: + /// Reynolds number + double Re; + /// Prandtl number + double Pr; + /// nondimensional dynamic viscosity + double mu; + /// vector of SA model parameters + mfem::Vector sacs; + + mfem::Vector qfs; +}; + + +/// Source-term integrator for a 2D SA MMS problem +/// \note For details on the MMS problem, see the file rans_mms.py +class SAMMSIntegrator : public MMSIntegrator +{ +public: + /// Construct an integrator for a 2D RANS SA MMS source + /// \param[in] diff_stack - for algorithmic differentiation + /// \param[in] Re_num - Reynolds number + /// \param[in] Pr_num - Prandtl number + /// \param[in] a - used to move residual to lhs (1.0) or rhs(-1.0) + SAMMSIntegrator(double Re_num, double Pr_num, double a = -1.0) + : MMSIntegrator(5, a), + Re(Re_num), Pr(Pr_num) {} + + /// Computes the MMS source term at a give point + /// \param[in] x - spatial location at which to evaluate the source + /// \param[out] src - source term evaluated at `x` + void calcSource(const mfem::Vector &x, mfem::Vector &src) + { + double mu = 1.0/Re; + calcSAMMS(mu, Pr, x.GetData(), src.GetData()); + } + +private: + /// Reynolds number + double Re; + /// Prandtl number + double Pr; +}; + +#include "rans_integ_def.hpp" + +} // namespace mach + +#endif \ No newline at end of file diff --git a/src/physics/fluidflow/rans_integ_def.hpp b/src/physics/fluidflow/rans_integ_def.hpp new file mode 100644 index 00000000..01cc5b70 --- /dev/null +++ b/src/physics/fluidflow/rans_integ_def.hpp @@ -0,0 +1,1286 @@ +//============================================================================== +// SASourceIntegrator methods +template +void SASourceIntegrator::AssembleElementVector( + const mfem::FiniteElement &fe, mfem::ElementTransformation &Trans, + const mfem::Vector &elfun, mfem::Vector &elvect) +{ + using namespace mfem; + const SBPFiniteElement &sbp = dynamic_cast(fe); + int num_nodes = sbp.GetDof(); +#ifdef MFEM_THREAD_SAFE + Vector ui, uci, xi, uj, grad_nu_i, curl_i; + DenseMatrix grad, curl; +#endif + elvect.SetSize(num_states * num_nodes); + ui.SetSize(num_states); + uci.SetSize(num_states); + xi.SetSize(dim); + grad_nu_i.SetSize(dim); + grad_rho_i.SetSize(dim); + curl_i.SetSize(3); + DenseMatrix uc(elfun.GetData(), num_nodes, num_states); + DenseMatrix u; u = uc; + u.Transpose(); //NOTE: Do not manipulate elfun directly + DenseMatrix Dui(num_states, dim); + Vector Duidi(num_states); + + // convert momentum to velocity + for (int nn = 0; nn < num_nodes; nn++) + { + for (int di = 0; di < dim; di++) + { + u(di+1, nn) *= 1.0/u(0, nn); + } + //u(dim+2, nn) *= 1.0/u(0, nn); + } + + // precompute certain values + elvect = 0.0; + for (int i = 0; i < num_nodes; ++i) + { + // get the Jacobian (Trans.Weight) and cubature weight (node.weight) + const IntegrationPoint &node = fe.GetNodes().IntPoint(i); + Trans.SetIntPoint(&node); + u.GetColumn(i, ui); + //uc.GetRow(i, uci); + Trans.Transform(node, xi); + + // get Dui + Duidi = 0.0; Dui = 0.0; + for (int di = 0; di < dim; di++) + { + sbp.multStrongOperator(di, i, u, Duidi); + Dui.SetCol(di, Duidi); + } + + // compute vorticity at node and take magnitude + calcVorticity(Dui.GetData(), Trans.InverseJacobian().GetData(), curl_i.GetData()); + double S = sqrt((curl_i(0)*curl_i(0) + curl_i(1)*curl_i(1) +curl_i(2)*curl_i(2))); + //S = abs(curl_i(2)); + // compute gradient of turbulent viscosity at node + calcGrad(dim+2, Dui.GetData(), Trans.InverseJacobian().GetData(), grad_nu_i.GetData()); + // compute gradient of density at node + calcGrad(0, Dui.GetData(), Trans.InverseJacobian().GetData(), grad_rho_i.GetData()); + + // get distance function value at node + double d = dist.GetValue(Trans.ElementNo, node); //temp solution, y distance from wall + + // accumulate source terms + double src = calcSAFullSource(ui.GetData(), + mu, d, S, Re, d0, grad_nu_i.GetData(), + grad_rho_i.GetData(), sacs.GetData(), + prod, dest); + + + //byNODES + elvect(i + num_nodes*(num_states-1)) = -alpha * Trans.Weight() * node.weight * src; + //elvect((dim+2) + i*(num_states)) = alpha * Trans.Weight() * node.weight * src; + } // loop over element nodes i +} + +template +void SASourceIntegrator::AssembleElementGrad( + const mfem::FiniteElement &fe, mfem::ElementTransformation &Trans, + const mfem::Vector &elfun, mfem::DenseMatrix &elmat) +{ + using namespace mfem; + const SBPFiniteElement &sbp = dynamic_cast(fe); + int num_nodes = sbp.GetDof(); +#ifdef MFEM_THREAD_SAFE + Vector ui, uci, xi, uj, grad_nu_i, curl_i; +#endif + ui.SetSize(num_states); + uci.SetSize(num_states); + xi.SetSize(dim); + grad_nu_i.SetSize(dim); + grad_rho_i.SetSize(dim); + curl_i.SetSize(3); + DenseMatrix uc(elfun.GetData(), num_nodes, num_states); + DenseMatrix u; u = uc; + u.Transpose(); + DenseMatrix Dui(num_states, dim); + Vector Duidi(num_states); + vector jac_curl(dim); + vector jac_grad_nu(dim); + vector jac_grad_rho(dim); + std::vector sacs_a(sacs.Size()); + adept::set_values(sacs_a.data(), sacs.Size(), sacs.GetData()); + + // each matrix represents d/dx, d/dy, d/dz + for (int d = 0; d < dim; ++d) + { + jac_curl[d].SetSize(3, num_states); + jac_grad_nu[d].SetSize(dim, num_states); + jac_grad_rho[d].SetSize(dim, num_states); + } + + // partial derivative terms + Vector dSdc(3); //partial vorticity mag (S) w.r.t. curl + Vector dSrcdu(num_states); //partial src w.r.t. u + Vector dSrcdS(1); //partial src w.r.t. S + Vector dSrcdgradnu(dim); //partial src w.r.t grad + Vector dSrcdgradrho(dim); //partial src w.r.t grad + Vector work1(num_nodes*num_states); + Vector work2s(num_states); + Vector work2gn(num_states); + Vector work2gr(num_states); + Vector work3(num_nodes); + Vector dnu(num_nodes*num_states); //total derivative + + // convert momentum to velocity + for (int nn = 0; nn < num_nodes; nn++) + { + for (int di = 0; di < dim; di++) + { + u(di+1, nn) *= 1.0/u(0, nn); + } + //u(dim+2, nn) *= 1.0/u(0, nn); + } + + elmat.SetSize(num_states*num_nodes); + elmat = 0.0; + for (int i = 0; i < num_nodes; ++i) + { + dnu = 0.0; + // get the Jacobian (Trans.Weight) and cubature weight (node.weight) + const IntegrationPoint &node = fe.GetNodes().IntPoint(i); + Trans.SetIntPoint(&node); + u.GetColumn(i, ui); + //uc.GetRow(i, uci); + Trans.Transform(node, xi); + + // get Dui, insert each spatial derivative into column matrix + Duidi = 0.0; Dui = 0.0; + for (int di = 0; di < dim; di++) + { + sbp.multStrongOperator(di, i, u, Duidi); + Dui.SetCol(di, Duidi); + } + + int Dui_size = Dui.Height() * Dui.Width(); + // set adept inputs + std::vector ui_a(ui.Size()); + //std::vector Dui_a(Dui_size); + std::vector curl_i_a(curl_i.Size()); + std::vector grad_nu_i_a(grad_nu_i.Size()); + std::vector grad_rho_i_a(grad_rho_i.Size()); + adouble mu_a = mu; + adouble Re_a = Re; + adept::set_values(ui_a.data(), ui.Size(), ui.GetData()); + + // compute vorticity at node + calcVorticity(Dui.GetData(), + Trans.InverseJacobian().GetData(), curl_i.GetData()); + calcVorticityJacDw(this->stack, Dui.GetData(), + Trans.InverseJacobian().GetData(), jac_curl); + + // compute gradient of turbulent viscosity at node + calcGrad(dim+2, Dui.GetData(), + Trans.InverseJacobian().GetData(), grad_nu_i.GetData()); + calcGradJacDw(this->stack, dim+2, Dui.GetData(), + Trans.InverseJacobian().GetData(), jac_grad_nu); + // compute gradient of density at node + calcGrad(0, Dui.GetData(), + Trans.InverseJacobian().GetData(), grad_rho_i.GetData()); //curl.GetRow(i, curl_i); + calcGradJacDw(this->stack, 0, Dui.GetData(), + Trans.InverseJacobian().GetData(), jac_grad_rho); + + // set adept outputs + adept::set_values(curl_i_a.data(), curl_i.Size(), curl_i.GetData()); + adept::set_values(grad_nu_i_a.data(), grad_nu_i.Size(), grad_nu_i.GetData()); + adept::set_values(grad_rho_i_a.data(), grad_rho_i.Size(), grad_rho_i.GetData()); + + // get distance function value at node + adouble d = dist.GetValue(Trans.ElementNo, node); //evaluate the distance gridfunction at the node + adouble d0_a = d0; + // if (d < 1e-12) + // d = d0; + + // vorticity magnitude deriv + this->stack.new_recording(); + adouble S = sqrt((curl_i_a[0]*curl_i_a[0] + curl_i_a[1]*curl_i_a[1] +curl_i_a[2]*curl_i_a[2])); + //S = abs(curl_i_a[2]); + this->stack.independent(curl_i_a.data(), curl_i.Size()); + this->stack.dependent(S); + this->stack.jacobian(dSdc.GetData()); + adouble src; + + dSrcdu = 0.0; + // ui differentiation + this->stack.new_recording(); + src = calcSAFullSource(ui_a.data(), + mu_a, d, S, Re_a, d0_a, grad_nu_i_a.data(), + grad_rho_i_a.data(), sacs_a.data(), + prod, dest); + //src = prod*(uci_a[dim+2]-3)*d; ///!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + this->stack.independent(ui_a.data(), ui.Size()); + this->stack.dependent(src); + this->stack.jacobian(dSrcdu.GetData()); + + dSrcdS = 0.0; + // S differentiation + this->stack.new_recording(); + src = calcSAFullSource(ui_a.data(), + mu_a, d, S, Re_a, d0_a, grad_nu_i_a.data(), + grad_rho_i_a.data(), sacs_a.data(), + prod, dest); + //src = prod*(uci_a[dim+2]-3)*d; ///!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + this->stack.independent(S); + this->stack.dependent(src); + this->stack.jacobian(dSrcdS.GetData()); + + dSrcdgradnu = 0.0; + // grad nu differentiation + this->stack.new_recording(); + src = calcSAFullSource(ui_a.data(), + mu_a, d, S, Re_a, d0_a, grad_nu_i_a.data(), + grad_rho_i_a.data(), sacs_a.data(), + prod, dest); + //src = prod*(uci_a[dim+2]-3)*d; ///!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + this->stack.independent(grad_nu_i_a.data(), grad_nu_i_a.size()); + this->stack.dependent(src); + this->stack.jacobian(dSrcdgradnu.GetData()); + + dSrcdgradrho = 0.0; + // grad rho differentiation + this->stack.new_recording(); + src = calcSAFullSource(ui_a.data(), + mu_a, d, S, Re_a, d0_a, grad_nu_i_a.data(), + grad_rho_i_a.data(), sacs_a.data(), + prod, dest); + //src = prod*(uci_a[dim+2]-3)*d; ///!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + this->stack.independent(grad_rho_i_a.data(), grad_rho_i_a.size()); + this->stack.dependent(src); + this->stack.jacobian(dSrcdgradrho.GetData()); + + // Assemble nu derivative + + // Dui sensitivities + for (int di = 0; di < dim; di++) + { + work1 = 0.0; work2s = 0.0; work2gn = 0.0; work2gr = 0.0; work3 = 0.0; + sbp.getStrongOperator(di, i, work3); + jac_curl[di].MultTranspose(dSdc, work2s); + work2s *= dSrcdS(0); + jac_grad_nu[di].MultTranspose(dSrcdgradnu, work2gn); + jac_grad_rho[di].MultTranspose(dSrcdgradrho, work2gr); + for (int nn = 0; nn < num_nodes; nn++) + { + // assuming ordering BY NODES + for (int ns = 0; ns < num_states; ns++) + { + work1(nn + ns*num_nodes) = work3(nn)*(work2s(ns) + work2gn(ns) + work2gr(ns)); + } + } + dnu += work1; + } + + // Direct ui sensitivities + for(int ns = 0; ns < num_states; ns ++) + { + dnu(i + ns*num_nodes) += dSrcdu(ns); + } + + // node weighting + dnu *= (Trans.Weight() * node.weight); + + // account for momentum conversion + for (int nn = 0; nn < num_nodes; nn++) + { + for (int di = 0; di < dim; di++) + { + dnu(nn) -= dnu(nn + (di+1)*num_nodes)*uc(nn, di+1)/(uc(nn, 0)*uc(nn, 0)); + dnu(nn + (di+1)*num_nodes) *= 1.0/uc(nn, 0); + } + //dnu(nn) -= dnu(nn + (dim+2)*num_nodes)*uc(nn, dim+2)/(uc(nn, 0)*uc(nn, 0)); + //dnu(nn + (dim+2)*num_nodes) *= 1.0/uc(nn, 0); + } + + // Set elmat entry + + //byNODES + //elmat.SetRow((num_states-1)+i*num_states, dnu); + elmat.SetRow(i+num_nodes*(num_states-1), dnu); + } // loop over element nodes i + elmat *= -alpha; +} + +//==================================================================================== +//SA Inviscid Integrator Methods + +template +void SAInviscidIntegrator::calcFlux(int di, const mfem::Vector &qL, + const mfem::Vector &qR, + mfem::Vector &flux) +{ + if (entvar) + { + calcIsmailRoeFluxUsingEntVars(di, qL.GetData(), qR.GetData(), + flux.GetData()); + } + else + { + calcIsmailRoeFlux(di, qL.GetData(), qR.GetData(), + flux.GetData()); + } + //add flux term for SA + //flux(dim+2) = 0.5*(qfs(di+1)/qfs(0) + qfs(di+1)/qfs(0))*0.5*(qL(dim+2) + qR(dim+2)); + flux(dim+2) = 0.5*(qL(di+1)/qL(0) + qR(di+1)/qR(0))*0.5*(qL(dim+2) + qR(dim+2)); +} + +template +void SAInviscidIntegrator::calcFluxJacStates( + int di, const mfem::Vector &qL, const mfem::Vector &qR, + mfem::DenseMatrix &jacL, mfem::DenseMatrix &jacR) +{ + // store the full jacobian in jac + mfem::DenseMatrix jac(dim + 3, 2 * (dim + 3)); + // vector of active input variables + std::vector qL_a(qL.Size()); + std::vector qR_a(qR.Size()); + // initialize adouble inputs + adept::set_values(qL_a.data(), qL.Size(), qL.GetData()); + adept::set_values(qR_a.data(), qR.Size(), qR.GetData()); + // start recording + this->stack.new_recording(); + // create vector of active output variables + std::vector flux_a(qL.Size()); + // run algorithm + if (entvar) + { + mach::calcIsmailRoeFluxUsingEntVars(di, qL_a.data(), + qR_a.data(), + flux_a.data()); + } + else + { + mach::calcIsmailRoeFlux(di, qL_a.data(), + qR_a.data(), flux_a.data()); + } + // identify independent and dependent variables + this->stack.independent(qL_a.data(), qL.Size()); + this->stack.independent(qR_a.data(), qR.Size()); + this->stack.dependent(flux_a.data(), qL.Size()); + // compute and store jacobian in jac + this->stack.jacobian_reverse(jac.GetData()); + // retrieve the jacobian w.r.t left state + jacL.CopyCols(jac, 0, dim + 2); + // retrieve the jacobian w.r.t right state + jacR.CopyCols(jac, dim + 3, 2 * (dim + 3) - 1); + // add SA variable contributions ///NOTE: no real need to do it this way + jacL(dim+2,0) = -0.25*(qL(di+1)/(qL(0)*qL(0)))*(qL(dim+2) + qR(dim+2)); + jacR(dim+2,0) = -0.25*(qR(di+1)/(qR(0)*qR(0)))*(qL(dim+2) + qR(dim+2)); + jacL(dim+2,di+1) = 0.25*(1.0/qL(0))*(qL(dim+2) + qR(dim+2)); + jacR(dim+2,di+1) = 0.25*(1.0/qR(0))*(qL(dim+2) + qR(dim+2)); + jacL(dim+2,dim+2) = 0.25*(qL(di+1)/qL(0) + qR(di+1)/qR(0)); + jacR(dim+2,dim+2) = 0.25*(qL(di+1)/qL(0) + qR(di+1)/qR(0)); + // jacL(dim+2,dim+2) = 0.5*(qfs(di+1)/qfs(0) + qfs(di+1)/qfs(0)); + // jacR(dim+2,dim+2) = 0.5*(qfs(di+1)/qfs(0) + qfs(di+1)/qfs(0)); +} + +//==================================================================================== +//SA Boundary Integrator Methods + +//============================================================================== +//SA No-Slip Adiabatic Wall Integrator methods +template +void SANoSlipAdiabaticWallBC::calcFlux(const mfem::Vector &x, + const mfem::Vector &dir, double jac, + const mfem::Vector &q, + const mfem::DenseMatrix &Dw, + mfem::Vector &flux_vec) +{ + // Step 1: apply the EC slip wall flux + calcSlipWallFlux(x.GetData(), dir.GetData(), q.GetData(), + flux_vec.GetData()); + flux_vec(dim+2) = 0.0; + // Step 2: evaluate the adiabatic flux + double mu_Re = mu; + if (mu < 0.0) + { + mu_Re = calcSutherlandViscosity(q.GetData()); + } + double fv1 = calcSACoefficient(q.GetData(), mu, + sacs.GetData()); + if (q(dim+2)<0) + fv1 = 0.0; + double mu_Re_SA = (mu_Re + q(0)*q(dim+2)*fv1)/Re; + calcAdiabaticWallFlux(dir.GetData(), mu_Re_SA, Pr, q.GetData(), + Dw.GetData(), work_vec.GetData()); + work_vec(dim+2) = 0.0; + flux_vec -= work_vec; // note the minus sign!!! + // evaluate wall normal eddy viscosity flux + double grad[dim]; + for (int di = 0; di < dim; di++) + grad[di] = Dw(dim+2, di); + double SAflux = dot(dir.GetData(), grad); + double fn = calcSANegativeCoefficient(q.GetData(), mu, + sacs.GetData()); + flux_vec(dim+2) -= (mu/q(0) + fn*q(dim+2))*SAflux/(sacs(2)*Re); + // Step 3: evaluate the no-slip penalty + calcNoSlipPenaltyFlux(dir.GetData(), jac, mu_Re_SA, Pr, qfs.GetData(), + q.GetData(), work_vec.GetData()); + work_vec(dim+2) = 0.0; + flux_vec += work_vec; + double dnu = q(dim+2); + double dnuflux = (mu/q(0) + fn*q(dim+2))*dnu/(sacs(2)*Re); + double fac = 1*sqrt(dot(dir, dir))/jac; + flux_vec(dim+2) += dnuflux*fac; +} + +template +void SANoSlipAdiabaticWallBC::calcFluxJacState( + const mfem::Vector &x, const mfem::Vector &dir, double jac, + const mfem::Vector &q, const mfem::DenseMatrix &Dw, + mfem::DenseMatrix &flux_jac) +{ + // create containers for active double objects for each input + int Dw_size = Dw.Height() * Dw.Width(); + std::vector x_a(x.Size()); + std::vector q_a(q.Size()); + std::vector Dw_a(Dw_size); + std::vector work_vec_a(work_vec.Size()); + std::vector dir_a(dir.Size()); + std::vector qfs_a(qfs.Size()); + std::vector sacs_a(sacs.Size()); + // initialize active double containers with data from inputs + adept::set_values(q_a.data(), q.Size(), q.GetData()); + adept::set_values(dir_a.data(), dir.Size(), dir.GetData()); + adept::set_values(Dw_a.data(), Dw_size, Dw.GetData()); + adept::set_values(qfs_a.data(), qfs.Size(), qfs.GetData()); + adept::set_values(x_a.data(), x.Size(), x.GetData()); + adept::set_values(sacs_a.data(), sacs.Size(), sacs.GetData()); + // start new stack recording + this->stack.new_recording(); + // create container for active double flux output + std::vector flux_a(q.Size()); + // Step 1: apply the EC slip wall flux + mach::calcSlipWallFlux(x_a.data(), dir_a.data(), q_a.data(), + flux_a.data()); + flux_a[dim+2] = 0.0; + // Step 2: evaluate the adiabatic flux + adouble mu_Re = mu; + if (mu < 0.0) + { + mu_Re = mach::calcSutherlandViscosity(q_a.data()); + } + //mu_Re /= Re; + adouble fv1 = calcSACoefficient(q_a.data(), mu, + sacs_a.data()); + if (q_a[dim+2]<0) + fv1 = 0.0; + adouble mu_Re_SA = (mu_Re + q_a[0]*q_a[dim+2]*fv1)/Re; + mach::calcAdiabaticWallFlux(dir_a.data(), mu_Re_SA, Pr, q_a.data(), + Dw_a.data(), work_vec_a.data()); + work_vec_a[dim+2] = 0.0; + for (int i = 0; i < flux_a.size(); ++i) + { + flux_a[i] -= work_vec_a[i]; // note the minus sign!!! + } + // evaluate wall normal eddy viscosity flux + adouble grad[dim]; + for (int di = 0; di < dim; di++) + grad[di] = Dw_a[dim+2 + di*(dim+3)]; + adouble SAflux = dot(dir_a.data(), grad); + adouble fn = calcSANegativeCoefficient(q_a.data(), mu, + sacs_a.data()); + flux_a[dim+2] -= (mu/q_a[0] + fn*q_a[dim+2])*SAflux/(sacs_a[2]*Re); + // Step 3: evaluate the no-slip penalty + mach::calcNoSlipPenaltyFlux(dir_a.data(), jac, mu_Re_SA, Pr, qfs_a.data(), + q_a.data(), work_vec_a.data()); + work_vec_a[dim+2] = 0.0; + for (int i = 0; i < flux_a.size(); ++i) + { + flux_a[i] += work_vec_a[i]; + } + adouble dnu = q_a[dim+2]; + adouble dnuflux = (mu/q_a[0] + fn*q_a[dim+2])*dnu/(sacs_a[2]*Re); + adouble fac = 1*sqrt(dot(dir_a.data(), dir_a.data()))/jac; + flux_a[dim+2] += dnuflux*fac; + + this->stack.independent(q_a.data(), q.Size()); + this->stack.dependent(flux_a.data(), q.Size()); + this->stack.jacobian(flux_jac.GetData()); +} + +template +void SANoSlipAdiabaticWallBC::calcFluxJacDw(const mfem::Vector &x, const mfem::Vector &dir, double jac, + const mfem::Vector &q, const mfem::DenseMatrix &Dw, + vector &flux_jac) +{ + // create containers for active double objects for each input + int Dw_size = Dw.Height() * Dw.Width(); + std::vector x_a(x.Size()); + std::vector q_a(q.Size()); + std::vector Dw_a(Dw_size); + std::vector work_vec_a(work_vec.Size()); + std::vector dir_a(dir.Size()); + std::vector qfs_a(qfs.Size()); + std::vector sacs_a(sacs.Size()); + // initialize active double containers with data from inputs + adept::set_values(q_a.data(), q.Size(), q.GetData()); + adept::set_values(dir_a.data(), dir.Size(), dir.GetData()); + adept::set_values(Dw_a.data(), Dw_size, Dw.GetData()); + adept::set_values(qfs_a.data(), qfs.Size(), qfs.GetData()); + adept::set_values(sacs_a.data(), sacs.Size(), sacs.GetData()); + // start new stack recording + this->stack.new_recording(); + // create container for active double flux output + std::vector flux_a(q.Size()); + // Step 1: apply the EC slip wall flux + mach::calcSlipWallFlux(x_a.data(), dir_a.data(), q_a.data(), + flux_a.data()); + flux_a[dim+2] = 0.0; + // Step 2: evaluate the adiabatic flux + adouble mu_Re = mu; + if (mu < 0.0) + { + mu_Re = mach::calcSutherlandViscosity(q_a.data()); + } + //mu_Re /= Re; + adouble fv1 = calcSACoefficient(q_a.data(), mu, + sacs_a.data()); + if (q_a[dim+2]<0) + fv1 = 0.0; + adouble mu_Re_SA = (mu_Re + q_a[0]*q_a[dim+2]*fv1)/Re; + mach::calcAdiabaticWallFlux(dir_a.data(), mu_Re_SA, Pr, q_a.data(), + Dw_a.data(), work_vec_a.data()); + work_vec_a[dim+2] = 0.0; + for (int i = 0; i < flux_a.size(); ++i) + { + flux_a[i] -= work_vec_a[i]; // note the minus sign!!! + } + // evaluate wall normal eddy viscosity flux + adouble grad[dim]; + for (int di = 0; di < dim; di++) + grad[di] = Dw_a[dim+2 + di*(dim+3)]; + adouble SAflux = dot(dir_a.data(), grad); + adouble fn = calcSANegativeCoefficient(q_a.data(), mu, + sacs_a.data()); + flux_a[dim+2] -= (mu/q_a[0] + fn*q_a[dim+2])*SAflux/(sacs_a[2]*Re); + // Step 3: evaluate the no-slip penalty + mach::calcNoSlipPenaltyFlux(dir_a.data(), jac, mu_Re_SA, Pr, qfs_a.data(), + q_a.data(), work_vec_a.data()); + work_vec_a[dim+2] = 0.0; + for (int i = 0; i < flux_a.size(); ++i) + { + flux_a[i] += work_vec_a[i]; + } + adouble dnu = q_a[dim+2]; + adouble dnuflux = (mu/q_a[0] + fn*q_a[dim+2])*dnu/(sacs_a[2]*Re); + adouble fac = 1*sqrt(dot(dir_a.data(), dir_a.data()))/jac; + flux_a[dim+2] += dnuflux*fac; + + this->stack.independent(Dw_a.data(),Dw_size); + this->stack.dependent(flux_a.data(), q.Size()); + // compute and store jacobian in CDw_jac + mfem::Vector work(dim*this->num_states*this->num_states); + this->stack.jacobian(work.GetData()); + for (int i = 0; i < dim; ++i) + { + flux_jac[i] = (work.GetData() + i*this->num_states*this->num_states); + } +} + +template +void SANoSlipAdiabaticWallBC::calcFluxDv(const mfem::Vector &x, + const mfem::Vector &dir, + const mfem::Vector &q, + mfem::DenseMatrix &flux_mat) +{ + double mu_Re = mu; + if (mu < 0.0) + mu_Re = calcSutherlandViscosity(q.GetData()); + //mu_Re /= Re; + double fv1 = calcSACoefficient(q.GetData(), mu, + sacs.GetData()); + if (q(dim+2)<0) + fv1 = 0.0; + double mu_Re_SA = (mu_Re + q(0)*q(dim+2)*fv1)/Re; + double fn = calcSANegativeCoefficient(q.GetData(), mu, + sacs.GetData()); + double mu_2 = (mu_Re/q(0) + q(dim+2)*fn)/Re; + calcNoSlipDualFluxSA(dir.GetData(), mu_Re_SA, mu_2, Pr, q.GetData(), + flux_mat.GetData()); +} + +template +void SANoSlipAdiabaticWallBC::calcFluxDvJacState( + const mfem::Vector &x, const mfem::Vector dir, const mfem::Vector &q, + std::vector &flux_jac) +{ + // create containers for active double objects for each input + int flux_size = dim*(dim+3); + std::vector q_a(q.Size()); + std::vector dir_a(dir.Size()); + std::vector sacs_a(sacs.Size()); + // initialize active double containers with data from inputs + adept::set_values(q_a.data(), q.Size(), q.GetData()); + adept::set_values(dir_a.data(), dir.Size(), dir.GetData()); + adept::set_values(sacs_a.data(), sacs.Size(), sacs.GetData()); + // start new stack recording + this->stack.new_recording(); + // create container for active double flux output + std::vector fluxes_a(flux_size); + // evaluate the fluxes + adouble mu_Re = mu; + if (mu < 0.0) + mu_Re = calcSutherlandViscosity(q_a.data()); + //mu_Re /= Re; + adouble fv1 = calcSACoefficient(q_a.data(), mu, + sacs_a.data()); + if (q_a[dim+2]<0) + fv1 = 0.0; + adouble mu_Re_SA = (mu_Re + q_a[0]*q_a[dim+2]*fv1)/Re; + adouble fn = calcSANegativeCoefficient(q_a.data(), mu, + sacs_a.data()); + adouble mu_2 = (mu_Re/q_a[0] + q_a[dim+2]*fn)/Re; + calcNoSlipDualFluxSA(dir_a.data(), mu_Re_SA, mu_2, Pr, q_a.data(), + fluxes_a.data()); + + this->stack.independent(q_a.data(), q.Size()); + this->stack.dependent(fluxes_a.data(), flux_size); + // compute and store jacobian in flux_jac + mfem::Vector work(flux_size*(dim+3)); + this->stack.jacobian(work.GetData()); + for (int s = 0; s < dim+3; ++s) + { + for (int i = 0; i < dim; ++i) + { + flux_jac[i].SetCol(s, work.GetData() + (s*dim + i)*(dim+3)); + } + } +} + + +//============================================================================== +//SA Viscous Slip-Wall Integrator methods +template +void SAViscousSlipWallBC::calcFlux(const mfem::Vector &x, + const mfem::Vector &dir, double jac, + const mfem::Vector &q, + const mfem::DenseMatrix &Dw, + mfem::Vector &flux_vec) +{ + // Part 1: apply the inviscid slip wall BCs + calcSlipWallFlux(x.GetData(), dir.GetData(), q.GetData(), + flux_vec.GetData()); + flux_vec(dim+2) = 0.0; + int Dw_size = Dw.Height() * Dw.Width(); + mfem::Vector Dw_work(Dw_size); + setZeroNormalDeriv(dir.GetData(), Dw.GetData(), + Dw_work.GetData()); + // SA treatment + double nrm[dim]; + double Dw_nrm = 0.0; + double fac = 1.0 / sqrt(dot(dir, dir)); + for (int i = 0; i < dim; ++i) + { + nrm[i] = dir(i) * fac; + Dw_nrm += Dw(dim+2, i)*nrm[i]; + } + for (int i = 0; i < dim; ++i) + { + Dw_work(dim+2 + i*(dim+3)) = Dw(dim+2, i) - nrm[i]*Dw_nrm; + } + + // Part 2: viscous BCs + double mu_Re = mu; + if (mu < 0.0) + { + mu_Re = calcSutherlandViscosity(q.GetData()); + } + //mu_Re /= Re; + double fv1 = calcSACoefficient(q.GetData(), mu, + sacs.GetData()); + if (q(dim+2)<0) + fv1 = 0.0; + double mu_Re_SA = (mu_Re + q(0)*q(dim+2)*fv1)/Re; + for (int d = 0; d < dim; ++d) + { + work_vec = 0.0; + applyViscousScalingSA(d, mu_Re_SA, Pr, q.GetData(), + Dw_work.GetData(), work_vec.GetData()); + work_vec *= dir(d); + flux_vec -= work_vec; + double fn = calcSANegativeCoefficient(q.GetData(), mu, + sacs.GetData()); + flux_vec(dim+2) -= (mu/q(0) + fn*q(dim+2))*Dw_work(dim+2 + d*(dim+3))/(sacs(2)*Re); + } +} + +template +void SAViscousSlipWallBC::calcFluxJacState( + const mfem::Vector &x, const mfem::Vector &dir, double jac, + const mfem::Vector &q, const mfem::DenseMatrix &Dw, + mfem::DenseMatrix &flux_jac) +{ + // create containers for active double objects for each input + int Dw_size = Dw.Height() * Dw.Width(); + std::vector x_a(x.Size()); + std::vector q_a(q.Size()); + std::vector Dw_a(Dw_size); + std::vector work_vec_a(work_vec.Size()); + std::vector dir_a(dir.Size()); + std::vector sacs_a(sacs.Size()); + // initialize active double containers with data from inputs + adept::set_values(q_a.data(), q.Size(), q.GetData()); + adept::set_values(dir_a.data(), dir.Size(), dir.GetData()); + adept::set_values(Dw_a.data(), Dw_size, Dw.GetData()); + adept::set_values(x_a.data(), x.Size(), x.GetData()); + adept::set_values(sacs_a.data(), sacs.Size(), sacs.GetData()); + // start new stack recording + this->stack.new_recording(); + // create container for active double flux output + std::vector flux_a(q.Size()); + // Part 1: apply the inviscid slip wall BCs + calcSlipWallFlux(x_a.data(), dir_a.data(), q_a.data(), + flux_a.data()); + flux_a[dim+2] = 0.0; + std::vector Dw_work_a(Dw_size); + setZeroNormalDeriv(dir_a.data(), Dw_a.data(), + Dw_work_a.data()); + // SA treatment + adouble nrm[dim]; + adouble Dw_nrm = 0.0; + adouble fac = 1.0 / sqrt(dot(dir_a.data(), dir_a.data())); + for (int i = 0; i < dim; ++i) + { + nrm[i] = dir_a[i] * fac; + Dw_nrm += Dw_a[dim+2 + i*(dim+3)]*nrm[i]; + } + for (int i = 0; i < dim; ++i) + { + Dw_work_a[dim+2 + i*(dim+3)] = Dw_a[dim+2 + i*(dim+3)] - nrm[i]*Dw_nrm; + } + + // Part 2: viscous BCs + adouble mu_Re = mu; + if (mu < 0.0) + { + mu_Re = calcSutherlandViscosity(q_a.data()); + } + //mu_Re /= Re; + adouble fv1 = calcSACoefficient(q_a.data(), mu, + sacs_a.data()); + if (q_a[dim+2]<0) + fv1 = 0.0; + adouble mu_Re_SA = (mu_Re + q_a[0]*q_a[dim+2]*fv1)/Re; + for (int d = 0; d < dim; ++d) + { + work_vec_a[dim+2] = 0.0; + applyViscousScalingSA(d, mu_Re_SA, Pr, q_a.data(), + Dw_work_a.data(), work_vec_a.data()); + for (int k = 0; k < dim+2; k++) + { + work_vec_a[k] *= dir_a[d]; + flux_a[k] -= work_vec_a[k]; + } + adouble fn = calcSANegativeCoefficient(q_a.data(), mu, + sacs_a.data()); + flux_a[dim+2] -= (mu/q_a[0] + fn*q_a[dim+2])*Dw_work_a[dim+2 + d*(dim+3)]/(sacs_a[2]*Re); + } + this->stack.independent(q_a.data(), q.Size()); + this->stack.dependent(flux_a.data(), q.Size()); + this->stack.jacobian(flux_jac.GetData()); +} + +template +void SAViscousSlipWallBC::calcFluxJacDw(const mfem::Vector &x, const mfem::Vector &dir, double jac, + const mfem::Vector &q, const mfem::DenseMatrix &Dw, + vector &flux_jac) +{ + // create containers for active double objects for each input + int Dw_size = Dw.Height() * Dw.Width(); + std::vector x_a(x.Size()); + std::vector q_a(q.Size()); + std::vector Dw_a(Dw_size); + std::vector work_vec_a(work_vec.Size()); + std::vector dir_a(dir.Size()); + std::vector sacs_a(sacs.Size()); + // initialize active double containers with data from inputs + adept::set_values(q_a.data(), q.Size(), q.GetData()); + adept::set_values(dir_a.data(), dir.Size(), dir.GetData()); + adept::set_values(Dw_a.data(), Dw_size, Dw.GetData()); + adept::set_values(sacs_a.data(), sacs.Size(), sacs.GetData()); + // start new stack recording + this->stack.new_recording(); + // create container for active double flux output + std::vector flux_a(q.Size()); + // Part 1: apply the inviscid slip wall BCs + calcSlipWallFlux(x_a.data(), dir_a.data(), q_a.data(), + flux_a.data()); + flux_a[dim+2] = 0.0; + std::vector Dw_work_a(Dw_size); + setZeroNormalDeriv(dir_a.data(), Dw_a.data(), + Dw_work_a.data()); + // SA treatment + adouble nrm[dim]; + adouble Dw_nrm = 0.0; + adouble fac = 1.0 / sqrt(dot(dir_a.data(), dir_a.data())); + for (int i = 0; i < dim; ++i) + { + nrm[i] = dir_a[i] * fac; + Dw_nrm += Dw_a[dim+2 + i*(dim+3)]*nrm[i]; + } + for (int i = 0; i < dim; ++i) + { + Dw_work_a[dim+2 + i*(dim+3)] = Dw_a[dim+2 + i*(dim+3)] - nrm[i]*Dw_nrm; + } + + // Part 2: viscous BCs + adouble mu_Re = mu; + if (mu < 0.0) + { + mu_Re = calcSutherlandViscosity(q_a.data()); + } + //mu_Re /= Re; + adouble fv1 = calcSACoefficient(q_a.data(), mu, + sacs_a.data()); + if (q_a[dim+2]<0) + fv1 = 0.0; + adouble mu_Re_SA = (mu_Re + q_a[0]*q_a[dim+2]*fv1)/Re; + for (int d = 0; d < dim; ++d) + { + work_vec_a[dim+2] = 0.0; + applyViscousScalingSA(d, mu_Re_SA, Pr, q_a.data(), + Dw_work_a.data(), work_vec_a.data()); + for (int k = 0; k < dim+2; k++) + { + work_vec_a[k] *= dir_a[d]; + flux_a[k] -= work_vec_a[k]; + } + adouble fn = calcSANegativeCoefficient(q_a.data(), mu, + sacs_a.data()); + flux_a[dim+2] -= (mu/q_a[0] + fn*q_a[dim+2])*Dw_work_a[dim+2 + d*(dim+3)]/(sacs_a[2]*Re); + } + + this->stack.independent(Dw_a.data(),Dw_size); + this->stack.dependent(flux_a.data(), q.Size()); + // compute and store jacobian in CDw_jac + mfem::Vector work(dim*this->num_states*this->num_states); + this->stack.jacobian(work.GetData()); + for (int i = 0; i < dim; ++i) + { + flux_jac[i] = (work.GetData() + i*this->num_states*this->num_states); + } +} + +//==================================================================================== +//SA Far Field Integrator Methods +template +void SAFarFieldBC::calcFlux(const mfem::Vector &x, const mfem::Vector &dir, double jac, + const mfem::Vector &q, const mfem::DenseMatrix &Dw, + mfem::Vector &flux_vec) +{ + calcBoundaryFlux(dir.GetData(), qfs.GetData(), q.GetData(), + work_vec.GetData(), flux_vec.GetData()); + +#if 0 + // Part 2: evaluate the adiabatic flux + double mu_Re = mu; + if (mu < 0.0) + { + mu_Re = calcSutherlandViscosity(q.GetData()); + } + //mu_Re /= Re; + adouble mu_Re_SA = (mu_Re + q_a[dim+2]*fv1)/Re; + for (int d = 0; d < dim; ++d) + { + work_vec = 0.0; + applyViscousScalingSA(d, mu_Re_SA, Pr, q.GetData(), Dw.GetData(), + work_vec.GetData()); + work_vec *= dir[d]; + flux_vec -= work_vec; + } +#endif + + //flux_vec(dim+2) = 0; + // handle SA variable + double Unrm = dot(dir.GetData(), qfs.GetData()+1); + double fac = 1.0;//sqrt(dot(dir, dir)); + if (Unrm > 0.0) + { + + flux_vec(dim+2) = Unrm*q(dim+2)/(q(0)*fac); + } + else + { + //double dq = q(dim+2) - qfs(dim+2); + flux_vec(dim+2) = Unrm*qfs(dim+2)/(qfs(0)*fac);//Unrm*dq/fac; + } +} + + +template +void SAFarFieldBC::calcFluxJacState( + const mfem::Vector &x, const mfem::Vector &dir, double jac, + const mfem::Vector &q, const mfem::DenseMatrix &Dw, + mfem::DenseMatrix &flux_jac) +{ + flux_jac = 0.0; +// mach::calcFluxJacState(x, dir, jac, q, Dw, qfs, work_vec, this->stack, flux_jac); + + int Dw_size = Dw.Height() * Dw.Width(); + // create containers for active double objects for each input + std::vector q_a(q.Size()); + std::vector dir_a(dir.Size()); + std::vector x_a(x.Size()); + std::vector Dw_a(Dw_size); + std::vector qfs_a(q.Size()); + std::vector work_vec_a(work_vec.Size()); + // initialize active double containers with data from inputs + adept::set_values(q_a.data(), q.Size(), q.GetData()); + adept::set_values(dir_a.data(), dir.Size(), dir.GetData()); + adept::set_values(x_a.data(), x.Size(), x.GetData()); + adept::set_values(Dw_a.data(), Dw_size, Dw.GetData()); + adept::set_values(qfs_a.data(), qfs.Size(), qfs.GetData()); + // start new stack recording + this->stack.new_recording(); + // create container for active double flux output + std::vector flux_a(q.Size()); + // Part 1: apply the inviscid inflow boundary condition + mach::calcBoundaryFlux(dir_a.data(), qfs_a.data(), q_a.data(), + work_vec_a.data(), flux_a.data()); +#if 0 + adouble mu_Re = mu; + if (mu < 0.0) + { + mu_Re = mach::calcSutherlandViscosity(q_a.data()); + } + mu_Re /= Re; + adouble mu_Re_SA = mu_Re + q_a[dim+2]*fv1/Re; + for (int d = 0; d < dim; ++d) + { + for (int i = 0; i < work_vec_a.size(); ++i) + { + work_vec_a[i] = 0.0; + } + applyViscousScalingSA(d, mu_Re_SA, Pr, q_a.data(), Dw_a.data(), + work_vec_a.data()); + for (int i = 0; i < flux_a.size(); ++i) + { + work_vec_a[i] *= dir_a[d]; + flux_a[i] -= work_vec_a[i]; // note the minus sign!!! + } + } +#endif + + adouble Unrm = dot(dir_a.data(), qfs_a.data()+1); + adouble fac = 1.0;//sqrt(dot(dir_a.data(), dir_a.data())); + if (Unrm > 0.0) + { + flux_a[dim+2] = Unrm*q_a[dim+2]/(q_a[0]*fac); + } + else + { + //adouble dq = q_a[dim+2] - qfs_a[dim+2]; + flux_a[dim+2] = Unrm*qfs_a[dim+2]/(qfs_a[0]*fac); //Unrm*dq/fac; + } + + this->stack.independent(q_a.data(), q.Size()); + this->stack.dependent(flux_a.data(), q.Size()); + this->stack.jacobian(flux_jac.GetData()); + + // handle SA variable + // double Unrm = dot(dir.GetData(), qfs.GetData()+1); + // double fac = sqrt(dot(dir, dir)); + // if (Unrm > 0.0) + // { + // flux_jac(dim+2, 0) = -Unrm*q(dim+2)/(q(0)*q(0)*fac); + // for(int di = 0; di < dim; di++) + // { + // flux_jac(dim+2, di+1) = dir(di)*q(dim+2)/(q(0)*fac); + // } + // flux_jac(dim+2, dim+2) = Unrm/(q(0)*fac); + // } + // else + // { + // // for(int di = 0; di < dim; di++) + // // { + // // flux_jac(dim+2, di+1) = dir(di)*qfs(dim+2)/qfs(0); + // // } + // } +} + +template +void SAFarFieldBC::calcFluxJacDw( + const mfem::Vector &x, const mfem::Vector &dir, double jac, + const mfem::Vector &q, const mfem::DenseMatrix &Dw, + vector &flux_jac) +{ + // Presently, this BC has no dependence on the derivative + for (int i = 0; i < dim; ++i) + { + flux_jac[i] = 0.0; + } +} + +//==================================================================================== +//SA LPS Integrator Methods + +template +void SALPSIntegrator::convertVars( + const mfem::Vector &q, mfem::Vector &w) +{ + // This conditional should have no overhead, if the compiler is good + if (entvar) + { + w = q; + } + else + { + calcEntropyVars(q.GetData(), w.GetData()); + // SA variable copy + w(dim+2) = q(dim+2); + } +} + +template +void SALPSIntegrator::convertVarsJacState( + const mfem::Vector &q, mfem::DenseMatrix &dwdu) +{ + if (entvar) + { + dwdu = 0.0; + for (int i = 0; i < dim+3; ++i) + { + dwdu(i,i) = 1.0; + } + } + else + { + dwdu = 0.0; + // vector of active input variables + std::vector q_a(q.Size()); + // initialize adouble inputs + adept::set_values(q_a.data(), q.Size(), q.GetData()); + // start recording + this->stack.new_recording(); + // create vector of active output variables + std::vector w_a(q.Size()); + // run algorithm + calcEntropyVars(q_a.data(), w_a.data()); + // identify independent and dependent variables + this->stack.independent(q_a.data(), q.Size()); + this->stack.dependent(w_a.data(), q.Size()); + // compute and store jacobian in dwdu + this->stack.jacobian(dwdu.GetData()); + // SA Variable Derivative + dwdu(dim+2, dim+2) = 1.0; + } +} + +template +void SALPSIntegrator::applyScaling( + const mfem::DenseMatrix &adjJ, const mfem::Vector &q, + const mfem::Vector &vec, mfem::Vector &mat_vec) +{ + if (entvar) + { + applyLPSScalingUsingEntVars(adjJ.GetData(), q.GetData(), + vec.GetData(), mat_vec.GetData()); + throw MachException("Entropy variables not yet supported"); + } + else + { + applyLPSScaling(adjJ.GetData(), q.GetData(), vec.GetData(), + mat_vec.GetData()); + // SA Variable + double U = sqrt(dot(q.GetData()+1, q.GetData()+1))/q(0); + mat_vec(dim+2) = U*vec(dim+2); + } +} + +template +void SALPSIntegrator::applyScalingJacState( + const mfem::DenseMatrix &adjJ, const mfem::Vector &q, + const mfem::Vector &vec, mfem::DenseMatrix &mat_vec_jac) +{ + // declare vectors of active input variables + int adjJ_a_size = adjJ.Height() * adjJ.Width(); + std::vector adjJ_a(adjJ_a_size); + std::vector q_a(q.Size()); + std::vector vec_a(vec.Size()); + // copy data from mfem::Vector + adept::set_values(adjJ_a.data(), adjJ_a_size, adjJ.GetData()); + adept::set_values(q_a.data(), q.Size(), q.GetData()); + adept::set_values(vec_a.data(), vec.Size(), vec.GetData()); + // start recording + this->stack.new_recording(); + // the dependent variable must be declared after the recording + std::vector mat_vec_a(q.Size()); + if (entvar) + { + applyLPSScalingUsingEntVars(adjJ_a.data(), q_a.data(), + vec_a.data(), mat_vec_a.data()); + } + else + { + applyLPSScaling(adjJ_a.data(), q_a.data(), + vec_a.data(), mat_vec_a.data()); + adouble U = sqrt(dot(q_a.data()+1, q_a.data()+1))/q_a[0]; + mat_vec_a[dim+2] = U*vec_a[dim+2]; + } + // set the independent and dependent variable + this->stack.independent(q_a.data(), q.Size()); + this->stack.dependent(mat_vec_a.data(), q.Size()); + // Calculate the jabobian + this->stack.jacobian(mat_vec_jac.GetData()); +} + +template +void SALPSIntegrator::applyScalingJacV( + const mfem::DenseMatrix &adjJ, const mfem::Vector &q, + mfem::DenseMatrix &mat_vec_jac) +{ + // declare vectors of active input variables + int adjJ_a_size = adjJ.Height() * adjJ.Width(); + std::vector adjJ_a(adjJ_a_size); + std::vector q_a(q.Size()); + std::vector vec_a(q.Size()); + // copy data from mfem::Vector + adept::set_values(adjJ_a.data(), adjJ_a_size, adjJ.GetData()); + adept::set_values(q_a.data(), q.Size(), q.GetData()); + // dependence on vec is linear, so any value is ok; use q + adept::set_values(vec_a.data(), q.Size(), q.GetData()); + // start recording + this->stack.new_recording(); + // the dependent variable must be declared after the recording + std::vector mat_vec_a(q.Size()); + if (entvar) + { + // applyLPSScalingUsingEntVars(adjJ_a.data(), q_a.data(), + // vec_a.data(), mat_vec_a.data()); + throw MachException("Entropy variables not yet supported"); + } + else + { + applyLPSScaling(adjJ_a.data(), q_a.data(), + vec_a.data(), mat_vec_a.data()); + adouble U = sqrt(dot(q_a.data()+1, q_a.data()+1))/q_a[0]; + mat_vec_a[dim+2] = U*vec_a[dim+2]; + } + // set the independent and dependent variable + this->stack.independent(vec_a.data(), q.Size()); + this->stack.dependent(mat_vec_a.data(), q.Size()); + // Calculate the jabobian + this->stack.jacobian(mat_vec_jac.GetData()); +} + +//============================================================================== +// SAViscousIntegrator methods +template +void SAViscousIntegrator::applyScaling(int d, const mfem::Vector &x, + const mfem::Vector &q, const mfem::DenseMatrix &Dw, + mfem::Vector &CDw) +{ + double mu_Re = mu; + if (mu < 0.0) + { + mu_Re = calcSutherlandViscosity(q.GetData()); + } + //mu_Re /= Re; + double fv1 = calcSACoefficient(q.GetData(), mu, + sacs.GetData()); + if (q(dim+2)<0) + fv1 = 0.0; + double mu_Re_SA = (mu_Re + q(0)*q(dim+2)*fv1)/Re; + applyViscousScalingSA(d, mu_Re_SA, Pr, q.GetData(), Dw.GetData(), + CDw.GetData()); + double fn = calcSANegativeCoefficient(q.GetData(), mu, + sacs.GetData()); + CDw(dim+2) = (mu/q(0) + fn*q(dim+2))*Dw(dim+2, d)/(sacs(2)*Re); + // if (mu + fn*q(dim+2) < 0) + // cout << "Negative diffusion!: " < +void SAViscousIntegrator::applyScalingJacState(int d, const mfem::Vector &x, + const mfem::Vector &q, const mfem::DenseMatrix &Dw, + mfem::DenseMatrix &CDw_jac) +{ + // vector of active input variables + int Dw_size = Dw.Height() * Dw.Width(); + std::vector q_a(q.Size()); + std::vector Dw_a(Dw_size); + std::vector sacs_a(sacs.Size()); + // initialize adouble inputs + adept::set_values(q_a.data(), q.Size(), q.GetData()); + adept::set_values(Dw_a.data(), Dw_size, Dw.GetData()); + adept::set_values(sacs_a.data(), sacs.Size(), sacs.GetData()); + // start recording + this->stack.new_recording(); + // create vector of active output variables + std::vector CDw_a(q.Size()); + // run algorithm + adouble mu_Re = mu; + if (mu < 0.0) + { + mu_Re = calcSutherlandViscosity(q_a.data()); + } + //mu_Re /= Re; + adouble fv1 = calcSACoefficient(q_a.data(), mu, + sacs_a.data()); + if (q_a[dim+2]<0) + fv1 = 0.0; + adouble mu_Re_SA = (mu_Re + q_a[0]*q_a[dim+2]*fv1)/Re; + applyViscousScalingSA(d, mu_Re_SA, Pr, q_a.data(), Dw_a.data(), + CDw_a.data()); + adouble fn = calcSANegativeCoefficient(q_a.data(), mu, + sacs_a.data()); + CDw_a[dim+2] = (mu/q_a[0] + fn*q_a[dim+2])*Dw_a[dim+2 + d*(dim+3)]/(sacs_a[2]*Re); + // identify independent and dependent variables + this->stack.independent(q_a.data(), q.Size()); + this->stack.dependent(CDw_a.data(), q.Size()); + // compute and store jacobian in CDw_jac + this->stack.jacobian(CDw_jac.GetData()); +} + +template +void SAViscousIntegrator::applyScalingJacDw( + int d, const mfem::Vector &x, const mfem::Vector &q, + const mfem::DenseMatrix &Dw, vector &CDw_jac) +{ + // vector of active input variables + int Dw_size = Dw.Height() * Dw.Width(); + std::vector q_a(q.Size()); + std::vector Dw_a(Dw_size); + std::vector sacs_a(sacs.Size()); + // initialize adouble inputs + adept::set_values(q_a.data(), q.Size(), q.GetData()); + adept::set_values(Dw_a.data(), Dw_size, Dw.GetData()); + adept::set_values(sacs_a.data(), sacs.Size(), sacs.GetData()); + // start recording + this->stack.new_recording(); + // create vector of active output variables + std::vector CDw_a(q.Size()); + // run algorithm + adouble mu_Re = mu; + if (mu < 0.0) + { + mu_Re = calcSutherlandViscosity(q_a.data()); + } + //mu_Re /= Re; + adouble fv1 = calcSACoefficient(q_a.data(), mu, + sacs_a.data()); + if (q_a[dim+2]<0) + fv1 = 0.0; + adouble mu_Re_SA = (mu_Re + q_a[0]*q_a[dim+2]*fv1)/Re; + applyViscousScalingSA(d, mu_Re_SA, Pr, q_a.data(), Dw_a.data(), + CDw_a.data()); + adouble fn = calcSANegativeCoefficient(q_a.data(), mu, + sacs_a.data()); + CDw_a[dim+2] = (mu/q_a[0] + fn*q_a[dim+2])*Dw_a[dim+2 + d*(dim+3)]/(sacs_a[2]*Re); + // identify independent and dependent variables + this->stack.independent(Dw_a.data(), Dw_size); + this->stack.dependent(CDw_a.data(), q.Size()); + // compute and store jacobian in CDw_jac + mfem::Vector work(dim*this->num_states*this->num_states); + this->stack.jacobian(work.GetData()); + for (int i = 0; i < dim; ++i) + { + CDw_jac[i] = (work.GetData() + i*this->num_states*this->num_states); + } +} \ No newline at end of file diff --git a/src/physics/fluidflow/rans_sens_integ.hpp b/src/physics/fluidflow/rans_sens_integ.hpp new file mode 100644 index 00000000..6a0b49e3 --- /dev/null +++ b/src/physics/fluidflow/rans_sens_integ.hpp @@ -0,0 +1,81 @@ +#ifndef MACH_RANS_SENS_INTEG +#define MACH_RANS_SENS_INTEG + +#include "adept.h" +#include "mfem.hpp" + +#include "mesh_sens_integ.hpp" +#include "rans_fluxes.hpp" + +using adept::adouble; + +namespace mach +{ + +/// Integrator for the mesh sensitivity of the SA inviscid integrators +/// \tparam dim - number of spatial dimensions (1, 2, or 3) +/// \tparam entvar - if true, the state variables are the entropy variables +/// \note This derived class uses the CRTP +template +class SAInviscidMeshSensIntegrator : public DyadicMeshSensIntegrator< + SAInviscidMeshSensIntegrator> +{ +public: + /// Construct an integrator for the Ismail-Roe flux over domains + /// \param[in] a - factor, usually used to move terms to rhs + SAInviscidMeshSensIntegrator(const mfem::GridFunction &state_vec, + const mfem::GridFunction &adjoint_vec, + int num_state_vars = 1, double a = 1.0) + : DyadicMeshSensIntegrator>( + state_vec, adjoint_vec, num_state_vars, a) {} + + /// SA Inviscid two-point (dyadic) entropy conservative flux function + /// \param[in] di - physical coordinate direction in which flux is wanted + /// \param[in] qL - state variables at "left" state + /// \param[in] qR - state variables at "right" state + /// \param[out] flux - fluxes in the direction `di` + /// \note This is simply a wrapper for the function in `euler_fluxes.hpp` + void calcFlux(int di, const mfem::Vector &qL, + const mfem::Vector &qR, mfem::Vector &flux); +}; + +/// Integrator for viscous SA slip-wall boundary condition +/// \tparam dim - number of spatial dimensions (1, 2, or 3) +/// \tparam entvar - if true, states = ent. vars; otherwise, states = conserv. +/// \note This derived class uses the CRTP +// template +// class SASlipWallBCMeshSens : public BoundaryMeshSensIntegrator< +// SASlipWallBCMeshSens> +// { +// public: +// /// Constructs an integrator for a SA slip-wall boundary flux +// /// \param[in] diff_stack - for algorithmic differentiation +// /// \param[in] fe_coll - used to determine the face elements +// /// \param[in] a - used to move residual to lhs (1.0) or rhs(-1.0) +// SASlipWallBCMeshSens(adept::Stack &diff_stack, +// const mfem::GridFunction &state_vec, +// const mfem::GridFunction &adjoint_vec, +// int num_state_vars = 1, double a = 1.0) +// : BoundaryMeshSensIntegrator>( +// state_vec, adjoint_vec, dim + 2, a), stack(diff_stack) {} + +// /// Compute ther derivative of flux_bar*flux w.r.t. `dir` +// /// \param[in] x - location at which the derivative is evaluated (not used) +// /// \param[in] dir - vector normal to the boundary at `x` +// /// \param[in] u - state at which to evaluate the flux +// /// \param[in] flux_bar - flux weighting (e.g. the adjoint) +// /// \param[out] dir_bar - derivative with respect to `dir` +// void calcFluxBar(const mfem::Vector &x, const mfem::Vector &dir, +// const mfem::Vector &u, const mfem::Vector &flux_bar, +// mfem::Vector &dir_bar); + +// protected: +// /// stack used for algorithmic differentiation +// adept::Stack &stack; +// }; + +#include "rans_sens_integ_def.hpp" + +} // namespace mach + +#endif \ No newline at end of file diff --git a/src/physics/fluidflow/rans_sens_integ_def.hpp b/src/physics/fluidflow/rans_sens_integ_def.hpp new file mode 100644 index 00000000..d5afa763 --- /dev/null +++ b/src/physics/fluidflow/rans_sens_integ_def.hpp @@ -0,0 +1,18 @@ +template +void SAInviscidMeshSensIntegrator::calcFlux( + int di, const mfem::Vector &qL, const mfem::Vector &qR, mfem::Vector &flux) +{ + if (entvar) + { + calcIsmailRoeFluxUsingEntVars(di, qL.GetData(), qR.GetData(), + flux.GetData()); + } + else + { + calcIsmailRoeFlux(di, qL.GetData(), qR.GetData(), + flux.GetData()); + } + //add flux term for SA + //flux(dim+2) = 0.5*(qfs(di+1)/qfs(0) + qfs(di+1)/qfs(0))*0.5*(qL(dim+2) + qR(dim+2)); + flux(dim+2) = 0.5*(qL(di+1)/qL(0) + qR(di+1)/qR(0))*0.5*(qL(dim+2) + qR(dim+2)); +} \ No newline at end of file diff --git a/src/physics/fluidflow/viscous_integ_def.hpp b/src/physics/fluidflow/viscous_integ_def.hpp index 2cd8f36c..d5bec70c 100644 --- a/src/physics/fluidflow/viscous_integ_def.hpp +++ b/src/physics/fluidflow/viscous_integ_def.hpp @@ -188,7 +188,6 @@ void SymmetricViscousIntegrator::AssembleElementGrad( jac_term2(sk, sj); } } - } // loop over space dimension d2 } // loop over the element nodes j } // loop over element nodes k diff --git a/src/physics/solver.cpp b/src/physics/solver.cpp index dd279b1d..8084affd 100644 --- a/src/physics/solver.cpp +++ b/src/physics/solver.cpp @@ -1241,11 +1241,22 @@ unique_ptr AbstractSolver::constructLinearSolver( cg->SetPrintLevel(ptl); cg->SetPreconditioner(dynamic_cast(_prec)); } +#ifdef MFEM_USE_SUPERLU + else if (solver_type == "superlu") + { + lin_solver.reset(new SuperLUSolver(comm)); + SuperLUSolver *slu = dynamic_cast(lin_solver.get()); + } +#endif else { throw MachException("Unsupported iterative solver type!\n" "\tavilable options are: HypreGMRES, HypreFGMRES, GMRESSolver,\n" - "\tHyprePCG, CGSolver"); + "\tHyprePCG, CGSolver" +#ifdef MFEM_USE_SUPERLU + ", SuperLU" +#endif + "\n"); } return lin_solver; } @@ -1318,10 +1329,14 @@ unique_ptr AbstractSolver::constructNonlinearSolver( { nonlin_solver.reset(new mfem::NewtonSolver(comm)); } + else if (solver_type == "inexact-newton") + { + nonlin_solver.reset(new InexactNewton(comm)); + } else { throw MachException("Unsupported nonlinear solver type!\n" - "\tavilable options are: newton\n"); + "\tavilable options are: Newton, InexactNewton\n"); } //double eta = 1e-1; //newton_solver.reset(new InexactNewton(comm, eta)); diff --git a/test/unit/CMakeLists.txt b/test/unit/CMakeLists.txt index 14601070..45a13e31 100644 --- a/test/unit/CMakeLists.txt +++ b/test/unit/CMakeLists.txt @@ -31,6 +31,7 @@ set(FLUID_TEST_SRCS test_inexact_newton test_navier_stokes_fluxes test_navier_stokes_integ + #test_navier_stokes_sens_integ test_navier_stokes_assemble test_reconstruction test_sbp_fe @@ -39,6 +40,8 @@ set(FLUID_TEST_SRCS test_utils test_thermal_integ test_res_integ + test_rans_fluxes + test_rans_integ ) # group EM tests diff --git a/test/unit/euler_test_data.cpp b/test/unit/euler_test_data.cpp index bbc59f23..29ffef88 100644 --- a/test/unit/euler_test_data.cpp +++ b/test/unit/euler_test_data.cpp @@ -98,6 +98,31 @@ template void randBaselineVectorPert<1, false>(const Vector &x, Vector &u); template void randBaselineVectorPert<2, false>(const Vector &x, Vector &u); template void randBaselineVectorPert<3, false>(const Vector &x, Vector &u); +template +void randBaselinePertSA(const mfem::Vector &x, mfem::Vector &u) +{ + const double scale = 0.01; + u(0) = rho * (1.0 + scale * uniform_rand(gen)); + u(dim + 1) = rhoe * (1.0 + scale * uniform_rand(gen)); + for (int di = 0; di < dim; ++di) + { + u(di + 1) = rhou[di]* (1.0 + scale * uniform_rand(gen)); + } + if (entvar) + { + mfem::Vector q(u); + mach::calcEntropyVars(q.GetData(), u.GetData()); + } + u(dim + 2) = 3.0 *(1.0 + scale*(uniform_rand(gen) - 0.1)); +} +// explicit instantiation of the templated function above +template void randBaselinePertSA<1, true>(const mfem::Vector &x, mfem::Vector &u); +template void randBaselinePertSA<2, true>(const mfem::Vector &x, mfem::Vector &u); +template void randBaselinePertSA<3, true>(const mfem::Vector &x, mfem::Vector &u); +template void randBaselinePertSA<1, false>(const mfem::Vector &x, mfem::Vector &u); +template void randBaselinePertSA<2, false>(const mfem::Vector &x, mfem::Vector &u); +template void randBaselinePertSA<3, false>(const mfem::Vector &x, mfem::Vector &u); + void randVectorState(const Vector &x, Vector &u) { for (int i = 0; i < u.Size(); ++i) diff --git a/test/unit/euler_test_data.hpp b/test/unit/euler_test_data.hpp index 165fbd96..cbc310a4 100644 --- a/test/unit/euler_test_data.hpp +++ b/test/unit/euler_test_data.hpp @@ -63,6 +63,9 @@ extern double adjJ_data[9]; // Use this for spatial derivatives of entropy-variables extern double delw_data[15]; +// Spalart-Allmaras model parameters +const double sa_params[13] = {0.1355, 0.622, 0.666666666666667, 0.41, 0.3, 2, 7.1, 1.2, 0.5, 10.0, 16.0, 0.7, 0.9}; + /// Returns a perturbed version of the baseline flow state /// \param[in] x - coordinates (not used) /// \param[out] u - pertrubed state variable @@ -71,6 +74,14 @@ extern double delw_data[15]; template void randBaselineVectorPert(const mfem::Vector &x, mfem::Vector &u); +/// Returns a perturbed version of the baseline flow state (Spalart-Allmaras) +/// \param[in] x - coordinates (not used) +/// \param[out] u - pertrubed state variable +/// \tparam dim - number of spatial dimensions (1, 2, or 3) +/// \tparam entvar - if true, returns entropy variables +template +void randBaselinePertSA(const mfem::Vector &x, mfem::Vector &u); + /// Returns a random state with entries uniformly distributed in [-1,1] /// \param[in] x - coordinates (not used) /// \param[out] u - rand state variable diff --git a/test/unit/test_navier_stokes_sens_integ.cpp b/test/unit/test_navier_stokes_sens_integ.cpp new file mode 100644 index 00000000..8aaf35fd --- /dev/null +++ b/test/unit/test_navier_stokes_sens_integ.cpp @@ -0,0 +1,90 @@ +#include + +#include "catch.hpp" +#include "mfem.hpp" +#include "euler_fluxes.hpp" +#include "euler_integ.hpp" +#include "euler_sens_integ.hpp" +#include "euler_test_data.hpp" +#include "navier_stokes_integ.hpp" +#include "navier_stokes_sens_integ.hpp" + +TEMPLATE_TEST_CASE_SIG("ESViscousMeshSensIntegrator::AssembleElementVector", + "[ESViscousMeshSens]", + ((bool entvar), entvar), false) +{ + using namespace mfem; + using namespace euler_data; + + const int dim = 2; // templating is hard here because mesh constructors + double Re_num = 1; + double Pr_num = 1; + double mu = -1; + int num_state = dim + 2; + adept::Stack diff_stack; // needed for forward problem + double delta = 1e-5; + + // generate a 8 element mesh + int num_edge = 2; + std::unique_ptr mesh(new Mesh(num_edge, num_edge, Element::TRIANGLE, + true /* gen. edges */, 1.0, 1.0, true)); + mesh->EnsureNodes(); + for (int p = 1; p <= 4; ++p) + { + DYNAMIC_SECTION("...for degree p = " << p) + { + // get the finite-element space for the state and adjoint + std::unique_ptr fec( + new SBPCollection(p, dim)); + std::unique_ptr fes(new FiniteElementSpace( + mesh.get(), fec.get(), num_state, Ordering::byVDIM)); + + // we use res for finite-difference approximation + NonlinearForm res(fes.get()); + res.AddDomainIntegrator( + new mach::ESViscousIntegrator<2>(diff_stack, Re_num, Pr_num, mu)); + + // initialize state and adjoint; here we randomly perturb a constant state + GridFunction state(fes.get()), adjoint(fes.get()); + VectorFunctionCoefficient pert(num_state, randBaselineVectorPert<2,entvar>); + state.ProjectCoefficient(pert); + adjoint.ProjectCoefficient(pert); + + // extract mesh nodes and get their finite-element space + GridFunction *x_nodes = mesh->GetNodes(); + FiniteElementSpace *mesh_fes = x_nodes->FESpace(); + + // build the linear form for d(psi^T R)/dx + LinearForm dfdx(mesh_fes); + dfdx.AddDomainIntegrator( + new mach::ESViscousMeshSensIntegrator<2,entvar>( + state, adjoint, Re_num, Pr_num, mu)); + dfdx.Assemble(); + + // initialize the vector that we use to perturb the mesh nodes + GridFunction v(mesh_fes); + VectorFunctionCoefficient v_rand(dim, randVectorState); + v.ProjectCoefficient(v_rand); + + // contract dfdx with v + double dfdx_v = dfdx * v; + + // now compute the finite-difference approximation... + GridFunction x_pert(*x_nodes); + GridFunction r(fes.get()); + x_pert.Add(delta, v); + mesh->SetNodes(x_pert); + res.Mult(state, r); + double dfdx_v_fd = adjoint * r; + x_pert.Add(-2 * delta, v); + mesh->SetNodes(x_pert); + res.Mult(state, r); + dfdx_v_fd -= adjoint * r; + dfdx_v_fd /= (2 * delta); + mesh->SetNodes(*x_nodes); // remember to reset the mesh nodes + + REQUIRE(dfdx_v == Approx(dfdx_v_fd).margin(1e-10)); + } + } +} + diff --git a/test/unit/test_rans_fluxes.cpp b/test/unit/test_rans_fluxes.cpp new file mode 100644 index 00000000..a339da06 --- /dev/null +++ b/test/unit/test_rans_fluxes.cpp @@ -0,0 +1,663 @@ +#include + +#include "catch.hpp" +#include "mfem.hpp" +#include "euler_integ.hpp" +#include "rans_fluxes.hpp" +#include "rans_integ.hpp" +#include "euler_test_data.hpp" + +//rotational vortex state +void uvort(const mfem::Vector &x, mfem::Vector& u); + +//exact vorticity magnitude field +double s_proj(const mfem::Vector &x); + +TEMPLATE_TEST_CASE_SIG("SA calcVorticity Accuracy", "[SAVorticity]", + ((bool entvar), entvar), false) +{ + using namespace mfem; + using namespace euler_data; + + const int dim = 2; // templating is hard here because mesh constructors + int num_state = dim + 3; + + // generate a mesh + int num_edge = 1; + std::unique_ptr mesh(new Mesh(num_edge, num_edge, Element::TRIANGLE, + true /* gen. edges */, 1.0, 1.0, true)); + + for (int p = 1; p <= 1; ++p) + { + DYNAMIC_SECTION("Vorticity correct for degree p = " << p) + { + std::unique_ptr fec( + new SBPCollection(p, dim)); + std::unique_ptr fes(new FiniteElementSpace( + mesh.get(), fec.get(), num_state, Ordering::byVDIM)); + std::unique_ptr fess(new FiniteElementSpace( + mesh.get(), fec.get())); + + // state gridfunction + GridFunction u(fes.get()); + + // vorticity magnitude gridfunctions + GridFunction s_exact(fess.get()); + GridFunction s_comp(fess.get()); + + FunctionCoefficient s_coeff(*s_proj); + VectorFunctionCoefficient u_coeff(dim+3, *uvort); + s_exact.ProjectCoefficient(s_coeff); + u.ProjectCoefficient(u_coeff); + + Array vdofs, vdofs2; + Vector el_x, el_y; + const FiniteElement *fe; + ElementTransformation *T; + + // element loop + for (int i = 0; i < fes->GetNE(); i++) + { + + fe = fes->GetFE(i); + fes->GetElementVDofs(i, vdofs); + fess->GetElementVDofs(i, vdofs2); + T = fes->GetElementTransformation(i); + const SBPFiniteElement &sbp = dynamic_cast(*fe); + int num_nodes = sbp.GetDof(); + int num_states = dim+3; + u.GetSubVector(vdofs, el_x); + s_comp.GetSubVector(vdofs2, el_y); + + Vector curl(3); + Vector Duidi(num_states); + DenseMatrix Dui(num_states, dim); + DenseMatrix uc(el_x.GetData(), num_nodes, num_states); + DenseMatrix ucc; ucc = uc; + ucc.Transpose(); + + // dof loop + for (int j = 0; j < num_nodes; j++) + { + const IntegrationPoint &node = fe->GetNodes().IntPoint(j); + T->SetIntPoint(&node); + + Duidi = 0.0; Dui = 0.0; + for (int di = 0; di < dim; di++) + { + sbp.multStrongOperator(di, j, ucc, Duidi); + Dui.SetCol(di, Duidi); + } + + // get vorticity at dof + mach::calcVorticity(Dui.GetData(), T->InverseJacobian().GetData(), + curl.GetData()); + + double S = curl.Norml2(); + + el_y(j) = S; + } + // insert into vorticity gridfunction + s_comp.AddElementVector(vdofs2, el_y); + + } + + ofstream sol_ofs("why.vtk"); + sol_ofs.precision(14); + mesh->PrintVTK(sol_ofs, 1); + s_comp.SaveVTK(sol_ofs, "vorticity_error", 1); + + // compare vorticity + for (int j = 0; j < s_comp.Size(); j++) + { + std::cout << "Computed: " << s_comp(j) << std::endl; + std::cout << "Exact: " << s_exact(j) << std::endl; + //REQUIRE(s_comp(j) == Approx(s_exact(j)).margin(1e-10)); + } + } + } +} + + +TEMPLATE_TEST_CASE_SIG("SA calcVorticity Jacobian", "[SAVorticity]", + ((int dim), dim), 2, 3) +{ + using namespace euler_data; + mfem::DenseMatrix Dw(delw_data, dim+3, dim); + mfem::DenseMatrix trans_jac(adjJ_data, dim, dim); + mfem::Vector curl_plus(3); + mfem::Vector curl_minus(3); + mfem::Vector v(dim+3); + mfem::Vector jac_v(3); + double delta = 1e-5; + + for (int di = 0; di < dim; ++di) + { + Dw(dim+2, di) = 0.7; + } + // create perturbation vector + for (int di = 0; di < dim + 3; ++di) + { + v(di) = vec_pert[di]; + } + adept::Stack stack; + + DYNAMIC_SECTION("Jacobians w.r.t Dw is correct") + { + std::vector mat_vec_jac(dim); + for (int d = 0; d < dim; ++d) + { + mat_vec_jac[d].SetSize(3, dim+3); + } + // compute the jacobian + mach::calcVorticityJacDw(stack, Dw.GetData(), trans_jac.GetData(), + mat_vec_jac); + + for (int d = 0; d < dim; ++d) + { + // perturb one column of delw everytime + mfem::DenseMatrix Dw_plus(Dw), Dw_minus(Dw); + for (int s = 0; s < dim+3; ++s) + { + Dw_plus.GetColumn(d)[s] += v(s) * delta; + Dw_minus.GetColumn(d)[s] -= v(s) * delta; + } + // get perturbed states conv vector + mach::calcVorticity(Dw_plus.GetData(), trans_jac.GetData(), + curl_plus.GetData()); + mach::calcVorticity(Dw_minus.GetData(), trans_jac.GetData(), + curl_minus.GetData()); + + mat_vec_jac[d].Mult(v, jac_v); + + // finite difference jacobian + mfem::Vector jac_v_fd(3); + jac_v_fd = curl_plus; + jac_v_fd -= curl_minus; + jac_v_fd /= 2.0 * delta; + //std::cout << "viscous applyScaling jac Dw:" << std::endl; + for (int j = 0; j < 3; ++j) + { + std::cout << "FD " << j << " Deriv: " << jac_v_fd[j] << std::endl; + std::cout << "AN " << j << " Deriv: " << jac_v(j) << std::endl; + REQUIRE(jac_v(j) == Approx(jac_v_fd[j]).margin(1e-10)); + } + } + } +} + +TEMPLATE_TEST_CASE_SIG("SA calcGrad Jacobian", "[SAVorticity]", + ((int dim), dim), 2, 3) +{ + using namespace euler_data; + mfem::DenseMatrix Dw(delw_data, dim+3, dim); + mfem::DenseMatrix trans_jac(adjJ_data, dim, dim); + mfem::Vector grad_plus(dim); + mfem::Vector grad_minus(dim); + mfem::Vector v(dim+3); + mfem::Vector jac_v(dim); + double delta = 1e-5; + + for (int di = 0; di < dim; ++di) + { + Dw(dim+2, di) = 0.7; + } + // create perturbation vector + for (int di = 0; di < dim + 3; ++di) + { + v(di) = vec_pert[di]; + } + adept::Stack stack; + + DYNAMIC_SECTION("Jacobians w.r.t Dw is correct") + { + std::vector mat_vec_jac(dim); + for (int d = 0; d < dim; ++d) + { + mat_vec_jac[d].SetSize(dim, dim+3); + } + // compute the jacobian + mach::calcGradJacDw(stack, dim+2, Dw.GetData(), trans_jac.GetData(), + mat_vec_jac); + + for (int d = 0; d < dim; ++d) + { + // perturb one column of delw everytime + mfem::DenseMatrix Dw_plus(Dw), Dw_minus(Dw); + for (int s = 0; s < dim+3; ++s) + { + Dw_plus.GetColumn(d)[s] += v(s) * delta; + Dw_minus.GetColumn(d)[s] -= v(s) * delta; + } + // get perturbed states conv vector + mach::calcGrad(dim+2, Dw_plus.GetData(), trans_jac.GetData(), + grad_plus.GetData()); + mach::calcGrad(dim+2, Dw_minus.GetData(), trans_jac.GetData(), + grad_minus.GetData()); + + mat_vec_jac[d].Mult(v, jac_v); + + // finite difference jacobian + mfem::Vector jac_v_fd(dim); + jac_v_fd = grad_plus; + jac_v_fd -= grad_minus; + jac_v_fd /= 2.0 * delta; + //std::cout << "viscous applyScaling jac Dw:" << std::endl; + for (int j = 0; j < dim; ++j) + { + std::cout << "FD " << j << " Deriv: " << jac_v_fd[j] << std::endl; + std::cout << "AN " << j << " Deriv: " << jac_v(j) << std::endl; + REQUIRE(jac_v(j) == Approx(jac_v_fd[j]).margin(1e-10)); + } + } + } +} + + +TEST_CASE("SA S derivatives", "[SASource]") +{ + using namespace euler_data; + const int dim = 2; + mfem::Vector q(dim+3); + mfem::Vector grad(dim); + double delta = 1e-5; + + grad(0) = 0.0; + grad(1) = 0.0; + + q(0) = rho; + q(dim + 1) = rhoe; + q(dim + 2) = 4; + for (int di = 0; di < dim; ++di) + { + q(di + 1) = rhou[di]; + } + + mfem::Vector sacs(13); + // create SA parameter vector + for (int di = 0; di < 13; ++di) + { + sacs(di) = sa_params[di]; + } + + adept::Stack stack; + std::vector sacs_a(sacs.Size()); + adept::set_values(sacs_a.data(), sacs.Size(), sacs.GetData()); + std::vector q_a(sacs.Size()); + adept::set_values(q_a.data(), q.Size(), q.GetData()); + std::vector grad_a(grad.Size()); + adept::set_values(grad_a.data(), grad.Size(), grad.GetData()); + + double sval = 0.0; + + DYNAMIC_SECTION("Jacobians w.r.t S is correct") + { + // compute the jacobian + stack.new_recording(); + mfem::Vector dSrcdu(1); + adouble S_a; + S_a.set_value(sval); + adouble d_a = 1.001; + adouble mu_a = 1.1; + adouble Re_a = 5000; + adouble src;// = mach::calcSADestruction(q_a.data(), S_a, mu_a, d_a, Re_a, sacs_a.data()); + src = mach::calcSASource( + q_a.data(), grad_a.data(), sacs_a.data())/Re_a; + // use negative model if turbulent viscosity is negative + // if (q_a[dim+2] < 0) + // { + // src += mach::calcSANegativeProduction( + // q_a.data(), S_a, sacs_a.data()); + // src += mach::calcSANegativeDestruction( + // q_a.data(), d_a, sacs_a.data())/Re_a; + // } + // else + { + src += mach::calcSAProduction( + q_a.data(), mu_a, d_a, S_a, Re_a, sacs_a.data())/Re_a; + src += mach::calcSADestruction( + q_a.data(), mu_a, d_a, S_a, Re_a, sacs_a.data())/Re_a; + } + + stack.independent(S_a); + stack.dependent(src); + stack.jacobian(dSrcdu.GetData()); + + double Sp = sval + delta; double Sm = sval - delta; + double d = 1.001; + double mu = 1.1; + double Re = 5000; + + // get perturbed states conv vector + //double srcp = mach::calcSADestruction(q.GetData(), Sp, mu, d, Re, sacs.GetData()); + //double srcm = mach::calcSADestruction(q.GetData(), Sm, mu, d, Re, sacs.GetData()); + + double srcp = mach::calcSASource( + q.GetData(), grad.GetData(), sacs.GetData())/Re; + // if (q(dim+2) < 0) + // { + // srcp += mach::calcSANegativeProduction( + // q.GetData(), Sp, sacs.GetData()); + // srcp += mach::calcSANegativeDestruction( + // q.GetData(), d, sacs.GetData())/Re; + // } + // else + { + srcp += mach::calcSAProduction( + q.GetData(), mu, d, Sp, Re, sacs.GetData())/Re; + srcp += mach::calcSADestruction( + q.GetData(), mu, d, Sp, Re, sacs.GetData())/Re; + } + + double srcm = mach::calcSASource( + q.GetData(), grad.GetData(), sacs.GetData())/Re; + // if (q(dim+2) < 0) + // { + // srcm += mach::calcSANegativeProduction( + // q.GetData(), Sm, sacs.GetData()); + // srcm += mach::calcSANegativeDestruction( + // q.GetData(), d, sacs.GetData())/Re; + // } + // else + { + srcm += mach::calcSAProduction( + q.GetData(), mu, d, Sm, Re, sacs.GetData())/Re; + srcm += mach::calcSADestruction( + q.GetData(), mu, d, Sm, Re, sacs.GetData())/Re; + } + + // finite difference jacobian + double jac_fd = srcp; + jac_fd -= srcm; + jac_fd /= 2.0 * delta; + + std::cout << "FD Deriv: " << jac_fd << std::endl; + std::cout << "AN Deriv: " << dSrcdu(0) << std::endl; + REQUIRE(dSrcdu(0) == Approx(jac_fd).margin(1e-10)); + + } +} + +TEST_CASE("SA grad derivatives", "[SAVorticity]") +{ + using namespace euler_data; + const int dim = 2; + mfem::Vector q(dim+3); + mfem::Vector grad(dim); + mfem::Vector v(dim); + double delta = 1e-5; + + grad(0) = 0.0; + grad(1) = 0.0; + + q(0) = rho; + q(dim + 1) = rhoe; + q(dim + 2) = 4; + for (int di = 0; di < dim; ++di) + { + q(di + 1) = rhou[di]; + } + + for (int di = 0; di < dim; ++di) + { + v(di) = vec_pert[di]; + } + + mfem::Vector sacs(13); + // create SA parameter vector + for (int di = 0; di < 13; ++di) + { + sacs(di) = sa_params[di]; + } + + adept::Stack stack; + DYNAMIC_SECTION("Jacobians w.r.t grad is correct") + { + // compute the jacobian + stack.new_recording(); + mfem::Vector dSrcdgrad(dim); + std::vector sacs_a(sacs.Size()); + adept::set_values(sacs_a.data(), sacs.Size(), sacs.GetData()); + std::vector q_a(sacs.Size()); + adept::set_values(q_a.data(), q.Size(), q.GetData()); + std::vector grad_a(grad.Size()); + adept::set_values(grad_a.data(), grad.Size(), grad.GetData()); + adouble S_a = 0.000; + adouble d_a = 10.1; + adouble mu_a = 1.1; + adouble Re_a = 5000000; + adouble src;// = mach::calcSADestruction(q_a.data(), S_a, mu_a, d_a, Re_a, sacs_a.data()); + src = mach::calcSASource( + q_a.data(), grad_a.data(), sacs_a.data())/Re_a; + // use negative model if turbulent viscosity is negative + if (q_a[dim+2] < 0) + { + src += mach::calcSANegativeProduction( + q_a.data(), S_a, sacs_a.data()); + src += mach::calcSANegativeDestruction( + q_a.data(), d_a, Re_a, sacs_a.data()); + } + else + { + src += mach::calcSAProduction( + q_a.data(), mu_a, d_a, S_a, Re_a, sacs_a.data()); + src += mach::calcSADestruction( + q_a.data(), mu_a, d_a, S_a, Re_a, sacs_a.data()); + } + + stack.independent(grad_a.data(), dim); + stack.dependent(src); + stack.jacobian(dSrcdgrad.GetData()); + + double S = 0.000; + double d = 10.1; + double mu = 1.1; + double Re = 5000000; + mfem::Vector gradp(grad); mfem::Vector gradm(grad); + gradp.Add(delta, v); gradm.Add(-delta, v); + + // get perturbed states conv vector + //double srcp = mach::calcSADestruction(q.GetData(), Sp, mu, d, Re, sacs.GetData()); + //double srcm = mach::calcSADestruction(q.GetData(), Sm, mu, d, Re, sacs.GetData()); + + double srcp = mach::calcSASource( + q.GetData(), gradp.GetData(), sacs.GetData())/Re; + if (q(dim+2) < 0) + { + srcp += mach::calcSANegativeProduction( + q.GetData(), S, sacs.GetData()); + srcp += mach::calcSANegativeDestruction( + q.GetData(), d, Re, sacs.GetData()); + } + else + { + srcp += mach::calcSAProduction( + q.GetData(), mu, d, S, Re, sacs.GetData()); + srcp += mach::calcSADestruction( + q.GetData(), mu, d, S, Re, sacs.GetData()); + } + + double srcm = mach::calcSASource( + q.GetData(), gradm.GetData(), sacs.GetData())/Re; + if (q(dim+2) < 0) + { + srcm += mach::calcSANegativeProduction( + q.GetData(), S, sacs.GetData()); + srcm += mach::calcSANegativeDestruction( + q.GetData(), d, Re, sacs.GetData()); + } + else + { + srcm += mach::calcSAProduction( + q.GetData(), mu, d, S, Re, sacs.GetData()); + srcm += mach::calcSADestruction( + q.GetData(), mu, d, S, Re, sacs.GetData()); + } + + // finite difference jacobian + double jac_fd = srcp; + jac_fd -= srcm; + jac_fd /= 2.0 * delta; + + double jac_v = dSrcdgrad*v; + + std::cout << "FD Deriv: " << jac_fd << std::endl; + std::cout << "AN Deriv: " << jac_v << std::endl; + REQUIRE(jac_v == Approx(jac_fd).margin(1e-10)); + + } +} + +TEST_CASE("SA q derivatives", "[SAVorticity]") +{ + using namespace euler_data; + const int dim = 2; + mfem::Vector q(dim+3); + mfem::Vector grad(dim); + mfem::Vector v(dim+3); + double delta = 1e-5; + + grad(0) = 0.0; + grad(1) = 0.0; + + q(0) = rho; + q(dim + 1) = rhoe; + q(dim + 2) = 4; + for (int di = 0; di < dim; ++di) + { + q(di + 1) = rhou[di]; + } + + for (int di = 0; di < dim+3; ++di) + { + v(di) = vec_pert[di]; + } + + mfem::Vector sacs(13); + // create SA parameter vector + for (int di = 0; di < 13; ++di) + { + sacs(di) = sa_params[di]; + } + + adept::Stack stack; + DYNAMIC_SECTION("Jacobians w.r.t grad is correct") + { + // compute the jacobian + stack.new_recording(); + mfem::Vector dSrcdq(dim+3); + std::vector sacs_a(sacs.Size()); + adept::set_values(sacs_a.data(), sacs.Size(), sacs.GetData()); + std::vector q_a(sacs.Size()); + adept::set_values(q_a.data(), q.Size(), q.GetData()); + std::vector grad_a(grad.Size()); + adept::set_values(grad_a.data(), grad.Size(), grad.GetData()); + adouble S_a = 0.000; + adouble d_a = 10.1; + adouble mu_a = 1.1; + adouble Re_a = 5000000; + adouble src;// = mach::calcSADestruction(q_a.data(), S_a, mu_a, d_a, Re_a, sacs_a.data()); + src = mach::calcSASource( + q_a.data(), grad_a.data(), sacs_a.data())/Re_a; + // use negative model if turbulent viscosity is negative + if (q_a[dim+2] < 0) + { + src += mach::calcSANegativeProduction( + q_a.data(), S_a, sacs_a.data()); + src += mach::calcSANegativeDestruction( + q_a.data(), d_a, Re_a, sacs_a.data()); + } + else + { + src += mach::calcSAProduction( + q_a.data(), mu_a, d_a, S_a, Re_a, sacs_a.data()); + src += mach::calcSADestruction( + q_a.data(), mu_a, d_a, S_a, Re_a, sacs_a.data()); + } + + stack.independent(q_a.data(), dim+3); + stack.dependent(src); + stack.jacobian(dSrcdq.GetData()); + + double S = 0.000; + double d = 10.1; + double mu = 1.1; + double Re = 5000000; + mfem::Vector qp(q); mfem::Vector qm(q); + qp.Add(delta, v); qm.Add(-delta, v); + + // get perturbed states conv vector + //double srcp = mach::calcSADestruction(q.GetData(), Sp, mu, d, Re, sacs.GetData()); + //double srcm = mach::calcSADestruction(q.GetData(), Sm, mu, d, Re, sacs.GetData()); + + double srcp = mach::calcSASource( + qp.GetData(), grad.GetData(), sacs.GetData())/Re; + if (qp(dim+2) < 0) + { + srcp += mach::calcSANegativeProduction( + qp.GetData(), S, sacs.GetData()); + srcp += mach::calcSANegativeDestruction( + qp.GetData(), d, Re, sacs.GetData()); + } + else + { + srcp += mach::calcSAProduction( + qp.GetData(), mu, d, S, Re, sacs.GetData()); + srcp += mach::calcSADestruction( + qp.GetData(), mu, d, S, Re, sacs.GetData()); + } + + double srcm = mach::calcSASource( + qm.GetData(), grad.GetData(), sacs.GetData())/Re; + if (qm(dim+2) < 0) + { + srcm += mach::calcSANegativeProduction( + qm.GetData(), S, sacs.GetData()); + srcm += mach::calcSANegativeDestruction( + qm.GetData(), d, Re, sacs.GetData()); + } + else + { + srcm += mach::calcSAProduction( + qm.GetData(), mu, d, S, Re, sacs.GetData()); + srcm += mach::calcSADestruction( + qm.GetData(), mu, d, S, Re, sacs.GetData()); + } + + // finite difference jacobian + double jac_fd = srcp; + jac_fd -= srcm; + jac_fd /= 2.0 * delta; + + double jac_v = dSrcdq*v; + + std::cout << "FD Deriv: " << jac_fd << std::endl; + std::cout << "AN Deriv: " << jac_v << std::endl; + REQUIRE(jac_v == Approx(jac_fd).margin(1e-10)); + + } +} + +void uvort(const mfem::Vector &x, mfem::Vector& q) +{ + // q.SetSize(4); + // Vector u(4); + q.SetSize(5); + mfem::Vector u(5); + double om = 4; + u = 0.0; + u(0) = 1.0; + u(1) = -om*x(1); + u(2) = om*x(0); + u(3) = 1.0/(1.4) + 0.5*1.5*1.5; + u(4) = u(0)*3; + + q = u; +} + +double s_proj(const mfem::Vector &x) +{ + double om = 4; + + return 2*om; +} \ No newline at end of file diff --git a/test/unit/test_rans_integ.cpp b/test/unit/test_rans_integ.cpp new file mode 100644 index 00000000..d9e7822c --- /dev/null +++ b/test/unit/test_rans_integ.cpp @@ -0,0 +1,1375 @@ +#include +#include +#include "catch.hpp" +#include "mfem.hpp" +#include "euler_integ.hpp" +#include "rans_fluxes.hpp" +#include "rans_integ.hpp" +#include "euler_test_data.hpp" + +void uinit(const mfem::Vector &x, mfem::Vector& u); + +TEMPLATE_TEST_CASE_SIG("SA Inviscid Flux Test", "[sa_inviscid_flux_test]", + ((int dim), dim), 1, 2, 3) +{ + using namespace euler_data; + mfem::Vector qL(dim + 3); + mfem::Vector qR(dim + 3); + mfem::Vector flux1(dim + 2); + mfem::Vector flux2(dim + 3); + mfem::Vector nrm(dim); + qL(0) = rho; qR(0) = rho*1.1; + qL(dim + 1) = rhoe; qR(dim + 1) = rhoe*1.1; + for (int di = 0; di < dim; ++di) + { + qL(di + 1) = rhou[di]; + qR(di + 1) = rhou[di]*1.1; + nrm(di) = dir[di]; + } + qL(dim + 2) = 5.0; qR(dim + 2) = 5.5; + adept::Stack stack; + mach::IsmailRoeIntegrator irinteg(stack); + mach::SAInviscidIntegrator sainteg(stack); + + SECTION("Check if SA integrator flux function matches the output for the conservative variables") + { + // calculate the flux + irinteg.calcFlux(0, qL, qR, flux1); + sainteg.calcFlux(0, qL, qR, flux2); + + // check if euler variables are the same as before + for(int n = 0; n < dim+2; n++) + { + // std::cout << "Ismail-Roe " << n << " Flux: " << flux1(n) << std::endl; + // std::cout << "Spalart-Allmaras " << n << " Flux: " << flux2(n) << std::endl; + REQUIRE(flux1(n) - flux2(n) == Approx(0.0)); + } + } +} + +TEMPLATE_TEST_CASE_SIG("SAInviscid Jacobian", "[SAInviscid]", + ((int dim), dim), 1, 2, 3) +{ + using namespace euler_data; + // copy the data into mfem vectors for convenience + mfem::Vector qL(dim + 3); + mfem::Vector qR(dim + 3); + mfem::Vector flux(dim + 3); + mfem::Vector flux_plus(dim + 3); + mfem::Vector flux_minus(dim + 3); + mfem::Vector v(dim + 3); + mfem::Vector jac_v(dim + 3); + mfem::DenseMatrix jacL(dim + 3, 2 * (dim + 3)); + mfem::DenseMatrix jacR(dim + 3, 2 * (dim + 3)); + double delta = 1e-5; + qL(0) = rho; + qL(dim + 1) = rhoe; + qL(dim + 2) = 4; + qR(0) = rho2; + qR(dim + 1) = rhoe2; + qR(dim + 2) = 4.5; + for (int di = 0; di < dim; ++di) + { + qL(di + 1) = rhou[di]; + qR(di + 1) = rhou2[di]; + } + // create perturbation vector + for (int di = 0; di < dim + 3; ++di) + { + v(di) = vec_pert[di]; + } + // perturbed vectors + mfem::Vector qL_plus(qL), qL_minus(qL); + mfem::Vector qR_plus(qR), qR_minus(qR); + adept::Stack diff_stack; + mach::SAInviscidIntegrator sainvinteg(diff_stack); + // +ve perturbation + qL_plus.Add(delta, v); + qR_plus.Add(delta, v); + // -ve perturbation + qL_minus.Add(-delta, v); + qR_minus.Add(-delta, v); + for (int di = 0; di < dim; ++di) + { + DYNAMIC_SECTION("Ismail-Roe flux jacsainvintegian is correct w.r.t left state ") + { + // get perturbed states flux vector + sainvinteg.calcFlux(di, qL_plus, qR, flux_plus); + sainvinteg.calcFlux(di, qL_minus, qR, flux_minus); + // compute the jacobian + sainvinteg.calcFluxJacStates(di, qL, qR, jacL, jacR); + jacL.Mult(v, jac_v); + // finite difference jacobian + mfem::Vector jac_v_fd(flux_plus); + jac_v_fd -= flux_minus; + jac_v_fd /= 2.0 * delta; + // compare each component of the matrix-vector products + for (int i = 0; i < dim + 2; ++i) + { + REQUIRE(jac_v[i] == Approx(jac_v_fd[i]).margin(1e-10)); + } + } + } +} + +TEMPLATE_TEST_CASE_SIG("SAFarFieldBC Jacobian", "[SAFarFieldBC]", + ((int dim), dim), 1, 2, 3) +{ + using namespace euler_data; + // copy the data into mfem vectors for convenience + mfem::Vector nrm(dim); mfem::Vector x(dim); + mfem::Vector q(dim + 3); + mfem::Vector qfar(dim + 3); + mfem::Vector flux(dim + 3); + mfem::Vector v(dim + 3); + mfem::Vector jac_v(dim + 3); + mfem::DenseMatrix Dw; + mfem::DenseMatrix jac(dim + 3, dim + 3); + double delta = 1e-5; + + x = 0.0; + q(0) = rho; + q(dim + 1) = rhoe; + q(dim + 2) = 4; + for (int di = 0; di < dim; ++di) + { + q(di + 1) = rhou[di]; + } + qfar.Set(1.1, q); + // create direction vector + for (int di = 0; di < dim; ++di) + { + nrm(di) = dir[di]; + } + // create perturbation vector + for (int di = 0; di < dim + 3; ++di) + { + v(di) = vec_pert[di]; + } + // perturbed vectors + mfem::Vector q_plus(q), q_minus(q); + mfem::Vector flux_plus(q), flux_minus(q); + adept::Stack diff_stack; + mfem::H1_FECollection fe_coll(1); //dummy + double Re = 1000.0; double Pr = 1.0; + mach::SAFarFieldBC safarfieldinteg(diff_stack, &fe_coll, Re, Pr, qfar, 1.0); + // +ve perturbation + q_plus.Add(delta, v); + // -ve perturbation + q_minus.Add(-delta, v); + for (int di = 0; di < 2; di++) //reverse direction to check both inflow and outflow + { + if(di == 1) + nrm *= -1.0; + + DYNAMIC_SECTION("SA Far-Field BC jacobian is correct w.r.t state when di = "< sanoslipinteg(diff_stack, &fe_coll, Re, Pr, sacs, qfar, 1.5); + // +ve perturbation + q_plus.Add(delta, v); + // -ve perturbation + q_minus.Add(-delta, v); + DYNAMIC_SECTION("SA No-Slip BC jacobian is correct w.r.t state") + { + + // get perturbed states flux vector + sanoslipinteg.calcFlux(x, nrm, 1.0, q_plus, Dw, flux_plus); + sanoslipinteg.calcFlux(x, nrm, 1.0, q_minus, Dw, flux_minus); + // compute the jacobian + sanoslipinteg.calcFluxJacState(x, nrm, 1.0, q, Dw, jac); + jac.Mult(v, jac_v); + // finite difference jacobian + mfem::Vector jac_v_fd(flux_plus); + jac_v_fd -= flux_minus; + jac_v_fd /= 2.0 * delta; + // compare each component of the matrix-vector products + for (int i = 0; i < dim + 3; ++i) + { + // std::cout << "FD " << i << " Deriv: " << jac_v_fd[i] << std::endl; + // std::cout << "AN " << i << " Deriv: " << jac_v[i] << std::endl; + REQUIRE(jac_v[i] == Approx(jac_v_fd[i]).margin(1e-10)); + } + } +} + +TEMPLATE_TEST_CASE_SIG("SANoSlipAdiabaticWallBC Dw Jacobian", "[SANoSlipAdiabaticWallBC]", + ((int dim), dim), 1, 2, 3) +{ + using namespace euler_data; + // copy the data into mfem vectors for convenience + mfem::Vector nrm(dim); mfem::Vector x(dim); + mfem::Vector q(dim + 3); + mfem::Vector qfar(dim + 3); + mfem::Vector flux(dim + 3); + mfem::Vector v(dim + 3); + mfem::Vector jac_v(dim + 3); + mfem::DenseMatrix Dw(delw_data, dim+3, dim); + mfem::DenseMatrix jac(dim + 3, dim + 3); + mfem::Vector sacs(13); + // create SA parameter vector + for (int di = 0; di < 13; ++di) + { + sacs(di) = sa_params[di]; + } + double delta = 1e-5; + + x = 0.0; + q(0) = rho; + q(dim + 1) = rhoe; + q(dim + 2) = 4; + for (int di = 0; di < dim; ++di) + { + q(di + 1) = rhou[di]; + Dw(dim+2, di) = 0.7; + } + qfar.Set(1.1, q); + // create direction vector + for (int di = 0; di < dim; ++di) + { + nrm(di) = dir[di]; + } + // create perturbation vector + for (int di = 0; di < dim + 3; ++di) + { + v(di) = vec_pert[di]; + } + // perturbed vectors + mfem::Vector flux_plus(q), flux_minus(q); + adept::Stack diff_stack; + mfem::H1_FECollection fe_coll(1); //dummy + double Re = 1000.0; double Pr = 1.0; + mach::SANoSlipAdiabaticWallBC sanoslipinteg(diff_stack, &fe_coll, Re, Pr, sacs, qfar, 1.5); + DYNAMIC_SECTION("SA No-Slip BC Dw jacobian is correct w.r.t state") + { + std::vector mat_vec_jac(dim); + for (int d = 0; d < dim; ++d) + { + mat_vec_jac[d].SetSize(dim+3); + } + // compute the jacobian + sanoslipinteg.calcFluxJacDw(x, nrm, 1.0, q, Dw, mat_vec_jac); + + for (int d = 0; d < dim; ++d) + { + // perturb one column of delw everytime + mfem::DenseMatrix Dw_plus(Dw), Dw_minus(Dw); + for (int s = 0; s < dim+3; ++s) + { + Dw_plus.GetColumn(d)[s] += v(s) * delta; + Dw_minus.GetColumn(d)[s] -= v(s) * delta; + } + // get perturbed states flux vector + sanoslipinteg.calcFlux(x, nrm, 1.0, q, Dw_plus, flux_plus); + sanoslipinteg.calcFlux(x, nrm, 1.0, q, Dw_minus, flux_minus); + + mat_vec_jac[d].Mult(v, jac_v); + // finite difference jacobian + mfem::Vector jac_v_fd(flux_plus); + jac_v_fd -= flux_minus; + jac_v_fd /= 2.0 * delta; + // compare each component of the matrix-vector products + for (int i = 0; i < dim + 3; ++i) + { + // std::cout << "FD " << i << " Deriv: " << jac_v_fd[i] << std::endl; + // std::cout << "AN " << i << " Deriv: " << jac_v[i] << std::endl; + REQUIRE(jac_v[i] == Approx(jac_v_fd[i]).margin(1e-10)); + } + } + } +} + +TEMPLATE_TEST_CASE_SIG("SAViscousSlipWallBC Jacobian", "[SASlipWallBC]", + ((int dim), dim), 1, 2) +{ + using namespace euler_data; + // copy the data into mfem vectors for convenience +// std::cout << "Slip-Wall State Jac" << std::endl; + mfem::Vector nrm(dim); mfem::Vector x(dim); + mfem::Vector q(dim + 3); + mfem::Vector flux(dim + 3); + mfem::Vector v(dim + 3); + mfem::Vector jac_v(dim + 3); + mfem::DenseMatrix Dw(delw_data, dim+3, dim); + mfem::DenseMatrix jac(dim + 3, dim + 3); + mfem::Vector sacs(13); + // create SA parameter vector + for (int di = 0; di < 13; ++di) + { + sacs(di) = sa_params[di]; + } + double delta = 1e-5; + + x = 0.0; + q(0) = rho; + q(dim + 1) = rhoe; + q(dim + 2) = 4; + for (int di = 0; di < dim; ++di) + { + q(di + 1) = rhou[di]; + Dw(dim+2, di) = 0.7+0.1*di; + } + // create direction vector + for (int di = 0; di < dim; ++di) + { + nrm(di) = dir[di]; + } + // create perturbation vector + for (int di = 0; di < dim + 3; ++di) + { + v(di) = vec_pert[di]; + } + // perturbed vectors + mfem::Vector q_plus(q), q_minus(q); + mfem::Vector flux_plus(q), flux_minus(q); + adept::Stack diff_stack; + mfem::H1_FECollection fe_coll(1); //dummy + double Re = 1000.0; double Pr = 0.75; + mach::SAViscousSlipWallBC saslipinteg(diff_stack, &fe_coll, Re, Pr, sacs, 1.5); + // +ve perturbation + q_plus.Add(delta, v); + // -ve perturbation + q_minus.Add(-delta, v); + DYNAMIC_SECTION("SA Slip-Wall BC jacobian is correct w.r.t state") + { + + // get perturbed states flux vector + saslipinteg.calcFlux(x, nrm, 1.0, q_plus, Dw, flux_plus); + saslipinteg.calcFlux(x, nrm, 1.0, q_minus, Dw, flux_minus); + // compute the jacobian + saslipinteg.calcFluxJacState(x, nrm, 1.0, q, Dw, jac); + jac.Mult(v, jac_v); + // finite difference jacobian + mfem::Vector jac_v_fd(flux_plus); + jac_v_fd -= flux_minus; + jac_v_fd /= 2.0 * delta; + // compare each component of the matrix-vector products + for (int i = 0; i < dim + 3; ++i) + { + // std::cout << "FD " << i << " Deriv: " << jac_v_fd[i] << std::endl; + // std::cout << "AN " << i << " Deriv: " << jac_v[i] << std::endl; + REQUIRE(jac_v[i] == Approx(jac_v_fd[i]).margin(1e-10)); + } + } +} + +TEMPLATE_TEST_CASE_SIG("SAViscousSlipWallBC Dw Jacobian", "[SASlipWallBC]", + ((int dim), dim), 1, 2, 3) +{ + using namespace euler_data; + // copy the data into mfem vectors for convenience +// std::cout << "Slip-Wall Dw Jac" << std::endl; + mfem::Vector nrm(dim); mfem::Vector x(dim); + mfem::Vector q(dim + 3); + mfem::Vector flux(dim + 3); + mfem::Vector v(dim + 3); + mfem::Vector jac_v(dim + 3); + mfem::DenseMatrix Dw(delw_data, dim+3, dim); + mfem::DenseMatrix jac(dim + 3, dim + 3); + mfem::Vector sacs(13); + // create SA parameter vector + for (int di = 0; di < 13; ++di) + { + sacs(di) = sa_params[di]; + } + double delta = 1e-5; + + x = 0.0; + q(0) = rho; + q(dim + 1) = rhoe; + q(dim + 2) = 4; + for (int di = 0; di < dim; ++di) + { + q(di + 1) = rhou[di]; + Dw(dim+2, di) = 0.7+0.1*di; + } + // create direction vector + for (int di = 0; di < dim; ++di) + { + nrm(di) = dir[di]; + } + // create perturbation vector + for (int di = 0; di < dim + 3; ++di) + { + v(di) = vec_pert[di]; + } + // perturbed vectors + adept::Stack diff_stack; + mfem::H1_FECollection fe_coll(1); //dummy + double Re = 1000.0; double Pr = 0.75; + mach::SAViscousSlipWallBC saslipinteg(diff_stack, &fe_coll, Re, Pr, sacs, 1.5); + DYNAMIC_SECTION("SA Slip-Wall BC Dw jacobian is correct w.r.t state") + { + std::vector mat_vec_jac(dim); + for (int d = 0; d < dim; ++d) + { + mat_vec_jac[d].SetSize(dim+3); + } + // compute the jacobian + saslipinteg.calcFluxJacDw(x, nrm, 1.0, q, Dw, mat_vec_jac); + + for (int d = 0; d < dim; ++d) + { + // perturb one column of delw everytime + mfem::DenseMatrix Dw_plus(Dw), Dw_minus(Dw); + mfem::Vector flux_plus(q.Size()), flux_minus(q.Size()); + flux_plus = 0.0; flux_minus = 0.0; + for (int s = 0; s < dim+3; ++s) + { + Dw_plus.GetColumn(d)[s] += v(s) * delta; + Dw_minus.GetColumn(d)[s] -= v(s) * delta; + } + // get perturbed states flux vector + saslipinteg.calcFlux(x, nrm, 1.0, q, Dw_plus, flux_plus); + saslipinteg.calcFlux(x, nrm, 1.0, q, Dw_minus, flux_minus); + + mat_vec_jac[d].Mult(v, jac_v); + // finite difference jacobian + mfem::Vector jac_v_fd(flux_plus); + jac_v_fd -= flux_minus; + jac_v_fd /= 2.0 * delta; + // compare each component of the matrix-vector products + for (int i = 0; i < dim + 3; ++i) + { + // std::cout << "FD " << i << " Deriv: " << jac_v_fd[i] << std::endl; + // std::cout << "AN " << i << " Deriv: " << jac_v[i] << std::endl; + REQUIRE(jac_v[i] == Approx(jac_v_fd[i]).margin(1e-10)); + } + } + } +} + + +TEMPLATE_TEST_CASE_SIG("SAViscous Jacobian", "[SAViscous]", + ((int dim), dim), 1, 2, 3) +{ + using namespace euler_data; + // copy the data into mfem vectors for convenience + mfem::Vector q(dim + 3); + mfem::Vector conv(dim + 3); + mfem::Vector conv_plus(dim + 3); + mfem::Vector conv_minus(dim + 3); + mfem::Vector scale(dim + 3); + mfem::Vector scale_plus(dim + 3); + mfem::Vector scale_minus(dim + 3); + mfem::Vector scale_plus_2(dim + 3); + mfem::Vector scale_minus_2(dim + 3); + mfem::DenseMatrix Dw(delw_data, dim+3, dim); + mfem::Vector v(dim + 3); + mfem::DenseMatrix vm(dim, dim + 3); + mfem::Vector jac_v(dim + 3); + mfem::DenseMatrix adjJ(dim, dim); + mfem::DenseMatrix jac_conv(dim + 3, dim + 3); + mfem::DenseMatrix jac_scale(dim + 3, dim + 3); + double delta = 1e-5; + double Re = 1000000; double Pr = 1; + mfem::Vector sacs(13); + q(0) = rho; + q(dim + 1) = rhoe; + q(dim + 2) = 4; + for (int di = 0; di < dim; ++di) + { + q(di + 1) = rhou[di]; + Dw(dim+2, di) = 0.7; + } + // create SA parameter vector + for (int di = 0; di < 13; ++di) + { + sacs(di) = sa_params[di]; + } + // create perturbation vector + for (int di = 0; di < dim + 3; ++di) + { + v(di) = vec_pert[di]; + } + // create perturbation matrix + for (int di = 0; di < dim + 3; ++di) + { + for (int di2 = 0; di2 < dim; ++di2) + { + vm(di2, di) = (1.1 + 0.1*di2)*vec_pert[di]; + } + } + // perturbed vectors + mfem::Vector q_plus(q), q_minus(q); + adept::Stack diff_stack; + mach::SAViscousIntegrator saviscousinteg(diff_stack, Re, Pr, sacs); + for (int di = 0; di < dim; ++di) + { + DYNAMIC_SECTION("Jacobians w.r.t state is correct, dir"< saviscousinteg(diff_stack, Re, Pr, sacs, 1.0); + mach::ESViscousIntegrator esviscousinteg(diff_stack, Re, Pr, 1.0); + + + for (int di = 0; di < dim; ++di) + { + DYNAMIC_SECTION("Jacobians w.r.t Dw is correct, dir"< mat_vec_jac(dim); + std::vector mat_vec_jac_ns(dim); + for (int d = 0; d < dim; ++d) + { + mat_vec_jac[d].SetSize(dim+3); + mat_vec_jac_ns[d].SetSize(dim+2); + } + // compute the jacobian + saviscousinteg.applyScalingJacDw(di, v, q, Dw, mat_vec_jac); + esviscousinteg.applyScalingJacDw(di, vr, qr, Dwr, mat_vec_jac_ns); + + stringstream nssolname; stringstream ranssolname; + nssolname << "visc_jac_rans"; ranssolname << "visc_jac_ns"; + std::ofstream matlabns(nssolname.str()); matlabns.precision(15); + std::ofstream matlabrans(ranssolname.str()); matlabrans.precision(15); + mat_vec_jac[0].PrintMatlab(matlabrans); + mat_vec_jac_ns[0].PrintMatlab(matlabns); + + for (int d = 0; d < dim; ++d) + { + // perturb one column of delw everytime + mfem::DenseMatrix Dw_plus(Dw), Dw_minus(Dw); + for (int s = 0; s < dim+3; ++s) + { + Dw_plus.GetColumn(d)[s] += v(s) * delta; + Dw_minus.GetColumn(d)[s] -= v(s) * delta; + } + // get perturbed states conv vector + saviscousinteg.applyScaling(di, v, q, Dw_plus, scale_plus_2); + saviscousinteg.applyScaling(di, v, q, Dw_minus, scale_minus_2); + + mat_vec_jac[d].Mult(v, jac_v); + // finite difference jacobian + mfem::Vector jac_v_fd(dim+3); + jac_v_fd = scale_plus_2; + jac_v_fd -= scale_minus_2; + jac_v_fd /= 2.0 * delta; + //std::cout << "viscous applyScaling jac Dw:" << std::endl; + for (int j = 0; j < dim + 3; ++j) + { + // std::cout << "FD " << j << " Deriv: " << jac_v_fd[j] << std::endl; + // std::cout << "AN " << j << " Deriv: " << jac_v(j) << std::endl; + REQUIRE(jac_v(j) == Approx(jac_v_fd[j]).margin(1e-10)); + } + } + } + } +} + + +TEMPLATE_TEST_CASE_SIG("SALPS Jacobian", "[SALPS]", + ((int dim), dim), 1, 2, 3) +{ + using namespace euler_data; + // copy the data into mfem vectors for convenience + mfem::Vector q(dim + 3); + mfem::Vector conv(dim + 3); + mfem::Vector conv_plus(dim + 3); + mfem::Vector conv_minus(dim + 3); + mfem::Vector scale(dim + 3); + mfem::Vector scale_plus(dim + 3); + mfem::Vector scale_minus(dim + 3); + mfem::Vector scale_plus_2(dim + 3); + mfem::Vector scale_minus_2(dim + 3); + mfem::Vector vec(vec_pert, dim + 3); + mfem::Vector v(dim + 3); + mfem::Vector jac_v(dim + 3); + mfem::DenseMatrix adjJ(dim, dim); + mfem::DenseMatrix jac_conv(dim + 3, dim + 3); + mfem::DenseMatrix jac_scale(dim + 3, dim + 3); + mfem::DenseMatrix jac_scale_2(dim + 3, dim + 3); + double delta = 1e-5; + adjJ = 0.7; + q(0) = rho; + q(dim + 1) = rhoe; + q(dim + 2) = 4; + for (int di = 0; di < dim; ++di) + { + q(di + 1) = rhou[di]; + } + // create perturbation vector + for (int di = 0; di < dim + 3; ++di) + { + v(di) = vec_pert[di]; + } + // perturbed vectors + mfem::Vector q_plus(q), q_minus(q); + mfem::Vector vec_plus(vec), vec_minus(vec); + adept::Stack diff_stack; + mach::SALPSIntegrator salpsinteg(diff_stack); + // +ve perturbation + q_plus.Add(delta, v); + vec_plus.Add(delta, v); + // -ve perturbation + q_minus.Add(-delta, v); + vec_minus.Add(-delta, v); + DYNAMIC_SECTION("LPS salpsinteg convertVars jacobian is correct w.r.t state ") + { + // get perturbed states conv vector + salpsinteg.convertVars(q_plus, conv_plus); + salpsinteg.convertVars(q_minus, conv_minus); + salpsinteg.applyScaling(adjJ, q_plus, vec, scale_plus); + salpsinteg.applyScaling(adjJ, q_minus, vec, scale_minus); + salpsinteg.applyScaling(adjJ, q, vec_plus, scale_plus_2); + salpsinteg.applyScaling(adjJ, q, vec_minus, scale_minus_2); + // compute the jacobian + salpsinteg.convertVarsJacState(q, jac_conv); + salpsinteg.applyScalingJacState(adjJ, q, vec, jac_scale); + salpsinteg.applyScalingJacV(adjJ, q, jac_scale_2); + jac_conv.Mult(v, jac_v); + // finite difference jacobian + mfem::Vector jac_v_fd(conv_plus); + jac_v_fd -= conv_minus; + jac_v_fd /= 2.0 * delta; + // compare each component of the matrix-vector products + //std::cout << "convertVars jac:" << std::endl; + for (int i = 0; i < dim + 3; ++i) + { + REQUIRE(jac_v[i] == Approx(jac_v_fd[i]).margin(1e-10)); + } + + jac_scale.Mult(v, jac_v); + // finite difference jacobian + jac_v_fd = scale_plus; + jac_v_fd -= scale_minus; + jac_v_fd /= 2.0 * delta; + //std::cout << "applyScaling jac:" << std::endl; + for (int i = 0; i < dim + 3; ++i) + { + REQUIRE(jac_v[i] == Approx(jac_v_fd[i]).margin(1e-10)); + } + + jac_scale_2.Mult(v, jac_v); + // finite difference jacobian + jac_v_fd = scale_plus_2; + jac_v_fd -= scale_minus_2; + jac_v_fd /= 2.0 * delta; + //std::cout << "applyScaling jac vec:" << std::endl; + for (int i = 0; i < dim + 3; ++i) + { + REQUIRE(jac_v[i] == Approx(jac_v_fd[i]).margin(1e-10)); + } + } +} + +TEST_CASE("SAViscousIntegrator::AssembleElementGrad", "[SAViscousIntegrator]") +{ + using namespace mfem; + using namespace euler_data; + + const int dim = 2; // templating is hard here because mesh constructors + int num_state = dim + 3; + adept::Stack diff_stack; + double delta = 1e-5; + + mfem::Vector sacs(13); + // create SA parameter vector + for (int di = 0; di < 13; ++di) + { + sacs(di) = sa_params[di]; + } + + // generate a 2 element mesh + int num_edge = 2; + std::unique_ptr mesh(new Mesh(num_edge, num_edge, Element::TRIANGLE, + true /* gen. edges */, 1.0, 1.0, true)); + for (int p = 1; p <= 2; ++p) + { + DYNAMIC_SECTION("...for degree p = " << p) + { + std::unique_ptr fec( + new SBPCollection(p, dim)); + std::unique_ptr fes(new FiniteElementSpace( + mesh.get(), fec.get(), num_state, Ordering::byVDIM)); + + GridFunction dist; + NonlinearForm res(fes.get()); + res.AddDomainIntegrator(new mach::SAViscousIntegrator(diff_stack, 1000.0, 1.0, sacs)); + + // initialize state; here we randomly perturb a constant state + GridFunction q(fes.get()); + VectorFunctionCoefficient pert(num_state, randBaselinePertSA<2>); + q.ProjectCoefficient(pert); + + // initialize the vector that the Jacobian multiplies + GridFunction v(fes.get()); + VectorFunctionCoefficient v_rand(num_state, randVectorState); + v.ProjectCoefficient(v_rand); + + // evaluate the Ja;cobian and compute its product with v + Operator &Jac = res.GetGradient(q); + GridFunction jac_v(fes.get()); + Jac.Mult(v, jac_v); + + // now compute the finite-difference approximation... + GridFunction q_pert(q), r(fes.get()), jac_v_fd(fes.get()); + q_pert.Add(-delta, v); + res.Mult(q_pert, r); + q_pert.Add(2 * delta, v); + res.Mult(q_pert, jac_v_fd); + jac_v_fd -= r; + jac_v_fd /= (2 * delta); + + for (int i = 0; i < jac_v.Size(); ++i) + { + // std::cout << "FD " << i << " Deriv: " << jac_v_fd[i] << std::endl; + // std::cout << "AN " << i << " Deriv: " << jac_v(i) << std::endl; + REQUIRE(jac_v(i) == Approx(jac_v_fd(i)).margin(1e-10)); + } + } + } +} + + +TEST_CASE("SANoSlipAdiabaticWallBC::AssembleElementGrad", "[SANoSlipAdiabaticWallBC]") +{ + using namespace mfem; + using namespace euler_data; + + const int dim = 2; // templating is hard here because mesh constructors + int num_state = dim + 3; + mfem::Vector qfar(dim + 3); + adept::Stack diff_stack; + double delta = 1e-5; + + qfar(0) = rho; + qfar(dim + 1) = rhoe; + qfar(dim + 2) = 4; + for (int di = 0; di < dim; ++di) + { + qfar(di + 1) = rhou[di]; + } + + mfem::Vector sacs(13); + // create SA parameter vector + for (int di = 0; di < 13; ++di) + { + sacs(di) = sa_params[di]; + } + + // generate a 2 element mesh + int num_edge = 2; + std::unique_ptr mesh(new Mesh(num_edge, num_edge, Element::TRIANGLE, + true /* gen. edges */, 1.0, 1.0, true)); + for (int p = 1; p <= 2; ++p) + { + DYNAMIC_SECTION("...for degree p = " << p) + { + std::unique_ptr fec( + new SBPCollection(p, dim)); + std::unique_ptr fes(new FiniteElementSpace( + mesh.get(), fec.get(), num_state, Ordering::byVDIM)); + + GridFunction dist; + NonlinearForm res(fes.get()); + res.AddBdrFaceIntegrator(new mach::SANoSlipAdiabaticWallBC(diff_stack, fec.get(), 1000.0, 1.0, sacs, qfar, 1.0)); + + // initialize state; here we randomly perturb a constant state + GridFunction q(fes.get()); + VectorFunctionCoefficient pert(num_state, randBaselinePertSA<2>); + q.ProjectCoefficient(pert); + + // initialize the vector that the Jacobian multiplies + GridFunction v(fes.get()); + VectorFunctionCoefficient v_rand(num_state, randVectorState); + v.ProjectCoefficient(v_rand); + + // evaluate the Ja;cobian and compute its product with v + Operator &Jac = res.GetGradient(q); + GridFunction jac_v(fes.get()); + Jac.Mult(v, jac_v); + + // now compute the finite-difference approximation... + GridFunction q_pert(q), r(fes.get()), jac_v_fd(fes.get()); + q_pert.Add(-delta, v); + res.Mult(q_pert, r); + q_pert.Add(2 * delta, v); + res.Mult(q_pert, jac_v_fd); + jac_v_fd -= r; + jac_v_fd /= (2 * delta); + + for (int i = 0; i < jac_v.Size(); ++i) + { + // std::cout << "FD " << i << " Deriv: " << jac_v_fd[i] << std::endl; + // std::cout << "AN " << i << " Deriv: " << jac_v(i) << std::endl; + REQUIRE(jac_v(i) == Approx(jac_v_fd(i)).margin(1e-10)); + } + } + } +} + +TEST_CASE("SAViscousSlipWallBC::AssembleElementGrad", "[SAViscousSlipWallBC]") +{ + using namespace mfem; + using namespace euler_data; + + const int dim = 2; // templating is hard here because mesh constructors + int num_state = dim + 3; + mfem::Vector qfar(dim + 3); + adept::Stack diff_stack; + double delta = 1e-5; + + qfar(0) = rho; + qfar(dim + 1) = rhoe; + qfar(dim + 2) = 4; + for (int di = 0; di < dim; ++di) + { + qfar(di + 1) = rhou[di]; + } + + mfem::Vector sacs(13); + // create SA parameter vector + for (int di = 0; di < 13; ++di) + { + sacs(di) = sa_params[di]; + } + + // generate a 2 element mesh + int num_edge = 2; + std::unique_ptr mesh(new Mesh(num_edge, num_edge, Element::TRIANGLE, + true /* gen. edges */, 1.0, 1.0, true)); + for (int p = 1; p <= 2; ++p) + { + DYNAMIC_SECTION("...for degree p = " << p) + { + std::unique_ptr fec( + new SBPCollection(p, dim)); + std::unique_ptr fes(new FiniteElementSpace( + mesh.get(), fec.get(), num_state, Ordering::byVDIM)); + + GridFunction dist; + NonlinearForm res(fes.get()); + res.AddBdrFaceIntegrator(new mach::SAViscousSlipWallBC(diff_stack, fec.get(), 1000.0, 0.75, sacs, 1.5)); + + // initialize state; here we randomly perturb a constant state + GridFunction q(fes.get()); + VectorFunctionCoefficient pert(num_state, randBaselinePertSA<2>); + q.ProjectCoefficient(pert); + + // initialize the vector that the Jacobian multiplies + GridFunction v(fes.get()); + VectorFunctionCoefficient v_rand(num_state, randVectorState); + v.ProjectCoefficient(v_rand); + + // evaluate the Ja;cobian and compute its product with v + Operator &Jac = res.GetGradient(q); + GridFunction jac_v(fes.get()); + Jac.Mult(v, jac_v); + + // now compute the finite-difference approximation... + GridFunction q_pert(q), r(fes.get()), jac_v_fd(fes.get()); + q_pert.Add(-delta, v); + res.Mult(q_pert, r); + q_pert.Add(2 * delta, v); + res.Mult(q_pert, jac_v_fd); + jac_v_fd -= r; + jac_v_fd /= (2 * delta); + + for (int i = 0; i < jac_v.Size(); ++i) + { + // std::cout << "FD " << i << " Deriv: " << jac_v_fd[i] << std::endl; + // std::cout << "AN " << i << " Deriv: " << jac_v(i) << std::endl; + REQUIRE(jac_v(i) == Approx(jac_v_fd(i)).margin(1e-10)); + } + } + } +} + + + +TEST_CASE("SAFarFieldBC::AssembleElementGrad", "[SAFarFieldBC]") +{ + using namespace mfem; + using namespace euler_data; + + const int dim = 2; // templating is hard here because mesh constructors + int num_state = dim + 3; + mfem::Vector qfar(dim + 3); + adept::Stack diff_stack; + double delta = 1e-5; + + qfar(0) = rho; + qfar(dim + 1) = rhoe; + qfar(dim + 2) = 4; + for (int di = 0; di < dim; ++di) + { + qfar(di + 1) = rhou[di]; + } + + mfem::Vector sacs(13); + // create SA parameter vector + for (int di = 0; di < 13; ++di) + { + sacs(di) = sa_params[di]; + } + + // generate a 2 element mesh + int num_edge = 2; + std::unique_ptr mesh(new Mesh(num_edge, num_edge, Element::TRIANGLE, + true /* gen. edges */, 1.0, 1.0, true)); + for (int p = 1; p <= 2; ++p) + { + DYNAMIC_SECTION("...for degree p = " << p) + { + std::unique_ptr fec( + new SBPCollection(p, dim)); + std::unique_ptr fes(new FiniteElementSpace( + mesh.get(), fec.get(), num_state, Ordering::byVDIM)); + + GridFunction dist; + NonlinearForm res(fes.get()); + res.AddBdrFaceIntegrator(new mach::SAFarFieldBC(diff_stack, fec.get(), 1000.0, 1.0, qfar, 1.0)); + + // initialize state; here we randomly perturb a constant state + GridFunction q(fes.get()); + VectorFunctionCoefficient pert(num_state, randBaselinePertSA<2>); + q.ProjectCoefficient(pert); + + // initialize the vector that the Jacobian multiplies + GridFunction v(fes.get()); + VectorFunctionCoefficient v_rand(num_state, randVectorState); + v.ProjectCoefficient(v_rand); + + // evaluate the Ja;cobian and compute its product with v + Operator &Jac = res.GetGradient(q); + GridFunction jac_v(fes.get()); + Jac.Mult(v, jac_v); + + // now compute the finite-difference approximation... + GridFunction q_pert(q), r(fes.get()), jac_v_fd(fes.get()); + q_pert.Add(-delta, v); + res.Mult(q_pert, r); + q_pert.Add(2 * delta, v); + res.Mult(q_pert, jac_v_fd); + jac_v_fd -= r; + jac_v_fd /= (2 * delta); + + for (int i = 0; i < jac_v.Size(); ++i) + { + // std::cout << "FD " << i << " Deriv: " << jac_v_fd[i] << std::endl; + // std::cout << "AN " << i << " Deriv: " << jac_v(i) << std::endl; + REQUIRE(jac_v(i) == Approx(jac_v_fd(i)).margin(1e-10)); + } + } + } +} + + +///NOTE: Not currently active +TEST_CASE("SASourceIntegrator::AssembleElementGrad", "[SASourceIntegrator]") +{ + using namespace mfem; + using namespace euler_data; + + const int dim = 2; // templating is hard here because mesh constructors + int num_state = dim + 3; + adept::Stack diff_stack; + double delta = 1e-5; + + mfem::Vector sacs(13); + // create SA parameter vector + for (int di = 0; di < 13; ++di) + { + sacs(di) = sa_params[di]; + } + + // generate a 8 element mesh + int num_edge = 3; + std::unique_ptr mesh(new Mesh(num_edge, num_edge, Element::TRIANGLE, + true /* gen. edges */, 1.0, 1.0, true)); + + + + for (int p = 1; p <= 1; ++p) + { + DYNAMIC_SECTION("...for degree p = " << p) + { + std::unique_ptr fec( + new SBPCollection(p, dim)); + std::unique_ptr fes(new FiniteElementSpace( + mesh.get(), fec.get(), num_state, Ordering::byVDIM)); + std::unique_ptr fesx(new FiniteElementSpace( + mesh.get(), fec.get(), 2, Ordering::byVDIM)); + + mesh->EnsureNodes(); + + GridFunction dist(fes.get()); + GridFunction x_nodes(fesx.get()); + auto walldist = [](const Vector &x) + { + return x(1); + }; + auto coord = [](const Vector &x, Vector &u) + { + u = x; + }; + FunctionCoefficient wall_coeff(walldist); + VectorFunctionCoefficient x_coord(2, coord); + dist.ProjectCoefficient(wall_coeff); + x_nodes.ProjectCoefficient(x_coord); + //dist = 1.0000; + NonlinearForm res(fes.get()); + double Re = 5000000.0; + res.AddDomainIntegrator(new mach::SASourceIntegrator(diff_stack, dist, Re, sacs, 1.0, -1.0, 1.0, 1.0)); + + // initialize state; here we randomly perturb a constant state + GridFunction q(fes.get()); + VectorFunctionCoefficient pert(num_state, randBaselinePertSA<2>); + q.ProjectCoefficient(pert); + + // initialize the vector that the Jacobian multiplies + GridFunction v(fes.get()); + VectorFunctionCoefficient v_rand(num_state, randVectorState); + v.ProjectCoefficient(v_rand); + + // evaluate the Jacobian and compute its product with v + Operator &Jac = res.GetGradient(q); + + // stringstream nssolname; + // nssolname << "jac_test"; + // std::ofstream matlab(nssolname.str()); matlab.precision(15); + // Jac.PrintMatlab(matlab); + + GridFunction jac_v(fes.get()); + Jac.Mult(v, jac_v); + + // now compute the finite-difference approximation... + GridFunction q_pert(q), r(fes.get()), jac_v_fd(fes.get()); + q_pert.Add(-delta, v); + res.Mult(q_pert, r); + q_pert.Add(2 * delta, v); + res.Mult(q_pert, jac_v_fd); + jac_v_fd -= r; + jac_v_fd /= (2 * delta); + + GridFunction error(fes.get()); + //jac_v.Size()); + error = jac_v; + error -= jac_v_fd; + + // res.Mult(q, r); + // ofstream sol_ofs("why.vtk"); + // sol_ofs.precision(14); + // mesh->PrintVTK(sol_ofs, 1); + // r.SaveVTK(sol_ofs, "jac_error", 1); + + std::cout.precision(17); + std::cout << "Error Norm: " << error.Norml2() << std::endl; + for (int i = 0; i < jac_v.Size(); ++i) + { + int n = i/5; + // std::cout << "Node Coord: "<< x_nodes(0+2*n) <<", "<< x_nodes(1+2*n) < mesh(new Mesh(num_edge, num_edge, Element::TRIANGLE, + true /* gen. edges */, 1.0, 1.0, true)); + for (int p = 1; p <= 4; ++p) + { + DYNAMIC_SECTION("...for degree p = " << p) + { + std::unique_ptr fec( + new SBPCollection(p, dim)); + std::unique_ptr fes(new FiniteElementSpace( + mesh.get(), fec.get(), num_state, Ordering::byVDIM)); + + NonlinearForm res(fes.get()); + res.AddDomainIntegrator(new mach::SALPSIntegrator(diff_stack)); + + // initialize state; here we randomly perturb a constant state + GridFunction q(fes.get()); + VectorFunctionCoefficient pert(num_state, randBaselinePertSA<2>); + q.ProjectCoefficient(pert); + + // initialize the vector that the Jacobian multiplies + GridFunction v(fes.get()); + VectorFunctionCoefficient v_rand(num_state, randVectorState); + v.ProjectCoefficient(v_rand); + + // evaluate the Jacobian and compute its product with v + Operator &Jac = res.GetGradient(q); + GridFunction jac_v(fes.get()); + Jac.Mult(v, jac_v); + + // now compute the finite-difference approximation... + GridFunction q_pert(q), r(fes.get()), jac_v_fd(fes.get()); + q_pert.Add(-delta, v); + res.Mult(q_pert, r); + q_pert.Add(2 * delta, v); + res.Mult(q_pert, jac_v_fd); + jac_v_fd -= r; + jac_v_fd /= (2 * delta); + + for (int i = 0; i < jac_v.Size(); ++i) + { + REQUIRE(jac_v(i) == Approx(jac_v_fd(i)).margin(1e-10)); + } + } + } +} + +#if 0 +TEMPLATE_TEST_CASE_SIG("SAInviscid Gradient", + "[SAInviscidGrad]", + ((bool entvar), entvar), false) +{ + using namespace mfem; + using namespace euler_data; + + const int dim = 2; // templating is hard here because mesh constructors + double delta = 1e-5; + + // generate a 8 element mesh + int num_edge = 2; + std::unique_ptr mesh(new Mesh(num_edge, num_edge, Element::TRIANGLE, + true /* gen. edges */, 1.0, 1.0, true)); + mesh->EnsureNodes(); + for (int p = 1; p <= 4; ++p) + { + DYNAMIC_SECTION("...for degree p = " << p) + { + // get the finite-element space for the state and adjoint + std::unique_ptr fec( + new SBP_FECollection(p, dim)); + std::unique_ptr fes(new FiniteElementSpace( + mesh.get(), fec.get())); + + // we use res for finite-difference approximation + NonLinearForm res(fes.get()); + res.AddDomainIntegrator( + new SAInviscidIntegrator( //Inviscid term + this->diff_stack, alpha)); + + // initialize state ; here we randomly perturb a constant state + GridFunction state(fes.get()) + VectorFunctionCoefficient v_rand(dim+3, randVectorState); + state.ProjectCoefficient(pert); + + // initialize the vector that we use to perturb the mesh nodes + GridFunction v(fes.get()); + VectorFunctionCoefficient v_rand(dim, randVectorState); + v.ProjectCoefficient(v_rand); + + // evaluate df/dx and contract with v + GridFunction dfdx(*x_nodes); + dfdx_form.Mult(*x_nodes, dfdx); + double dfdx_v = dfdx * v; + + // now compute the finite-difference approximation... + GridFunction x_pert(*x_nodes); + GridFunction r(fes.get()); + x_pert.Add(delta, v); + mesh->SetNodes(x_pert); + res.Assemble(); + double dfdx_v_fd = adjoint * res; + x_pert.Add(-2 * delta, v); + mesh->SetNodes(x_pert); + res.Assemble(); + dfdx_v_fd -= adjoint * res; + dfdx_v_fd /= (2 * delta); + mesh->SetNodes(*x_nodes); // remember to reset the mesh nodes + + REQUIRE(dfdx_v == Approx(dfdx_v_fd).margin(1e-10)); + } + } +} +#endif + +void uinit(const mfem::Vector &x, mfem::Vector& q) +{ + // q.SetSize(4); + // Vector u(4); + q.SetSize(5); + mfem::Vector u(5); + + u = 0.0; + u(0) = 1.0; + u(1) = 1.5;//u(0)*1.5*cos(0.3*M_PI/180.0); + u(2) = 0;//u(0)*1.5*sin(0.3*M_PI/180.0); + u(3) = 1.0/(1.4) + 0.5*1.5*1.5; + u(4) = u(0)*3; + + q = u; +} \ No newline at end of file