diff --git a/src/M2ulPhyS.hpp b/src/M2ulPhyS.hpp index ebbc9d91..6f9b67b0 100644 --- a/src/M2ulPhyS.hpp +++ b/src/M2ulPhyS.hpp @@ -295,9 +295,6 @@ class M2ulPhyS : public TPS::PlasmaSolver { // for use if we want to write a serial file Mesh *serial_mesh; - // I/O organizer - IODataOrganizer ioData; - std::list tableHost_; #ifdef HAVE_MASA @@ -470,6 +467,8 @@ class M2ulPhyS : public TPS::PlasmaSolver { int getMaximumIterations() const { return MaxIters; } int getCurrentIterations() const { return iter; } void setMaximumIterations(int value) { MaxIters = value; } + + double getCurrentTime() const { return time; } }; #endif // M2ULPHYS_HPP_ diff --git a/src/io.hpp b/src/io.hpp index 4685132b..ea1b3b76 100644 --- a/src/io.hpp +++ b/src/io.hpp @@ -217,6 +217,9 @@ class IOFamily { * @param read_order Polynomial order of function in HDF5 file */ void readChangeOrder(hid_t file, int read_order); + + /** @brief Returns pointer to ParGridFunction for this IOFamily */ + mfem::ParGridFunction *getParGridFunction() { return pfunc_; } }; /** @@ -312,6 +315,13 @@ class IODataOrganizer { * @param read_order Polynomial order for data in file (optional, only required for variable order read) */ void read(hid_t file, bool serial, int read_order = -1); + + /** + * @brief Get reference to vector of IOFamily objects + * + * This supports operations on all IO fields (e.g., interpolation between grids) + */ + std::vector &getIOFamilies() { return families_; } }; /** diff --git a/src/loMach.hpp b/src/loMach.hpp index ba5a9066..9726b157 100644 --- a/src/loMach.hpp +++ b/src/loMach.hpp @@ -195,9 +195,8 @@ class LoMachSolver : public TPS::PlasmaSolver { StopWatch sw_setup_, sw_step_, sw_turb_, sw_thermChem_, sw_flow_, sw_press_; double tlast_; - // I/O helpers + // Viz helper ParaViewDataCollection *pvdc_ = nullptr; // visualization - IODataOrganizer ioData; // restart /// Update the EXTk/BDF time integration coefficient. void SetTimeIntegrationCoefficients(int step); diff --git a/src/lte_thermo_chem.cpp b/src/lte_thermo_chem.cpp index 30ae287e..901c87f7 100644 --- a/src/lte_thermo_chem.cpp +++ b/src/lte_thermo_chem.cpp @@ -172,6 +172,7 @@ LteThermoChem::LteThermoChem(mfem::ParMesh *pmesh, LoMachOptions *loMach_opts, t LteThermoChem::~LteThermoChem() { // allocated in initializeOperators + delete rad_kap_gradT_coeff_; delete rad_radiation_sink_coeff_; delete rad_jh_coeff_; delete rad_un_next_coeff_; @@ -183,29 +184,40 @@ LteThermoChem::~LteThermoChem() { delete sfes_filter_; delete sfec_filter_; + delete MqInv_; + delete MqInvPC_; delete HtInv_; delete HtInvPC_; delete MsInv_; delete MsInvPC_; delete MrhoInv_; delete MrhoInvPC_; + delete jh_form_; + delete Mq_form_; delete Ht_form_; delete M_rho_form_; delete M_rho_Cp_form_; delete Ms_form_; + delete A_rho_form_; delete At_form_; + delete LQ_form_; + delete LQ_bdry_; delete rho_Cp_u_coeff_; delete un_next_coeff_; delete kap_gradT_coeff_; delete gradT_coeff_; delete mut_coeff_; + delete kapt_coeff_; delete mult_coeff_; delete thermal_diff_coeff_; delete thermal_diff_sum_coeff_; delete thermal_diff_total_coeff_; delete rho_Cp_over_dt_coeff_; delete rho_Cp_coeff_; + delete Cp_coeff_; delete rho_coeff_; + delete jh_coeff_; + delete radiation_sink_coeff_; delete umag_coeff_; delete gscale_coeff_; @@ -224,6 +236,14 @@ LteThermoChem::~LteThermoChem() { // allocated in initializeSelf delete sfes_; delete sfec_; + + delete radiation_; + + delete mu_table_; + delete kappa_table_; + delete sigma_table_; + delete Rgas_table_; + delete Cp_table_; } void LteThermoChem::initializeSelf() { @@ -550,7 +570,7 @@ void LteThermoChem::initializeOperators() { kap_gradT_coeff_ = new ScalarVectorProductCoefficient(*thermal_diff_total_coeff_, *gradT_coeff_); un_next_coeff_ = new VectorGridFunctionCoefficient(flow_interface_->velocity); - rhon_next_coeff_ = new GridFunctionCoefficient(&rn_gf_); + // rhon_next_coeff_ = new GridFunctionCoefficient(&rn_gf_); rho_Cp_u_coeff_ = new ScalarVectorProductCoefficient(*rho_Cp_coeff_, *un_next_coeff_); diff --git a/src/solver.hpp b/src/solver.hpp index 87c0923d..6f396287 100644 --- a/src/solver.hpp +++ b/src/solver.hpp @@ -32,6 +32,7 @@ #ifndef SOLVER_HPP_ #define SOLVER_HPP_ +#include "io.hpp" #include "tps_mfem_wrap.hpp" #include "utils.hpp" @@ -119,6 +120,9 @@ class Solver { // Base class for flow solver implementations class PlasmaSolver : public Solver { + protected: + IODataOrganizer ioData; + public: virtual ~PlasmaSolver() {} @@ -141,6 +145,14 @@ class PlasmaSolver : public Solver { exit(1); return nullptr; } + + virtual void restart_files_hdf5(string mode, string inputFileName = std::string()) { + cout << "ERROR: " << __func__ << " remains unimplemented" << endl; + exit(1); + } + + /// Fetch reference to IODataOrganizer object + IODataOrganizer &getIODataOrganizer() { return ioData; } }; } // end namespace TPS diff --git a/utils/pfield_interpolate.cpp b/utils/pfield_interpolate.cpp index 99952894..d312175a 100644 --- a/utils/pfield_interpolate.cpp +++ b/utils/pfield_interpolate.cpp @@ -1,155 +1,301 @@ /* A utility to transfer tps solutions from one mesh to another */ -#include "tps.hpp" +#include "loMach.hpp" #include "M2ulPhyS.hpp" #include "mfem.hpp" +#include "tps.hpp" #include using namespace mfem; using namespace std; -int main (int argc, char *argv[]) -{ +int main(int argc, char *argv[]) { mfem::Mpi::Init(argc, argv); TPS::Tps tps(MPI_COMM_WORLD); - + TPS::Tps tar_tps(MPI_COMM_WORLD); // Set the method's default parameters. - const char *src_input_file = "coarse.run"; - const char *tar_input_file = "fine.run"; + const char *src_input_file = "coarse.ini"; + const char *tar_input_file = ""; + const char *tar_mesh_h1 = ""; + const char *pv_output_dir = "pv_output"; + int order = 1; // Parse command-line options. OptionsParser args(argc, argv); args.AddOption(&src_input_file, "-r1", "--runFile1", - "runFile for source case"); + "TPS input file for source case---i.e., the original mesh and " + "solution (required)"); args.AddOption(&tar_input_file, "-r2", "--runFile2", - "runFile for target case"); + "TPS input file for target case (optional--must use -r2 or " + "-mh1, but not both)"); + args.AddOption(&tar_mesh_h1, "-mh1", "--mesh-h1", + "Mesh to interpolate to (using H1 space) (optional--must use " + "-mh1 or -r2, but not both)"); + args.AddOption(&pv_output_dir, "-out", "--paraview-output", + "Paraview output directory (-mh1 only)"); + args.AddOption(&order, "-p", "--order-h1", "FEM order (H1 only)"); args.Parse(); - if (!args.Good()) - { - args.PrintUsage(cout); - return 1; - } + if (!args.Good()) { + args.PrintUsage(cout); + return 1; + } args.PrintOptions(cout); - // Instantiate M2ulPhyS classes for *both* the coarse and fine - - // NB: the M2ulPhyS ctor calls M2ulPhyS::initVariables, which reads - // the restart files, assuming that RESTART_CYCLE is set in the - // input file. So, we require that RESTART_CYCLE is set in the - // *source* run file.... string srcFileName(src_input_file); + string tarFileName(tar_input_file); + string tarMeshH1(tar_mesh_h1); + + // We have to have set one of these... + if (tarFileName.empty() && tarMeshH1.empty()) { + args.PrintUsage(cout); + return 1; + } + // but not both + if (!tarFileName.empty() && !tarMeshH1.empty()) { + args.PrintUsage(cout); + return 1; + } + + // Instantiate M2ulPhyS class for the "coarse" (i.e., source) case + + // Note that the M2ulPhyS ctor is responsible for reading the + // restart file, assuming that io/enableRestart = True in the tps + // input file. So, you need to set this option if you want to + // interpolate from an existing restart file. Otherwise, the + // M2ulPhyS ctor will initialize the field (to uniform by default) + // and then the code below will interpolate using that field. tps.parseInputFile(srcFileName); tps.chooseDevices(); - M2ulPhyS srcField( srcFileName, &tps ); - RunConfiguration& srcConfig = srcField.GetConfig(); - assert(srcConfig.GetRestartCycle()>0); - + const std::string &src_solver_type = tps.getSolverType(); + + TPS::PlasmaSolver *srcField = nullptr; + if (src_solver_type == "flow") { + srcField = new M2ulPhyS(srcFileName, &tps); + //RunConfiguration &srcConfig = srcField->GetConfig(); + //assert(srcConfig.GetRestartCycle() > 0); + } else if (src_solver_type == "loMach") { + srcField = new LoMachSolver(&tps); + srcField->parseSolverOptions(); + srcField->initialize(); + } tps.closeInputFile(); - ParMesh* mesh_1 = srcField.getMesh(); + ParMesh *mesh_1 = srcField->getMesh(); const int dim = mesh_1->Dimension(); + // Instantiate the "fine" (i.e., target) mesh. + // This must be done via the tps input file if the goal is to + // generate a tps restart on a new mesh. Alternatively, if just a + // mesh is specified, using the option --mesh-h1, that mesh is read + // and decomposed here. + TPS::PlasmaSolver *tarField = nullptr; + ParMesh *mesh_2 = nullptr; + if (!tarFileName.empty()) { + // Instantiate M2ulPhyS class for the "fine" (i.e., target) case + // Note that the M2ulPhyS ctor is responsible for reading the + // restart file, assuming that io/enableRestart = True in the tps + // input file. So, assuming there is not an existing restart file + // for the target mesh, you should not set this option. If + // it is set, M2ulPhyS will try to read a file that doesn't exist + // and give an error. + tar_tps.parseInputFile(tarFileName); + + const std::string &tar_solver_type = tar_tps.getSolverType(); + assert(src_solver_type == tar_solver_type); + + if (src_solver_type == "flow") { + tarField = new M2ulPhyS(tarFileName, &tar_tps); + } else if (src_solver_type == "loMach") { + tarField = new LoMachSolver(&tar_tps); + tarField->parseSolverOptions(); + tarField->initialize(); + } + tar_tps.closeInputFile(); - // but, we require that RESTART_CYCLE is *not* set in the - // target run file, since the target restart files do not exist yet. - string tarFileName(tar_input_file); - - tps.parseInputFile(tarFileName); - - M2ulPhyS tarField( tarFileName, &tps ); - RunConfiguration& tarConfig = tarField.GetConfig(); - assert(tarConfig.GetRestartCycle()==0); - - tps.closeInputFile(); + mesh_2 = tarField->getMesh(); + } else { + // Otherwise, read the mesh directly - // Get meshes - ParMesh* mesh_2 = tarField.getMesh(); + // All mpi ranks read the serial mesh + Mesh serial_mesh(tarMeshH1.c_str()); - //const int dim = mesh_1->Dimension(); + // Decompose + mesh_2 = new ParMesh(MPI_COMM_WORLD, serial_mesh); + } // GSLIB only works in 2 and 3 dimensions - assert(dim>1); + assert(dim > 1); // Input meshes must have same dimension - assert(mesh_2->Dimension()==dim); + assert(mesh_2->Dimension() == dim); - if (mesh_1->GetNodes() == NULL) { mesh_1->SetCurvature(1); } - if (mesh_2->GetNodes() == NULL) { mesh_2->SetCurvature(1); } + if (mesh_1->GetNodes() == NULL) { + mesh_1->SetCurvature(1); + } + if (mesh_2->GetNodes() == NULL) { + mesh_2->SetCurvature(1); + } - cout << "Source mesh curvature: " - << mesh_1->GetNodes()->OwnFEC()->Name() << endl - << "Target mesh curvature: " - << mesh_2->GetNodes()->OwnFEC()->Name() << endl; + cout << "Source mesh curvature: " << mesh_1->GetNodes()->OwnFEC()->Name() + << endl; + cout << "Target mesh curvature: " << mesh_2->GetNodes()->OwnFEC()->Name() + << endl; // 2) Set up source field const FiniteElementCollection *src_fec = NULL; ParGridFunction *func_source = NULL; - src_fec = srcField.getFEC(); - func_source = srcField.GetSolutionGF(); + src_fec = srcField->getFEC(); // 3) Some checks const Geometry::Type gt = mesh_2->GetNodalFESpace()->GetFE(0)->GetGeomType(); - MFEM_VERIFY(gt != Geometry::PRISM, "Wedge elements are not currently " - "supported."); - MFEM_VERIFY(mesh_2->GetNumGeometries(mesh_2->Dimension()) == 1, "Mixed meshes" - "are not currently supported."); + MFEM_VERIFY(gt != Geometry::PRISM, + "Wedge elements are not currently supported."); + MFEM_VERIFY(mesh_2->GetNumGeometries(mesh_2->Dimension()) == 1, + "Mixed meshes are not currently supported."); std::cout << "Source FE collection: " << src_fec->Name() << std::endl; - // Setup the FiniteElementSpace and GridFunction on the target mesh. - const FiniteElementCollection* tar_fec = tarField.getFEC(); - ParFiniteElementSpace* tar_fes = tarField.getFESpace(); - ParGridFunction* func_target = tarField.GetSolutionGF(); - std::cout << "Target FE collection: " << tar_fec->Name() << std::endl; - - const int NE = mesh_2->GetNE(); - const int nsp = tar_fes->GetFE(0)->GetNodes().GetNPoints(); - const int tar_ncomp = func_target->VectorDim(); - cout << "tar_ncomp = " << tar_ncomp << std::endl; - - // Generate list of points where the grid function will be evaluated. - Vector vxyz; - - // NB: We assume DG here! - vxyz.SetSize(nsp*NE*dim); - for (int i = 0; i < NE; i++) { - const FiniteElement *fe = tar_fes->GetFE(i); - const IntegrationRule ir = fe->GetNodes(); - ElementTransformation *et = tar_fes->GetElementTransformation(i); - - DenseMatrix pos; - et->Transform(ir, pos); - Vector rowx(vxyz.GetData() + i*nsp, nsp); - Vector rowy(vxyz.GetData() + i*nsp + NE*nsp, nsp); - Vector rowz; - if (dim == 3) { - rowz.SetDataAndSize(vxyz.GetData() + i*nsp + 2*NE*nsp, nsp); + // For every IOFamily that has been registered with the 'source', + // interpolate it to the 'target' mesh + std::vector& src_io_families = srcField->getIODataOrganizer().getIOFamilies(); + + for (size_t ifield = 0; ifield < src_io_families.size(); ifield++) { + //func_source = srcField->GetSolutionGF(); + func_source = src_io_families[ifield].getParGridFunction(); + + + // Setup the FiniteElementSpace and GridFunction on the target mesh. + const FiniteElementCollection *tar_fec = nullptr; + ParFiniteElementSpace *tar_fes = nullptr; + ParGridFunction *func_target = nullptr; + + if (!tarFileName.empty()) { + tar_fec = tarField->getFEC(); + tar_fes = tarField->getFESpace(); + std::vector& tar_io_families = tarField->getIODataOrganizer().getIOFamilies(); + func_target = tar_io_families[ifield].getParGridFunction(); + tar_fes = func_target->ParFESpace();// = tar_io_families[ifield].getParGridFunction(); + //func_target = tarField->GetSolutionGF(); + } else { + const int num_variables = func_source->VectorDim(); + tar_fec = new H1_FECollection(order, dim); + tar_fes = new ParFiniteElementSpace(mesh_2, tar_fec, num_variables); + func_target = new ParGridFunction(tar_fes); + } + + std::cout << "Target FE collection: " << tar_fec->Name() << std::endl; + + const int NE = mesh_2->GetNE(); + const int nsp = tar_fes->GetFE(0)->GetNodes().GetNPoints(); + const int tar_ncomp = func_target->VectorDim(); + cout << "tar_ncomp = " << tar_ncomp << std::endl; + + // Generate list of points where the grid function will be evaluated. + Vector vxyz; + + // The method for getting the necessary points here works for both + // H1 and L2 fields. However, for H1, it will contain duplicate + // points, which means we must take care when setting the + // corresponding ParGridFunction below. + vxyz.SetSize(nsp * NE * dim); + for (int i = 0; i < NE; i++) { + const FiniteElement *fe = tar_fes->GetFE(i); + const IntegrationRule ir = fe->GetNodes(); + ElementTransformation *et = tar_fes->GetElementTransformation(i); + + DenseMatrix pos; + et->Transform(ir, pos); + Vector rowx(vxyz.GetData() + i * nsp, nsp); + Vector rowy(vxyz.GetData() + i * nsp + NE * nsp, nsp); + Vector rowz; + if (dim == 3) { + rowz.SetDataAndSize(vxyz.GetData() + i * nsp + 2 * NE * nsp, nsp); + } + pos.GetRow(0, rowx); + pos.GetRow(1, rowy); + if (dim == 3) { + pos.GetRow(2, rowz); + } + } + + const int nodes_cnt = vxyz.Size() / dim; + + // Evaluate source grid function. + Vector interp_vals(nodes_cnt * tar_ncomp); + FindPointsGSLIB finder(MPI_COMM_WORLD); + finder.Setup(*mesh_1); + finder.Interpolate(vxyz, *func_source, interp_vals); + + // FIXME: Ensure this functions correctly for H1 + if (tar_fec->GetContType() == FiniteElementCollection::DISCONTINUOUS) { + func_target->SetFromTrueDofs(interp_vals); + } else if (tar_fec->GetContType() == FiniteElementCollection::CONTINUOUS) { + // Fill solution element-by-element + Array vdofs; + Vector elem_dof_vals(nsp * tar_ncomp); + + for (int i = 0; i < mesh_2->GetNE(); i++) { + tar_fes->GetElementVDofs(i, vdofs); + for (int j = 0; j < nsp; j++) { + for (int d = 0; d < tar_ncomp; d++) { + // Arrange values byNodes + int idx = d * nsp * NE + i * nsp + j; + elem_dof_vals(j + d * nsp) = interp_vals(idx); + } + } + func_target->SetSubVector(vdofs, elem_dof_vals); + } + func_target->SetFromTrueVector(); + } else { + std::cout << "FEC not understood" << std::endl; + exit(1); + } + + if (tarFileName.empty()) { + // Dump paraview for visualization + ParaViewDataCollection pvc(pv_output_dir, mesh_2); + pvc.SetLevelsOfDetail(order); + pvc.SetHighOrderOutput(true); + pvc.SetPrecision(8); + + pvc.SetCycle(0); + pvc.SetTime(0.0); + + ParFiniteElementSpace fes(mesh_2, tar_fec, 1); + int ndofs = fes.GetNDofs(); + for (int ivar = 0; ivar < tar_ncomp; ivar++) { + string U("U_"); + pvc.RegisterField(U + to_string(ivar), + new ParGridFunction(&fes, func_target->HostReadWrite() + + ivar * ndofs)); + } + pvc.Save(); + } + + // Free the internal gslib data. + finder.FreeData(); + + // If we own them, delete fe collection, etc + if (tarFileName.empty()) { + delete tar_fes; + delete tar_fec; + delete func_target; } - pos.GetRow(0, rowx); - pos.GetRow(1, rowy); - if (dim == 3) { pos.GetRow(2, rowz); } } - const int nodes_cnt = vxyz.Size() / dim; - // Evaluate source grid function. - Vector interp_vals(nodes_cnt*tar_ncomp); - FindPointsGSLIB finder(MPI_COMM_WORLD); - finder.Setup(*mesh_1); - finder.Interpolate(vxyz, *func_source, interp_vals); + if (!tarFileName.empty()) { + tarField->restart_files_hdf5("write"); + } - // Project the interpolated values to the target FiniteElementSpace. - // NB: Assume DG - func_target->SetFromTrueDofs(interp_vals); - // Write restart files - tarField.writeHDF5(); - // Free the internal gslib data. - finder.FreeData(); + // delete the target M2ulPhyS class + delete tarField; + delete srcField; return 0; }