Skip to content

Commit

Permalink
extended cases-treatment to specific vel bc and temp ic/bc
Browse files Browse the repository at this point in the history
  • Loading branch information
Sigfried Haering committed Jun 1, 2024
1 parent a9201d6 commit f845e1b
Show file tree
Hide file tree
Showing 5 changed files with 113 additions and 134 deletions.
32 changes: 8 additions & 24 deletions src/calorically_perfect.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -38,6 +38,7 @@
#include <iomanip>
#include <mfem/general/forall.hpp>

#include "cases.hpp"
#include "loMach.hpp"
#include "loMach_options.hpp"
#include "logger.hpp"
Expand All @@ -47,8 +48,7 @@
using namespace mfem;
using namespace mfem::common;

/// forward declarations
double temp_rt3d(const Vector &x, double t);
/// forward declarations (if any)

MFEM_HOST_DEVICE double Sutherland(const double T, const double mu_star, const double T_star, const double S_star) {
const double T_rat = T / T_star;
Expand Down Expand Up @@ -265,12 +265,18 @@ void CaloricallyPerfectThermoChem::initializeSelf() {

// set IC if we have one at this point
if (!ic_string_.empty()) {
/*
if (ic_string_ == "rt3D") {
if (rank0_) std::cout << "Setting rt3D IC..." << std::endl;
FunctionCoefficient t_excoeff(temp_rt3d);
t_excoeff.SetTime(0.0);
Tn_gf_.ProjectCoefficient(t_excoeff);
}
*/
sfptr user_func = temp_ic(ic_string_);
FunctionCoefficient t_excoeff(user_func);
t_excoeff.SetTime(0.0);
Tn_gf_.ProjectCoefficient(t_excoeff);
} else {
ConstantCoefficient t_ic_coef;
t_ic_coef.constant = T_ic_;
Expand Down Expand Up @@ -1455,25 +1461,3 @@ double temp_inlet(const Vector &coords, double t) {
return temp;
}
#endif

double temp_rt3d(const Vector &x, double t) {
double CC = 0.05;
double twoPi = 6.28318530718;
double yWidth = 0.1;
double yInt, dy, wt;
double temp, dT;
double Tlo = 100.0;
double Thi = 1500.0;

yInt = std::cos(twoPi * x[0]) + std::cos(twoPi * x[2]);
yInt *= CC;
yInt += 4.0;

dy = x[1] - yInt;
dT = Thi - Tlo;

wt = 0.5 * (tanh(-dy / yWidth) + 1.0);
temp = Tlo + wt * dT;

return temp;
}
82 changes: 65 additions & 17 deletions src/cases.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -36,23 +36,23 @@

#include "cases.hpp"

#include <grvy.h>
#include <sys/stat.h>
#include <unistd.h>

#include <array>
#include <cstdio>
#include <cstdlib>
#include <iostream>
#include <memory>
#include <stdexcept>
#include <string>
#include <grvy.h>

#include "tps_mfem_wrap.hpp"
#include "utils.hpp"

using namespace mfem;


/// Used to set the velocity IC (and to check error)
void vel_exact_tgv2d(const Vector &x, double t, Vector &u) {
const double nu = 1.0;
Expand Down Expand Up @@ -116,20 +116,68 @@ void vel_channel(const Vector &x, double t, Vector &u) {
}
}

/// Add ic cases to selection here
vfptr vel_ic(std::string ic_string_) {
if (ic_string_ == "tgv2d") {
return vel_exact_tgv2d;
} else if (ic_string_ == "tgv2d_uniform") {
return vel_tgv2d_uniform;
} else if (ic_string_ == "channel") {
return vel_channel;
} else {
grvy_printf(GRVY_ERROR, "Attempting to use unknown vel_ic");
exit(ERROR);
return vel_exact_tgv2d;
}
}

/// Used to for pipe flow test case
void vel_exact_pipe(const Vector &x, double t, Vector &u) {
u(0) = 0.0;
u(1) = 2.0 * (1 - x[0] * x[0]);
}

/// Add bc cases to selection here
vfptr vel_bc(std::string type) {
if (type == "fully-developed-pipe") {
return vel_exact_pipe;
} else {
grvy_printf(GRVY_ERROR, "Attempting to use unknown vel_bc");
exit(ERROR);
return vel_exact_pipe;
}
}

/// Rayleigh-Taylor ic
double temp_rt3d(const Vector &x, double t) {
double CC = 0.05;
double twoPi = 6.28318530718;
double yWidth = 0.1;
double yInt, dy, wt;
double temp, dT;
double Tlo = 100.0;
double Thi = 1500.0;

yInt = std::cos(twoPi * x[0]) + std::cos(twoPi * x[2]);
yInt *= CC;
yInt += 4.0;

dy = x[1] - yInt;
dT = Thi - Tlo;

// typedef int (*fptr)();

/// Add cases to selection here
fptr vel_ic(std::string ic_string_) {
if (ic_string_ == "tgv2d") {
return vel_exact_tgv2d;
} else if (ic_string_ == "tgv2d_uniform") {
return vel_tgv2d_uniform;
} else if (ic_string_ == "channel") {
return vel_channel;
} else {
grvy_printf(GRVY_ERROR, "Attempting to use unknown vel_ic");
exit(ERROR);
return vel_exact_tgv2d;
}
wt = 0.5 * (tanh(-dy / yWidth) + 1.0);
temp = Tlo + wt * dT;

return temp;
}

/// Add temp ic cases to selection here
sfptr temp_ic(std::string ic_string_) {
if (ic_string_ == "rt3D") {
return temp_rt3d;
} else {
grvy_printf(GRVY_ERROR, "Attempting to use unknown temp_ic");
exit(ERROR);
return temp_rt3d;
}
}
18 changes: 14 additions & 4 deletions src/cases.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -33,20 +33,30 @@
#define CASES_HPP_

/** @file
* @brief A place for case-specific bc and ic.
* @brief A place for case-specific bc and ic. All case names specified in the input
* file must be added to if-list in cases.cpp along with the actual ic/bc.
*/

#include <assert.h>
#include <hdf5.h>

#include <string>

#include "tps_mfem_wrap.hpp"

// typedef std::function<void(const Vector &, double, Vector &)> (*fptr)();
typedef std::function<void(const mfem::Vector &, double, mfem::Vector &)> fptr;
fptr vel_ic(std::string ic_string_);
typedef std::function<void(const mfem::Vector &, double, mfem::Vector &)> vfptr;
vfptr vel_ic(std::string ic_string_);
vfptr vel_bc(std::string type);

typedef std::function<double(const mfem::Vector &, double)> sfptr;
sfptr temp_ic(std::string ic_string_);
sfptr temp_bc(std::string type);

void vel_exact_tgv2d(const mfem::Vector &x, double t, mfem::Vector &u);
void vel_tgv2d_uniform(const mfem::Vector &x, double t, mfem::Vector &u);
void vel_channel(const mfem::Vector &x, double t, mfem::Vector &u);
void vel_exact_pipe(const mfem::Vector &x, double t, mfem::Vector &u);

double temp_rt3d(const mfem::Vector &x, double t);

#endif // CASES_HPP_
112 changes: 24 additions & 88 deletions src/tomboulides.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -35,22 +35,17 @@
#include <mfem/general/forall.hpp>

#include "algebraicSubgridModels.hpp"
#include "cases.hpp"
#include "externalData_base.hpp"
#include "io.hpp"
#include "loMach.hpp"
#include "thermo_chem_base.hpp"
#include "tps.hpp"
#include "utils.hpp"
#include "cases.hpp"

using namespace mfem;

/// forward declarations
// void vel_exact_tgv2d(const Vector &x, double t, Vector &u);
// void vel_tgv2d_uniform(const Vector &x, double t, Vector &u);
void vel_exact_pipe(const Vector &x, double t, Vector &u);
// void vel_channel(const Vector &x, double t, Vector &u);

static double radius(const Vector &pos) { return pos[0]; }
FunctionCoefficient radius_coeff(radius);

Expand Down Expand Up @@ -388,14 +383,13 @@ void Tomboulides::initializeSelf() {

// set IC if we have one at this point
if (!ic_string_.empty()) {
// std::function *user_func;
fptr user_func = vel_ic(ic_string_);
vfptr user_func = vel_ic(ic_string_);
VectorFunctionCoefficient u_excoeff(nvel_, user_func);
u_excoeff.SetTime(0.0);
u_curr_gf_->ProjectCoefficient(u_excoeff);
}

/*
/*
if (ic_string_ == "tgv2d") {
if (rank0_) std::cout << "Setting tgv2d IC..." << std::endl;
VectorFunctionCoefficient u_excoeff(2, vel_exact_tgv2d);
Expand All @@ -413,7 +407,7 @@ void Tomboulides::initializeSelf() {
u_curr_gf_->ProjectCoefficient(u_excoeff);
}
}
*/
*/

// Boundary conditions
// number of BC regions defined
Expand Down Expand Up @@ -468,6 +462,20 @@ void Tomboulides::initializeSelf() {
// axisymmetric not testsed with interpolation BCs yet. For now, just stop.
assert(!axisym_);

} else {
Array<int> inlet_attr(pmesh_->bdr_attributes.Max());
inlet_attr = 0;
inlet_attr[patch - 1] = 1;

vfptr user_func = vel_bc(type);
addVelDirichletBC(user_func, inlet_attr);

if (axisym_) {
addSwirlDirichletBC(0.0, inlet_attr);
}
}

/*
} else if (type == "fully-developed-pipe") {
if (pmesh_->GetMyRank() == 0) {
std::cout << "Tomboulides: Setting uniform inlet velocity on patch = " << patch << std::endl;
Expand All @@ -477,7 +485,6 @@ void Tomboulides::initializeSelf() {
inlet_attr[patch - 1] = 1;
addVelDirichletBC(vel_exact_pipe, inlet_attr);

if (axisym_) {
addSwirlDirichletBC(0.0, inlet_attr);
}
Expand All @@ -489,6 +496,7 @@ void Tomboulides::initializeSelf() {
assert(false);
exit(1);
}
*/
}

// Wall Bcs
Expand Down Expand Up @@ -1549,12 +1557,11 @@ void Tomboulides::step() {
}
}


double Tomboulides::computeL2Error() const {
double err = -1.0;
if (ic_string_ == "tgv2d") {
std::cout << "Evaluating TGV2D error..." << std::endl;
fptr user_func = vel_ic(ic_string_);
vfptr user_func = vel_ic(ic_string_);
VectorFunctionCoefficient u_excoeff(nvel_, user_func);
u_excoeff.SetTime(coeff_.time);
err = u_curr_gf_->ComputeL2Error(u_excoeff);
Expand Down Expand Up @@ -1619,6 +1626,10 @@ void Tomboulides::addVelDirichletBC(void (*f)(const Vector &, double, Vector &),
addVelDirichletBC(new VectorFunctionCoefficient(dim_, f), attr);
}

void Tomboulides::addVelDirichletBC(std::function<void(const Vector &, double, Vector &)> f, Array<int> &attr) {
addVelDirichletBC(new VectorFunctionCoefficient(dim_, f), attr);
}

/// Add a Dirichlet boundary condition to the pressure field.
void Tomboulides::addPresDirichletBC(double p, Array<int> &attr) {
pres_dbcs_.emplace_back(attr, new ConstantCoefficient(p));
Expand Down Expand Up @@ -1685,78 +1696,3 @@ void Tomboulides::evaluateVelocityGradient() {
gradV_gf_->SetFromTrueDofs(gradV_);
gradW_gf_->SetFromTrueDofs(gradW_);
}

// Non-class functions that are only used in this file below here

/*
/// Used to set the velocity IC (and to check error)
void vel_exact_tgv2d(const Vector &x, double t, Vector &u) {
const double nu = 1.0;
const double F = std::exp(-2 * nu * t);
u(0) = F * std::sin(x[0]) * std::cos(x[1]);
u(1) = -F * std::cos(x[0]) * std::sin(x[1]);
}
*/

/// Used to for pipe flow test case
void vel_exact_pipe(const Vector &x, double t, Vector &u) {
u(0) = 0.0;
u(1) = 2.0 * (1 - x[0] * x[0]);
}

/*
/// Used to set the velocity IC with TG field and uniform
void vel_tgv2d_uniform(const Vector &x, double t, Vector &u) {
const double u0 = 1.0;
const double F = 0.1;
const double PI = 3.14159265359;
double twoPi = 2.0 * PI;
u(0) = u0;
u(1) = 0.0;
u(0) += +F * std::sin(twoPi * x[0]) * std::cos(twoPi * x[1]);
u(1) += -F * std::cos(twoPi * x[0]) * std::sin(twoPi * x[1]);
}
/// Used to set the channel IC
void vel_channel(const Vector &x, double t, Vector &u) {
double PI = 3.14159265359;
double Lx = 25.0;
double Ly = 2.0;
double Lz = 9.4;
double Umean = 1.0;
double uInt = 0.1;
int nModes = 4;
double uM;
double ax, by, cz;
double AA, BB, CC;
double wall;
// expects channel height (-1,1)
wall = (1.0 - std::pow(x(1), 8.0));
u(0) = Umean * wall;
u(1) = 0.0;
u(2) = 0.0;
for (int n = 1; n <= nModes; n++) {
ax = 4.0 * PI / Lx * (double)n;
by = 2.0 * PI / Ly * (double)n;
cz = 2.0 * PI / Lz * (double)n;
AA = 1.0;
BB = 1.0;
CC = -(AA * ax + BB * by) / cz;
uM = uInt / (double)n;
u(0) += uM * AA * cos(ax * (x(0) + (double)(n - 1) * Umean)) * sin(by * x(1)) *
sin(cz * (x(2) + 0.5 * (double)(n - 1) * Umean)) * wall;
u(1) += uM * BB * sin(ax * (x(0) + (double)(n - 1) * Umean)) * cos(by * x(1)) *
sin(cz * (x(2) + 0.5 * (double)(n - 1) * Umean)) * wall;
u(2) += uM * CC * sin(ax * (x(0) + (double)(n - 1) * Umean)) * sin(by * x(1)) *
cos(cz * (x(2) + 0.5 * (double)(n - 1) * Umean)) * wall;
}
}
*/
Loading

0 comments on commit f845e1b

Please sign in to comment.