diff --git a/.gitattributes b/.gitattributes index 6e9ba2359..283e1a89d 100644 --- a/.gitattributes +++ b/.gitattributes @@ -1,3 +1,5 @@ test/ref_solns/cyl3d.mflow.2iters.bulkVisc.dtconst.h5 filter=lfs diff=lfs merge=lfs -text test/ref_solns/cyl.r0.p2.iter4.h5 filter=lfs diff=lfs merge=lfs -text test/ref_solns/cyl.r1.p2.iter4.h5 filter=lfs diff=lfs merge=lfs -text +test/ref_solns/forcing-periodic-square.h5 filter=lfs diff=lfs merge=lfs -text +test/meshes/periodic-square.mesh filter=lfs diff=lfs merge=lfs -text diff --git a/src/dataStructures.hpp b/src/dataStructures.hpp index e08e56e7f..6461db6cd 100644 --- a/src/dataStructures.hpp +++ b/src/dataStructures.hpp @@ -42,6 +42,10 @@ using namespace mfem; // normal equal for both planes. The visc ratio will vary linearly // from 0 to viscRatio from the plane defined by pointInit and // the plane defined by point0. +// |-----------------------------------| +// |------------------------------<====| +// |-----------------------------------| +// p_init normal p0 struct linearlyVaryingVisc { Vector normal; Vector point0; @@ -49,6 +53,33 @@ struct linearlyVaryingVisc { double viscRatio; }; +enum SpongeZoneSolution { + USERDEF, // User defined target solution + MIXEDOUT, // Mixed-out values will be computed at the pInit plane and + // used as target solution in the sponge zone + NONE = -1 // Invalid value +}; + +struct SpongeZoneData { + Vector normal; // These planes are defined in the same manner as + Vector point0; // in the linearlyVaryingVisc struct case. + Vector pointInit; + + SpongeZoneSolution szType; + + double tol; // Tolerance for finding nodes at pInit plane + // These points are used to compute mixed-out values. + // Required if szType==MIXEDOUT + + Vector targetUp; // Target user-defined solution defined with primitive + // variables: rho, V, p. Required if szType==USERDEF + + double multFactor; // Factor multiplying the product of sigma*(U-U*) where + // U* is the target solution. Currently sigam=a/l where + // a is the mixed-out speed of sound and l is the lenght + // of the spnge zone. [OPTIONAL] +}; + struct volumeFaceIntegrationArrays { // nodes IDs and indirection array Array nodesIDs; diff --git a/src/em_options.hpp b/src/em_options.hpp index 9b3258d0a..1f1a8d4ab 100644 --- a/src/em_options.hpp +++ b/src/em_options.hpp @@ -59,54 +59,48 @@ class ElectromagneticOptions { double yinterp_max; /**< End value for uniformly spaced points */ ElectromagneticOptions() - : - mesh_file("hello.msh"), order(1), ref_levels(0), - max_iter(100), rtol(1e-6), atol(1e-10), - top_only(false), bot_only(false), - By_file("By.h5"), nBy(0), - yinterp_min(0.0), yinterp_max(1.0) - { } + : mesh_file("hello.msh"), + order(1), + ref_levels(0), + max_iter(100), + rtol(1e-6), + atol(1e-10), + top_only(false), + bot_only(false), + By_file("By.h5"), + nBy(0), + yinterp_min(0.0), + yinterp_max(1.0) {} void AddElectromagneticOptions(mfem::OptionsParser &args) { args.AddOption(&mesh_file, "-m", "--mesh", "Mesh file (for EM-only simulation)"); - args.AddOption(&order, "-o", "--order", - "Finite element order (polynomial degree) (for EM-only)."); - args.AddOption(&ref_levels, "-r", "--ref", - "Number of uniform refinements (for EM-only)."); - args.AddOption(&max_iter, "-i", "--maxiter", - "Maximum number of iterations (for EM-only)."); - args.AddOption(&rtol, "-t", "--rtol", - "Solver relative tolerance (for EM-only)."); - args.AddOption(&atol, "-a", "--atol", - "Solver absolute tolerance (for EM-only)."); - args.AddOption(&top_only, "-top", "--top-only", "-ntop", "--no-top-only", - "Run current through top branch only"); - args.AddOption(&bot_only, "-bot", "--bot-only", "-nbot", "--no-bot-only", - "Run current through bottom branch only"); - args.AddOption(&By_file, "-by", "--byfile", - "File for By interpolant output."); - args.AddOption(&nBy, "-ny", "--nyinterp", - "Number of interpolation points."); - args.AddOption(&yinterp_min, "-y0", "--yinterpMin", - "Minimum y interpolation value"); - args.AddOption(&yinterp_max, "-y1", "--yinterpMax", - "Maximum y interpolation value"); + args.AddOption(&order, "-o", "--order", "Finite element order (polynomial degree) (for EM-only)."); + args.AddOption(&ref_levels, "-r", "--ref", "Number of uniform refinements (for EM-only)."); + args.AddOption(&max_iter, "-i", "--maxiter", "Maximum number of iterations (for EM-only)."); + args.AddOption(&rtol, "-t", "--rtol", "Solver relative tolerance (for EM-only)."); + args.AddOption(&atol, "-a", "--atol", "Solver absolute tolerance (for EM-only)."); + args.AddOption(&top_only, "-top", "--top-only", "-ntop", "--no-top-only", "Run current through top branch only"); + args.AddOption(&bot_only, "-bot", "--bot-only", "-nbot", "--no-bot-only", "Run current through bottom branch only"); + args.AddOption(&By_file, "-by", "--byfile", "File for By interpolant output."); + args.AddOption(&nBy, "-ny", "--nyinterp", "Number of interpolation points."); + args.AddOption(&yinterp_min, "-y0", "--yinterpMin", "Minimum y interpolation value"); + args.AddOption(&yinterp_max, "-y1", "--yinterpMax", "Maximum y interpolation value"); } void print(std::ostream &out) { out << std::endl; out << "Electromagnetics options:" << std::endl; - out << " mesh_file = " << mesh_file << std::endl; - out << " order = " << order << std::endl; - out << " ref_levels = " << ref_levels << std::endl; - out << " ref_levels = " << ref_levels << std::endl; - out << " max_iter = " << max_iter << std::endl; - out << " rtol = " << rtol << std::endl; - out << " atol = " << atol << std::endl; - out << " top_only = " << top_only << std::endl; - out << " bot_only = " << top_only << std::endl; - out << " By_file = " << By_file << std::endl; - out << " nBy = " << nBy << std::endl; + out << " mesh_file = " << mesh_file << std::endl; + out << " order = " << order << std::endl; + out << " ref_levels = " << ref_levels << std::endl; + out << " ref_levels = " << ref_levels << std::endl; + out << " max_iter = " << max_iter << std::endl; + out << " rtol = " << rtol << std::endl; + out << " atol = " << atol << std::endl; + out << " top_only = " << top_only << std::endl; + out << " bot_only = " << top_only << std::endl; + out << " By_file = " << By_file << std::endl; + out << " nBy = " << nBy << std::endl; out << " yinterp_min = " << yinterp_min << std::endl; out << " yinterp_max = " << yinterp_max << std::endl; out << std::endl; diff --git a/src/equation_of_state.cpp b/src/equation_of_state.cpp index c3f2a13f2..a72735c37 100644 --- a/src/equation_of_state.cpp +++ b/src/equation_of_state.cpp @@ -109,6 +109,28 @@ double EquationOfState::ComputeMaxCharSpeed(const Vector& state, const int dim) return vel + sound; } +void EquationOfState::GetConservativesFromPrimitives(const Vector& primit, Vector& conserv, const int& dim, + const int& num_equations) { + conserv = primit; + + double v2 = 0.; + for (int d = 0; d < dim; d++) { + v2 += primit[1 + d]; + conserv[1 + d] *= primit[0]; + } + conserv[num_equations - 1] = primit[num_equations - 1] / (specific_heat_ratio - 1.) + 0.5 * primit[0] * v2; +} + +void EquationOfState::GetPrimitivesFromConservatives(const Vector& conserv, Vector& primit, const int& dim, + const int& num_equations) { + double p = ComputePressure(conserv, dim); + primit = conserv; + + for (int d = 0; d < dim; d++) primit[1 + d] /= conserv[0]; + + primit[num_equations - 1] = p; +} + // GPU FUNCTIONS /*#ifdef _GPU_ double EquationOfState::pressure( double *state, diff --git a/src/equation_of_state.hpp b/src/equation_of_state.hpp index 26c8b5dc3..07339e4c0 100644 --- a/src/equation_of_state.hpp +++ b/src/equation_of_state.hpp @@ -71,6 +71,9 @@ class EquationOfState { double ComputePressure(const Vector &state, int dim); + void GetPrimitivesFromConservatives(const Vector &conserv, Vector &primit, const int &dim, const int &num_equations); + void GetConservativesFromPrimitives(const Vector &primit, Vector &conserv, const int &dim, const int &num_equations); + // Compute the maximum characteristic speed. double ComputeMaxCharSpeed(const Vector &state, const int dim); diff --git a/src/forcing_terms.cpp b/src/forcing_terms.cpp index 3559315a8..bbe9cecdf 100644 --- a/src/forcing_terms.cpp +++ b/src/forcing_terms.cpp @@ -32,6 +32,7 @@ #include "forcing_terms.hpp" #include +#include ForcingTerms::ForcingTerms(const int &_dim, const int &_num_equation, const int &_order, const int &_intRuleType, IntegrationRules *_intRules, ParFiniteElementSpace *_vfes, ParGridFunction *_Up, @@ -213,6 +214,159 @@ void ConstantPressureGradient::updateTerms_gpu(const int numElems, const int off #endif +SpongeZone::SpongeZone(const int &_dim, const int &_num_equation, const int &_order, const int &_intRuleType, + Fluxes *_fluxClass, EquationOfState *_eqState, IntegrationRules *_intRules, + ParFiniteElementSpace *_vfes, ParGridFunction *_Up, ParGridFunction *_gradUp, + const volumeFaceIntegrationArrays &gpuArrays, RunConfiguration &_config) + : ForcingTerms(_dim, _num_equation, _order, _intRuleType, _intRules, _vfes, _Up, _gradUp, gpuArrays), + fluxes(_fluxClass), + eqState(_eqState), + szData(_config.GetSpongeZoneData()) { + targetU.SetSize(num_equation); + if (szData.szType == SpongeZoneSolution::USERDEF) { + Vector Up(num_equation); + Up[0] = szData.targetUp[0]; + for (int d = 0; d < dim; d++) Up[1 + d] = szData.targetUp[1 + d]; + Up[1 + dim] = szData.targetUp[4]; + eqState->GetConservativesFromPrimitives(Up, targetU, dim, num_equation); + } + + meanNormalFluxes.SetSize(num_equation + 1); + + ParMesh *mesh = vfes->GetParMesh(); + const FiniteElementCollection *fec = vfes->FEColl(); + ParFiniteElementSpace dfes(mesh, fec, dim, Ordering::byNODES); + ParFiniteElementSpace fes(mesh, fec); + + sigma = new ParGridFunction(&fes); + *sigma = 0.; + double *hSigma = sigma->HostWrite(); + + ParGridFunction coords(&dfes); + mesh->GetNodes(coords); + int ndofs = vfes->GetNDofs(); + + vector nodesVec; + nodesVec.clear(); + for (int n = 0; n < ndofs; n++) { + Vector Xn(dim); + for (int d = 0; d < dim; d++) Xn[d] = coords[n + d * ndofs]; + + // distance to the mix-out plane + double dist = 0.; + for (int d = 0; d < dim; d++) dist += szData.normal[d] * (Xn[d] - szData.pointInit[d]); + + if (fabs(dist) < szData.tol) nodesVec.push_back(n); + + // dist end plane + double distF = 0.; + for (int d = 0; d < dim; d++) distF += szData.normal[d] * (Xn[d] - szData.point0[d]); + + if (dist < 0. && distF > 0.) { + double planeDistance = distF - dist; + hSigma[n] = -dist / planeDistance / planeDistance; + } + } + + // find plane nodes + if (szData.szType == SpongeZoneSolution::MIXEDOUT) { + nodesInMixedOutPlane.SetSize(nodesVec.size()); + for (int n = 0; n < nodesInMixedOutPlane.Size(); n++) nodesInMixedOutPlane[n] = nodesVec[n]; + + cout << "__________________________________________________" << endl; + cout << "Running with Mixed-Out Sponge Zone" << endl; + cout << " Nodes searched with tolerance: " << szData.tol << endl; + cout << " Num. of nodes in mix-out plane: " << nodesInMixedOutPlane.Size() << endl; + cout << "^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^" << endl; + } +} + +SpongeZone::~SpongeZone() { delete sigma; } + +void SpongeZone::updateTerms(Vector &in) { + if (szData.szType == SpongeZoneSolution::MIXEDOUT) computeMixedOutValues(); + + addSpongeZoneForcing(in); +} + +void SpongeZone::addSpongeZoneForcing(Vector &in) { + const double *ds = sigma->HostRead(); + double *dataIn = in.HostReadWrite(); + const double *dataUp = Up->HostRead(); + + int nnodes = vfes->GetNDofs(); + + Vector Un(num_equation), Up(num_equation); + + // compute speed of sound + double gamma = eqState->GetSpecificHeatRatio(); + eqState->GetPrimitivesFromConservatives(targetU, Up, dim, num_equation); + double speedSound = sqrt(gamma * Up[num_equation - 1] / Up[0]); + + // add forcing to RHS, i.e., @in + for (int n = 0; n < nnodes; n++) { + double s = ds[n]; + if (s > 0.) { + s *= szData.multFactor; + for (int eq = 0; eq < num_equation; eq++) Up[eq] = dataUp[n + eq * nnodes]; + eqState->GetConservativesFromPrimitives(Up, Un, dim, num_equation); + + for (int eq = 0; eq < num_equation; eq++) dataIn[n + eq * nnodes] -= speedSound * s * (Un[eq] - targetU[eq]); + } + } +} + +void SpongeZone::computeMixedOutValues() { + int nnodes = vfes->GetNDofs(); + const double *dataUp = Up->HostRead(); + double gamma = eqState->GetSpecificHeatRatio(); + + // compute mean normal fluxes + meanNormalFluxes = 0.; + + Vector Un(num_equation), Up(num_equation); + DenseMatrix f(num_equation, dim); + for (int n = 0; n < nodesInMixedOutPlane.Size(); n++) { + int node = nodesInMixedOutPlane[n]; + + for (int eq = 0; eq < num_equation; eq++) Up[eq] = dataUp[node + eq * nnodes]; + eqState->GetConservativesFromPrimitives(Up, Un, dim, num_equation); + fluxes->ComputeConvectiveFluxes(Un, f); + + for (int eq = 0; eq < num_equation; eq++) + for (int d = 0; d < dim; d++) meanNormalFluxes[eq] += szData.normal[d] * f(eq, d); + } + meanNormalFluxes[num_equation] = double(nodesInMixedOutPlane.Size()); + + // communicate different partitions + Vector sum(num_equation + 1); + MPI_Allreduce(meanNormalFluxes.HostRead(), sum.HostWrite(), num_equation + 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD); + meanNormalFluxes = sum; + + // normalize + for (int eq = 0; eq < num_equation; eq++) meanNormalFluxes[eq] /= double(meanNormalFluxes[num_equation]); + + // compute mixed-out variables considering the outlet pressure + double temp = 0.; + for (int d = 0; d < dim; d++) temp += meanNormalFluxes[1 + d] * szData.normal[d]; + double A = 1. - 2. * gamma / (gamma - 1.); + double B = 2 * temp / (gamma - 1.); + double C = -2. * meanNormalFluxes[0] * meanNormalFluxes[num_equation - 1]; + for (int d = 0; d < dim; d++) C += meanNormalFluxes[1 + d] * meanNormalFluxes[1 + d]; + // double p = (-B+sqrt(B*B-4.*A*C))/(2.*A); + double p = (-B - sqrt(B * B - 4. * A * C)) / (2. * A); // real solution + + Up[0] = meanNormalFluxes[0] * meanNormalFluxes[0] / (temp - p); + Up[num_equation - 1] = p; + + for (int d = 0; d < dim; d++) Up[1 + d] = (meanNormalFluxes[1 + d] - p * szData.normal[d]) / meanNormalFluxes[0]; + + double v2 = 0.; + for (int d = 0; d < dim; d++) v2 += Up[1 + d] * Up[1 + d]; + + eqState->GetConservativesFromPrimitives(Up, targetU, dim, num_equation); +} + #ifdef _MASA_ MASA_forcings::MASA_forcings(const int &_dim, const int &_num_equation, const int &_order, const int &_intRuleType, IntegrationRules *_intRules, ParFiniteElementSpace *_vfes, ParGridFunction *_Up, diff --git a/src/forcing_terms.hpp b/src/forcing_terms.hpp index 00e0dde4c..42b359560 100644 --- a/src/forcing_terms.hpp +++ b/src/forcing_terms.hpp @@ -37,6 +37,8 @@ #include #include "dataStructures.hpp" +#include "equation_of_state.hpp" +#include "fluxes.hpp" #include "run_configuration.hpp" #ifdef _MASA_ @@ -100,6 +102,32 @@ class ConstantPressureGradient : public ForcingTerms { #endif }; +class SpongeZone : public ForcingTerms { + private: + Fluxes *fluxes; + EquationOfState *eqState; + + SpongeZoneData &szData; + Vector targetU; + + Array nodesInMixedOutPlane; + + ParGridFunction *sigma; // linearly varying factor + + Vector meanNormalFluxes; + + void computeMixedOutValues(); + void addSpongeZoneForcing(Vector &in); + + public: + SpongeZone(const int &_dim, const int &_num_equation, const int &_order, const int &_intRuleType, Fluxes *_fluxClass, + EquationOfState *_eqState, IntegrationRules *_intRules, ParFiniteElementSpace *_vfes, ParGridFunction *_Up, + ParGridFunction *_gradUp, const volumeFaceIntegrationArrays &gpuArrays, RunConfiguration &_config); + ~SpongeZone(); + + virtual void updateTerms(Vector &in); +}; + #ifdef _MASA_ // Manufactured Solution using MASA class MASA_forcings : public ForcingTerms { diff --git a/src/gradients.hpp b/src/gradients.hpp index 814a976bf..e5141e6fa 100644 --- a/src/gradients.hpp +++ b/src/gradients.hpp @@ -35,6 +35,7 @@ // Class to manage gradients of primitives #include + #include #include diff --git a/src/main.cpp b/src/main.cpp index 4569e4905..7baa43817 100644 --- a/src/main.cpp +++ b/src/main.cpp @@ -29,14 +29,16 @@ // OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE // OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. // -----------------------------------------------------------------------------------el- +#include + #include #include -#include -#include #include -#include "tps.hpp" +#include + #include "em_options.hpp" #include "quasimagnetostatic.hpp" +#include "tps.hpp" int main(int argc, char *argv[]) { MPI_Session mpi(argc, argv); @@ -52,12 +54,11 @@ int main(int argc, char *argv[]) { OptionsParser args(argc, argv); args.AddOption(&flow_only, "-flow", "--flow-only", "-nflow", "--not-flow-only", "Perform flow only simulation"); args.AddOption(&inputFile, "-run", "--runFile", "Name of the input file with run options."); - args.AddOption(&showVersion, "-v", "--version", "" , "--no-version", "Print code version and exit"); + args.AddOption(&showVersion, "-v", "--version", "", "--no-version", "Print code version and exit"); // Add options for EM ElectromagneticOptions em_opt; - args.AddOption(&em_only, "-em", "--em-only", "-nem", "--not-em-only", - "Perform electromagnetics only simulation"); + args.AddOption(&em_only, "-em", "--em-only", "-nem", "--not-em-only", "Perform electromagnetics only simulation"); em_opt.AddElectromagneticOptions(args); // device_config inferred from build setup @@ -82,8 +83,7 @@ int main(int argc, char *argv[]) { // version info tps.PrintHeader(); - if (showVersion) - exit(0); + if (showVersion) exit(0); if (mpi.Root()) args.PrintOptions(cout); diff --git a/src/quasimagnetostatic.cpp b/src/quasimagnetostatic.cpp index 9845c525f..db1258fcf 100644 --- a/src/quasimagnetostatic.cpp +++ b/src/quasimagnetostatic.cpp @@ -30,12 +30,14 @@ // OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. // -----------------------------------------------------------------------------------el- +#include "quasimagnetostatic.hpp" + #include #include + #include "../utils/mfem_extras/pfem_extras.hpp" #include "logger.hpp" #include "utils.hpp" -#include "quasimagnetostatic.hpp" using namespace mfem; using namespace mfem::common; @@ -43,27 +45,25 @@ using namespace mfem::common; void JFun(const Vector &x, Vector &f); QuasiMagnetostaticSolver::QuasiMagnetostaticSolver(MPI_Session &mpi, ElectromagneticOptions em_opts) - : - _mpi(mpi), - _em_opts(em_opts) { + : _mpi(mpi), _em_opts(em_opts) { // dump options to screen for user inspection if (_mpi.Root()) { _em_opts.print(std::cout); } - _pmesh = NULL; - _hcurl = NULL; - _h1 = NULL; - _hdiv = NULL; + _pmesh = NULL; + _hcurl = NULL; + _h1 = NULL; + _hdiv = NULL; _Aspace = NULL; _pspace = NULL; _Bspace = NULL; - _K = NULL; - _A = NULL; - _B = NULL; + _K = NULL; + _A = NULL; + _B = NULL; _operator_initialized = false; - _current_initialized = false; + _current_initialized = false; } QuasiMagnetostaticSolver::~QuasiMagnetostaticSolver() { @@ -116,12 +116,12 @@ void QuasiMagnetostaticSolver::Initialize() { // 2) Prepare the required finite elements //----------------------------------------------------- _hcurl = new ND_FECollection(_em_opts.order, _dim); - _h1 = new H1_FECollection(_em_opts.order, _dim); - _hdiv = new RT_FECollection(_em_opts.order-1, _dim); + _h1 = new H1_FECollection(_em_opts.order, _dim); + _hdiv = new RT_FECollection(_em_opts.order - 1, _dim); _Aspace = new ParFiniteElementSpace(_pmesh, _hcurl); - _pspace = new ParFiniteElementSpace(_pmesh, _h1 ); - _Bspace = new ParFiniteElementSpace(_pmesh, _hdiv ); + _pspace = new ParFiniteElementSpace(_pmesh, _h1); + _Bspace = new ParFiniteElementSpace(_pmesh, _hdiv); //----------------------------------------------------- // 3) Get BC dofs (everything is essential---i.e., PEC) @@ -184,18 +184,17 @@ void QuasiMagnetostaticSolver::InitializeCurrent() { // 2) Build a discretely divergence-free approximation of the source // current that lives in the Nedelec FE space defined in // Initialize() - ParGridFunction* Jorig = new ParGridFunction(_Aspace); - ParGridFunction* Jproj = new ParGridFunction(_Aspace); + ParGridFunction *Jorig = new ParGridFunction(_Aspace); + ParGridFunction *Jproj = new ParGridFunction(_Aspace); _r = new ParLinearForm(_Aspace); int irOrder = _pspace->GetElementTransformation(0)->OrderW() + 2 * _em_opts.order; - ParDiscreteGradOperator* grad = new ParDiscreteGradOperator(_pspace, _Aspace); + ParDiscreteGradOperator *grad = new ParDiscreteGradOperator(_pspace, _Aspace); grad->Assemble(); grad->Finalize(); - DivergenceFreeProjector* div_free = - new DivergenceFreeProjector(*_pspace, *_Aspace, irOrder, NULL, NULL, grad); + DivergenceFreeProjector *div_free = new DivergenceFreeProjector(*_pspace, *_Aspace, irOrder, NULL, NULL, grad); // This call (i.e., GlobalProjectDiscCoefficient) replaces the // functionality of Jorig->ProjectCoefficient(current) in a way that @@ -218,7 +217,7 @@ void QuasiMagnetostaticSolver::InitializeCurrent() { delete Jorig; // 3) Multiply by the mass matrix to get the RHS vector - ParBilinearForm* mass = new ParBilinearForm(_Aspace); + ParBilinearForm *mass = new ParBilinearForm(_Aspace); mass->AddDomainIntegrator(new VectorFEMassIntegrator); mass->Assemble(); mass->Finalize(); @@ -239,7 +238,6 @@ void QuasiMagnetostaticSolver::Solve() { // Solve for the magnetic vector potential and magnetic field in 3 steps - // 1) Solve the curl(curl(A)) = J for magnetic vector potential // using operators set up by Initialize() and InitializeCurrent() // fcns @@ -309,11 +307,11 @@ void QuasiMagnetostaticSolver::InterpolateToYAxis() const { // Set up array of points (on all ranks) DenseMatrix phys_points(_dim, _em_opts.nBy); - const double dy = (_em_opts.yinterp_max - _em_opts.yinterp_min)/(_em_opts.nBy-1); + const double dy = (_em_opts.yinterp_max - _em_opts.yinterp_min) / (_em_opts.nBy - 1); for (int ipt = 0; ipt < _em_opts.nBy; ipt++) { phys_points(0, ipt) = 0.0; - phys_points(1, ipt) = _em_opts.yinterp_min + ipt*dy; + phys_points(1, ipt) = _em_opts.yinterp_min + ipt * dy; phys_points(2, ipt) = 0.0; } @@ -395,16 +393,18 @@ void QuasiMagnetostaticSolver::InterpolateToYAxis() const { delete By; } -void JFun(const Vector & x, Vector & J) { - Vector a(3); a(0) = 0.0; a(1) = 1.0; a(2) = 0.0; +void JFun(const Vector &x, Vector &J) { + Vector a(3); + a(0) = 0.0; + a(1) = 1.0; + a(2) = 0.0; Vector axx(3); - axx(0) = a(1)*x(2) - a(2)*x(1); - axx(1) = a(2)*x(0) - a(0)*x(2); - axx(2) = a(0)*x(1) - a(1)*x(0); + axx(0) = a(1) * x(2) - a(2) * x(1); + axx(1) = a(2) * x(0) - a(0) * x(2); + axx(2) = a(0) * x(1) - a(1) * x(0); axx /= axx.Norml2(); J = axx; } - diff --git a/src/quasimagnetostatic.hpp b/src/quasimagnetostatic.hpp index 543d794b6..78062a75a 100644 --- a/src/quasimagnetostatic.hpp +++ b/src/quasimagnetostatic.hpp @@ -34,6 +34,7 @@ #define QUASIMAGNETOSTATIC_HPP_ #include + #include #include diff --git a/src/rhs_operator.cpp b/src/rhs_operator.cpp index 9221b1f04..9ea505d29 100644 --- a/src/rhs_operator.cpp +++ b/src/rhs_operator.cpp @@ -88,6 +88,10 @@ RHSoperator::RHSoperator(int &_iter, const int _dim, const int &_num_equations, forcing.Append(new ConstantPressureGradient(dim, num_equation, _order, intRuleType, intRules, vfes, Up, gradUp, gpuArrays, _config)); } + if (_config.GetSpongeZoneData().szType != SpongeZoneSolution::NONE) { + forcing.Append(new SpongeZone(dim, num_equation, _order, intRuleType, fluxClass, eqState, intRules, vfes, Up, + gradUp, gpuArrays, _config)); + } #ifdef _MASA_ forcing.Append( new MASA_forcings(dim, num_equation, _order, intRuleType, intRules, vfes, Up, gradUp, gpuArrays, _config)); diff --git a/src/run_configuration.cpp b/src/run_configuration.cpp index 2072e0ed7..3a15b3e06 100644 --- a/src/run_configuration.cpp +++ b/src/run_configuration.cpp @@ -85,6 +85,8 @@ RunConfiguration::RunConfiguration() { linViscData.pointInit = 0.; linViscData.viscRatio = 0.; + initSpongeData(); + isForcing = false; for (int ii = 0; ii < 3; ii++) gradPress[ii] = 0.; @@ -97,6 +99,28 @@ RunConfiguration::RunConfiguration() { RunConfiguration::~RunConfiguration() {} +void RunConfiguration::initSpongeData() { + spongeData.multFactor = 1.; + + spongeData.normal.UseDevice(true); + spongeData.point0.UseDevice(true); + spongeData.pointInit.UseDevice(true); + spongeData.targetUp.UseDevice(true); + + spongeData.normal.SetSize(3); + spongeData.point0.SetSize(3); + spongeData.pointInit.SetSize(3); + spongeData.targetUp.SetSize(5); + + spongeData.normal = 0.; + spongeData.point0 = 0.; + spongeData.pointInit = 0.; + spongeData.targetUp = 0.; + + spongeData.tol = 1e-5; + spongeData.szType = SpongeZoneSolution::NONE; +} + void RunConfiguration::readInputFile(std::string inpuFileName) { // runfile ifstream runFile(inpuFileName); @@ -330,7 +354,7 @@ void RunConfiguration::readInputFile(std::string inpuFileName) { gradPress[2] = stod(word); isForcing = true; - } else if (word.compare("PLANE_NORM") == 0) { + } else if (word.compare("LV_PLANE_NORM") == 0) { auto normal = linViscData.normal.HostWrite(); ss >> word; normal[0] = stod(word); @@ -344,7 +368,7 @@ void RunConfiguration::readInputFile(std::string inpuFileName) { modulus = sqrt(modulus); for (int d = 0; d < 3; d++) normal[d] /= modulus; - } else if (word.compare("PLANE_P0") == 0) { + } else if (word.compare("LV_PLANE_P0") == 0) { auto point0 = linViscData.point0.HostWrite(); ss >> word; point0[0] = stod(word); @@ -353,7 +377,7 @@ void RunConfiguration::readInputFile(std::string inpuFileName) { ss >> word; point0[2] = stod(word); - } else if (word.compare("PLANE_PINIT") == 0) { + } else if (word.compare("LV_PLANE_PINIT") == 0) { auto pointInit = linViscData.pointInit.HostWrite(); ss >> word; pointInit[0] = stod(word); @@ -362,10 +386,67 @@ void RunConfiguration::readInputFile(std::string inpuFileName) { ss >> word; pointInit[2] = stod(word); - } else if (word.compare("VISC_RATIO") == 0) { + } else if (word.compare("LV_VISC_RATIO") == 0) { ss >> word; linViscData.viscRatio = stod(word); + } else if (word.compare("SZ_PLANE_NORM") == 0) { + auto hnorm = spongeData.normal.HostWrite(); + ss >> word; + hnorm[0] = stod(word); + ss >> word; + hnorm[1] = stod(word); + ss >> word; + hnorm[2] = stod(word); + + } else if (word.compare("SZ_PLANE_P0") == 0) { + auto hp0 = spongeData.point0.HostWrite(); + ss >> word; + hp0[0] = stod(word); + ss >> word; + hp0[1] = stod(word); + ss >> word; + hp0[2] = stod(word); + + } else if (word.compare("SZ_PLANE_PINIT") == 0) { + auto hpi = spongeData.pointInit.HostWrite(); + ss >> word; + hpi[0] = stod(word); + ss >> word; + hpi[1] = stod(word); + ss >> word; + hpi[2] = stod(word); + + } else if (word.compare("SZ_TYPE") == 0) { + ss >> word; + + auto hup = spongeData.targetUp.HostWrite(); + + int szTypeInt = stoi(word); + switch (szTypeInt) { + case 0: + spongeData.szType = SpongeZoneSolution::USERDEF; + ss >> word; + hup[0] = stod(word); // rho + ss >> word; + hup[1] = stod(word); // u + ss >> word; + hup[2] = stod(word); // v + ss >> word; + hup[3] = stod(word); // w + ss >> word; + hup[4] = stod(word); // p + break; + case 1: + spongeData.szType = SpongeZoneSolution::MIXEDOUT; + ss >> word; + spongeData.tol = stod(word); + break; + } + } else if (word.compare("SZ_MULT") == 0) { + ss >> word; + spongeData.multFactor = stod(word); + } else if (word.compare("INLET") == 0) { ss >> word; std::pair patchANDtype; diff --git a/src/run_configuration.hpp b/src/run_configuration.hpp index 0fde79319..aa79c4af1 100644 --- a/src/run_configuration.hpp +++ b/src/run_configuration.hpp @@ -139,6 +139,9 @@ class RunConfiguration { // the plane defined by point0. linearlyVaryingVisc linViscData; + SpongeZoneData spongeData; + void initSpongeData(); + // Reference length double refLength; @@ -212,6 +215,7 @@ class RunConfiguration { double* GetConstantInitialCondition() { return &initRhoRhoVp[0]; } linearlyVaryingVisc& GetLinearVaryingData() { return linViscData; } + SpongeZoneData& GetSpongeZoneData() { return spongeData; } // resource manager controls bool isAutoRestart() { return rm_enableMonitor_; } diff --git a/src/runfileExample.run b/src/runfileExample.run index 43ffdcc11..87315b4e2 100644 --- a/src/runfileExample.run +++ b/src/runfileExample.run @@ -127,10 +127,22 @@ GRAD_PRESSURE 8 0 0 # from 0 to VISC_RATIO from the plane defined by PLANE_PINIT and # the plane defined by PLANE_P0. This is intended to help in # problematic boundary conditions. -PLANE_NORM n1 n2 n3 -PLANE_P0 x0 y0 z0 -PLANE_PINIT xi yi zi -VISC_RATIO 50. +# |-----------------------------------| +# |------------------------------<====| +# |-----------------------------------| +# p_init normal p0 +LV_PLANE_NORM n1 n2 n3 +LV_PLANE_P0 x0 y0 z0 +LV_PLANE_PINIT xi yi zi +LV_VISC_RATIO 50. + +# Sponge zone input. See dataStructures.hpp for more info +SZ_PLANE_NORM n1 n2 n3 +SZ_PLANE_P0 x0 y0 z0 +SZ_PLANE_PINIT xi yi zi +SZ_TYPE 0 rho u v w p +#SZ_TYPE 1 tol +#SZ_MULT 0.5 # Inlet Boundary Conditions # types: diff --git a/src/tps.cpp b/src/tps.cpp index a26ea6257..a8450936a 100644 --- a/src/tps.cpp +++ b/src/tps.cpp @@ -33,7 +33,7 @@ tps::tps(MPI_Session &mpi, int &argc, char **&argv) { nprocs_ = mpi.WorldSize(); - rank_ = mpi.WorldRank(); + rank_ = mpi.WorldRank(); if (rank_ == 0) isRank0_ = true; else diff --git a/src/tps.hpp b/src/tps.hpp index c71b3bb91..97c7ae3f6 100644 --- a/src/tps.hpp +++ b/src/tps.hpp @@ -33,6 +33,7 @@ #define TPS_HPP_ #include + #include "M2ulPhyS.hpp" // Solver parent class shared by all implementations diff --git a/src/utils.cpp b/src/utils.cpp index 287981295..3da1e9157 100644 --- a/src/utils.cpp +++ b/src/utils.cpp @@ -29,18 +29,21 @@ // OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE // OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. // -----------------------------------------------------------------------------------el- +#include "utils.hpp" + #include #include + #include #include #include #include #include +#include #include #include -#include + #include "M2ulPhyS.hpp" -#include "utils.hpp" #ifdef HAVE_SLURM #include @@ -178,9 +181,7 @@ std::string systemCmd(const char *cmd) { return result; } -void LocalProjectDiscCoefficient(GridFunction &gf, - VectorCoefficient &coeff, - Array &dof_attr) { +void LocalProjectDiscCoefficient(GridFunction &gf, VectorCoefficient &coeff, Array &dof_attr) { Array vdofs; Vector vals; @@ -213,7 +214,7 @@ void LocalProjectDiscCoefficient(GridFunction &gf, // Mimic how negative indices are handled in Vector::SetSubVector int ind = vdofs[j]; if (ind < 0) { - ind = -1-ind; + ind = -1 - ind; vals[j] = -vals[j]; } @@ -225,8 +226,7 @@ void LocalProjectDiscCoefficient(GridFunction &gf, } } -void GlobalProjectDiscCoefficient(ParGridFunction &gf, - VectorCoefficient &coeff) { +void GlobalProjectDiscCoefficient(ParGridFunction &gf, VectorCoefficient &coeff) { // local maximal element attribute for each dof Array ldof_attr; diff --git a/src/utils.hpp b/src/utils.hpp index 1b54f8916..9a87dfefa 100644 --- a/src/utils.hpp +++ b/src/utils.hpp @@ -32,22 +32,24 @@ #ifndef UTILS_HPP_ #define UTILS_HPP_ -#include #include -#include +#include + #include +#include // Misc. utilities bool file_exists(const std::string &name); std::string systemCmd(const char *cmd); // HDF5 convenience utilities -inline hid_t h5_getType( int) { return(H5T_NATIVE_INT); } -inline hid_t h5_getType(double) { return(H5T_NATIVE_DOUBLE); } +inline hid_t h5_getType(int) { return (H5T_NATIVE_INT); } +inline hid_t h5_getType(double) { return (H5T_NATIVE_DOUBLE); } -template void h5_save_attribute(hid_t dest, std::string attribute, T value) { +template +void h5_save_attribute(hid_t dest, std::string attribute, T value) { hid_t attr, status; - hid_t attrType = h5_getType(value); + hid_t attrType = h5_getType(value); hid_t dataspaceId = H5Screate(H5S_SCALAR); attr = H5Acreate(dest, attribute.c_str(), attrType, dataspaceId, H5P_DEFAULT, H5P_DEFAULT); @@ -58,13 +60,14 @@ template void h5_save_attribute(hid_t dest, std::string attribute, H5Sclose(dataspaceId); } -template void h5_read_attribute(hid_t source, std::string attribute, T &value) { +template +void h5_read_attribute(hid_t source, std::string attribute, T &value) { herr_t status; hid_t attr; hid_t attrType = h5_getType(value); // T temp; - attr = H5Aopen_name(source, attribute.c_str()); + attr = H5Aopen_name(source, attribute.c_str()); status = H5Aread(attr, attrType, &value); assert(status >= 0); H5Aclose(attr); @@ -83,9 +86,7 @@ template void h5_read_attribute(hid_t source, std::string attribute * mfem::GridFunction::ProjectDiscCoefficient but has fixes s.t. it * will work with Nedelec elements. */ -void LocalProjectDiscCoefficient(mfem::GridFunction &gf, - mfem::VectorCoefficient &coeff, - mfem::Array &dof_attr); +void LocalProjectDiscCoefficient(mfem::GridFunction &gf, mfem::VectorCoefficient &coeff, mfem::Array &dof_attr); /** Project discontinous function onto FE space * @@ -97,7 +98,6 @@ void LocalProjectDiscCoefficient(mfem::GridFunction &gf, * mfem::ParGridFunction::ProjectDiscCoefficient but calls * LocalProjectDiscCoefficient. */ -void GlobalProjectDiscCoefficient(mfem::ParGridFunction &gf, - mfem::VectorCoefficient &coeff); +void GlobalProjectDiscCoefficient(mfem::ParGridFunction &gf, mfem::VectorCoefficient &coeff); #endif // UTILS_HPP_ diff --git a/test/Makefile.am b/test/Makefile.am index ab234db79..a5a98db50 100644 --- a/test/Makefile.am +++ b/test/Makefile.am @@ -4,7 +4,7 @@ EXTRA_DIST = tap-driver.sh soln_differ inputs meshes ref_solns/*.h5 \ vpath.sh die.sh cyl3d.gpu.test cyl3d.mflow.gpu.test wedge.gpu.test \ averages.gpu.test cyl3d.test cyl3d.mflow.test cyl3d.dtconst.test \ die.test averages.test wedge.test cyl3d.interp.test qms.rings.test \ - valgrind.test + sponge_zone.test valgrind.test TESTS = vpath.sh @@ -20,7 +20,8 @@ TESTS += cyl3d.test \ die.test \ averages.test \ wedge.test \ - qms.rings.test + qms.rings.test \ + sponge_zone.test endif if GSLIB_ENABLED diff --git a/test/inputs/input.sponge_zone.periodic b/test/inputs/input.sponge_zone.periodic new file mode 100644 index 000000000..0c1c7060b --- /dev/null +++ b/test/inputs/input.sponge_zone.periodic @@ -0,0 +1,38 @@ +MESH meshes/periodic-square.mesh + +OUTPUT_NAME sponge_zone + +POL_ORDER 1 + +INT_RULE 1 + +BASIS_TYPE 1 + +CFL 0.12 +#DT_CONSTANT + +NMAX 3000 +ITERS_OUT 100 + +IS_SBP 0 + +TIME_INTEGRATOR 4 + +FLUID 0 + +#VISC_MULT . + +EQ_SYSTEM 1 + +INIT_RHO 1.2 +INIT_RHOVX 100. +INIT_RHOVY 100 +INIT_RHOVZ 0. +INIT_P 101300 + +SZ_PLANE_NORM -1 0 0 +SZ_PLANE_P0 1 0 0 +SZ_PLANE_PINIT -1 0 0 +SZ_TYPE 0 1.2 100 100 0 101300 +SZ_MULT 1 + diff --git a/test/meshes/periodic-square.mesh b/test/meshes/periodic-square.mesh new file mode 100644 index 000000000..5739c7d3e --- /dev/null +++ b/test/meshes/periodic-square.mesh @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:b8419c487e65241336c6d4e83a4829588d40c2951c742347c2527335d903ae11 +size 1182 diff --git a/test/ref_solns/forcing-periodic-square.h5 b/test/ref_solns/forcing-periodic-square.h5 new file mode 100644 index 000000000..07029c620 --- /dev/null +++ b/test/ref_solns/forcing-periodic-square.h5 @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:0bf8372514bd29f11888551e76a82b97d3f57f0eba6045af385046021da41ea4 +size 5768 diff --git a/test/sponge_zone.test b/test/sponge_zone.test new file mode 100755 index 000000000..f9ce57c34 --- /dev/null +++ b/test/sponge_zone.test @@ -0,0 +1,28 @@ +#!./bats +# -*- mode: sh -*- + +TEST="sponge_zone.test" +RUNFILE="inputs/input.sponge_zone.periodic" +SOLN_FILE=restart_sponge_zone.sol.h5 +REF_FILE=ref_solns/forcing-periodic-square.h5 + +#setup() { +#} + +@test "[$TEST] check for input file $RUNFILE" { + test -s $RUNFILE +} + +@test "[$TEST] check for reference solution -> $REF_FILE" { + test -s $REF_FILE +} + + +@test "[$TEST] verify sponge-zone output -> $SOLN_FILE" { + rm -f $SOLN_FILE + + run ../src/tps --runFile $RUNFILE + + test -s $SOLN_FILE + ./soln_differ -d 2 $SOLN_FILE $REF_FILE +}