From cd8508d3ceb5c44321889dbfcf430a17b699408d Mon Sep 17 00:00:00 2001 From: Seung Whan Chung Date: Thu, 24 Oct 2024 16:25:35 -0700 Subject: [PATCH 1/6] online poisson rom --- src/Control.cc | 3 + src/Potentials.cc | 15 ++--- src/Potentials.h | 14 ++++- src/read_config.cc | 2 + src/rom_Control.h | 3 + src/rom_workflows.cc | 144 +++++++++++++++++++++++++++++++++++++++++-- 6 files changed, 168 insertions(+), 13 deletions(-) diff --git a/src/Control.cc b/src/Control.cc index 3f8a11cc..cae4f08b 100644 --- a/src/Control.cc +++ b/src/Control.cc @@ -2065,6 +2065,8 @@ void Control::setROMOptions(const boost::program_options::variables_map& vm) rom_pri_option.num_orbbasis = vm["ROM.basis.number_of_orbital_basis"].as(); rom_pri_option.num_potbasis = vm["ROM.basis.number_of_potential_basis"].as(); rom_pri_option.pot_rom_file = vm["ROM.potential_rom_file"].as(); + + rom_pri_option.test_restart_file = vm["ROM.online.test_restart_file"].as(); } // onpe0 // synchronize all processors @@ -2081,6 +2083,7 @@ void Control::syncROMOptions() mmpi.bcast(rom_pri_option.restart_file_fmt, comm_global_); mmpi.bcast(rom_pri_option.basis_file, comm_global_); mmpi.bcast(rom_pri_option.pot_rom_file, comm_global_); + mmpi.bcast(rom_pri_option.test_restart_file, comm_global_); auto bcast_check = [](int mpirc) { if (mpirc != MPI_SUCCESS) diff --git a/src/Potentials.cc b/src/Potentials.cc index a54f0c47..8ec568a7 100644 --- a/src/Potentials.cc +++ b/src/Potentials.cc @@ -899,8 +899,9 @@ void Potentials::initBackground(Ions& ions) if (fabs(background_charge_) < 1.e-10) background_charge_ = 0.; } +#ifdef MGMOL_HAS_LIBROM void Potentials::evalIonDensityOnSamplePts( - Ions& ions, const std::vector &local_idx, std::vector &sampled_rhoc) + Ions& ions, const std::vector &local_idx, CAROM::Vector &sampled_rhoc) { Mesh* mymesh = Mesh::instance(); MGmol_MPI& mmpi = *(MGmol_MPI::instance()); @@ -918,9 +919,8 @@ void Potentials::evalIonDensityOnSamplePts( } // initialize output vector - sampled_rhoc.resize(local_idx.size()); - for (int d = 0; d < sampled_rhoc.size(); d++) - sampled_rhoc[d] = 0.0; + sampled_rhoc.setSize(local_idx.size()); + sampled_rhoc = 0.0; // Loop over ions for (auto& ion : ions.overlappingVL_ions()) @@ -936,9 +936,9 @@ void Potentials::evalIonDensityOnSamplePts( } void Potentials::initializeRadialDataOnSampledPts( - const Vector3D& position, const Species& sp, const std::vector &local_idx, std::vector &sampled_rhoc) + const Vector3D& position, const Species& sp, const std::vector &local_idx, CAROM::Vector &sampled_rhoc) { - assert(local_idx.size() == sampled_rhoc.size()); + assert(local_idx.size() == sampled_rhoc.dim()); Control& ct = *(Control::instance()); @@ -984,11 +984,12 @@ void Potentials::initializeRadialDataOnSampledPts( /* accumulate ion species density */ const double r = position.minimage(point, lattice, ct.bcPoisson); if (r < lrad) - sampled_rhoc[k] += sp.getRhoComp(r); + sampled_rhoc(k) += sp.getRhoComp(r); } return; } +#endif template void Potentials::setVxc( const double* const vxc, const int iterativeIndex); diff --git a/src/Potentials.h b/src/Potentials.h index 0ddd10fa..d580f7a0 100644 --- a/src/Potentials.h +++ b/src/Potentials.h @@ -18,6 +18,10 @@ #include #include +#ifdef MGMOL_HAS_LIBROM +#include "librom.h" +#endif + class Ions; class Species; template @@ -95,8 +99,10 @@ class Potentials void initializeSupersampledRadialDataOnMesh( const Vector3D& position, const Species& sp); +#ifdef MGMOL_HAS_LIBROM void initializeRadialDataOnSampledPts( - const Vector3D& position, const Species& sp, const std::vector &local_idx, std::vector &sampled_rhoc); + const Vector3D& position, const Species& sp, const std::vector &local_idx, CAROM::Vector &sampled_rhoc); +#endif public: Potentials(const bool vh_frozen = false); @@ -164,6 +170,8 @@ class Potentials const double getBackgroundCharge() const { return background_charge_; } + const double getIonicCharge() const { return ionic_charge_; } + /*! * initialize total potential as local pseudopotential */ @@ -201,7 +209,9 @@ class Potentials void initBackground(Ions& ions); void addBackgroundToRhoComp(); - void evalIonDensityOnSamplePts(Ions& ions, const std::vector &local_idx, std::vector &sampled_rhoc); +#ifdef MGMOL_HAS_LIBROM + void evalIonDensityOnSamplePts(Ions& ions, const std::vector &local_idx, CAROM::Vector &sampled_rhoc); +#endif #ifdef HAVE_TRICUBIC void readExternalPot(const string filename, const char type); diff --git a/src/read_config.cc b/src/read_config.cc index 884aa8d3..02c0aa23 100644 --- a/src/read_config.cc +++ b/src/read_config.cc @@ -438,5 +438,7 @@ void setupROMConfigOption(po::options_description &rom_cfg) "Number of potential POD basis to build Hartree potential ROM operator.") ("ROM.potential_rom_file", po::value()->default_value(""), "File name to save/load potential ROM operators."); + ("ROM.online.test_restart_file", po::value()->default_value(""), + "File name to test online ROM."); } #endif diff --git a/src/rom_Control.h b/src/rom_Control.h index 200061cb..5653267f 100644 --- a/src/rom_Control.h +++ b/src/rom_Control.h @@ -56,6 +56,9 @@ struct ROMPrivateOptions int num_orbbasis = -1; int num_potbasis = -1; std::string pot_rom_file = ""; + + /* options for online Poisson ROM */ + std::string test_restart_file = ""; }; #endif // ROM_CONTROL_H diff --git a/src/rom_workflows.cc b/src/rom_workflows.cc index 52ad6e47..49aea1d4 100644 --- a/src/rom_workflows.cc +++ b/src/rom_workflows.cc @@ -207,6 +207,13 @@ void buildROMPoissonOperator(MGmolInterface *mgmol_) std::vector global_sampled_row(num_pot_basis), sampled_rows_per_proc(nprocs); DEIM(pot_basis, num_pot_basis, global_sampled_row, sampled_rows_per_proc, pot_rhs_rom, rank, nprocs); + if (rank == 0) + { + int num_sample_rows = 0; + for (int k = 0; k < sampled_rows_per_proc.size(); k++) + num_sample_rows += sampled_rows_per_proc[k]; + printf("number of sampled row: %d\n", num_sample_rows); + } /* ROM rescaleTotalCharge operator */ CAROM::Vector fom_ones(pot_basis->numRows(), true); @@ -230,6 +237,14 @@ void buildROMPoissonOperator(MGmolInterface *mgmol_) h5_helper.putDoubleArray("potential_rom_inverse", pot_rom.getData(), num_pot_basis * num_pot_basis, false); + /* save right-hand side hyper-reduction sample index */ + h5_helper.putIntegerArray("potential_rhs_hr_idx", global_sampled_row.data(), + global_sampled_row.size(), false); + h5_helper.putIntegerArray("potential_rhs_hr_idcs_per_proc", sampled_rows_per_proc.data(), + sampled_rows_per_proc.size(), false); + h5_helper.putInteger("potential_rhs_nprocs", nprocs); + h5_helper.putInteger("potential_rhs_hr_idx_size", global_sampled_row.size()); + /* save right-hand side hyper-reduction operator */ h5_helper.putDoubleArray("potential_rhs_rom_inverse", pot_rhs_rom.getData(), num_pot_basis * num_pot_basis, false); @@ -272,6 +287,26 @@ void runPoissonROM(MGmolInterface *mgmol_) CAROM::Matrix pot_rom_inv(num_pot_basis, num_pot_basis, false); CAROM::Matrix pot_rhs_rom(num_pot_basis, num_pot_basis, false); CAROM::Vector pot_rhs_rescaler(num_pot_basis, false); + + int nprocs0 = -1, hr_idx_size = -1; + if (MPIdata::onpe0) + { + std::string rom_oper = rom_options.pot_rom_file; + CAROM::HDFDatabase h5_helper; + h5_helper.open(rom_oper, "r"); + + h5_helper.getInteger("potential_rhs_nprocs", nprocs0); + h5_helper.getInteger("potential_rhs_hr_idx_size", hr_idx_size); + if (nprocs0 != nprocs) + { + std::cerr << "runPoissonROM error: same number of processors must be used as for buildPoissonROM!\n" << std::endl; + MPI_Abort(MPI_COMM_WORLD, 0); + } + + h5_helper.close(); + } + mmpi.bcastGlobal(&hr_idx_size); + std::vector global_sampled_row(hr_idx_size), sampled_rows_per_proc(nprocs); /* Load ROM operator */ // read the file from PE0 only @@ -291,6 +326,12 @@ void runPoissonROM(MGmolInterface *mgmol_) h5_helper.getDoubleArray("potential_rom_inverse", pot_rom_inv.getData(), num_pot_basis * num_pot_basis, false); + /* load right-hand side hyper-reduction sample index */ + h5_helper.getIntegerArray("potential_rhs_sample_idx", global_sampled_row.data(), + global_sampled_row.size(), false); + h5_helper.getIntegerArray("potential_rhs_sampled_rows_per_proc", sampled_rows_per_proc.data(), + sampled_rows_per_proc.size(), false); + /* load right-hand side hyper-reduction operator */ h5_helper.getDoubleArray("potential_rhs_rom_inverse", pot_rhs_rom.getData(), num_pot_basis * num_pot_basis, false); @@ -301,6 +342,101 @@ void runPoissonROM(MGmolInterface *mgmol_) h5_helper.close(); } + // broadcast over all processes + mmpi.bcastGlobal(pot_rom.getData(), num_pot_basis * num_pot_basis, 0); + mmpi.bcastGlobal(pot_rom_inv.getData(), num_pot_basis * num_pot_basis, 0); + mmpi.bcastGlobal(pot_rhs_rom.getData(), num_pot_basis * num_pot_basis, 0); + mmpi.bcastGlobal(pot_rhs_rescaler.getData(), num_pot_basis, 0); + mmpi.bcastGlobal(global_sampled_row.data(), hr_idx_size, 0); + mmpi.bcastGlobal(sampled_rows_per_proc.data(), nprocs, 0); + + /* get local sampled row */ + std::vector offsets, sampled_row(sampled_rows_per_proc[rank]); + int num_global_sample = CAROM::get_global_offsets(sampled_rows_per_proc[rank], offsets); + for (int s = 0, gs = offsets[rank]; gs < offsets[rank+1]; gs++, s++) + sampled_row[s] = global_sampled_row[gs]; + + /* Load MGmol pointer and Potentials */ + MGmol *mgmol = static_cast *>(mgmol_); + Poisson *poisson = mgmol->electrostat_->getPoissonSolver(); + Potentials& pot = mgmol->getHamiltonian()->potential(); + std::shared_ptr> rho = mgmol->getRho(); + const int dim = pot.size(); + + /* ROM currently support only nspin=1 */ + CAROM_VERIFY(rho->rho_.size() == 1); + + /* load the test restart file */ + mgmol->loadRestartFile(rom_options.test_restart_file); + + /* get mgmol orbitals and copy to librom */ + OrbitalsType *orbitals = mgmol->getOrbitals(); + const int chrom_num = orbitals->chromatic_number(); + CAROM::Matrix psi(dim, chrom_num, true); + for (int c = 0; c < chrom_num; c++) + { + ORBDTYPE *d_psi = orbitals->getPsi(c); + for (int d = 0; d < dim ; d++) + psi.item(d, c) = *(d_psi + d); + } + + /* + get rom coefficients of orbitals based on orbital POD basis. + At this point, we take the FOM orbital as it is, + thus rom coefficients become the identity matrix. + */ + CAROM::Matrix rom_psi(chrom_num, chrom_num, false); + for (int i = 0; i < chrom_num; i++) + for (int j = 0; j < chrom_num; j++) + rom_psi(i, j) = (i == j) ? 1 : 0; + + /* get mgmol density matrix */ + ProjectedMatricesInterface *proj_matrices = orbitals->getProjMatrices(); + proj_matrices->updateSubMatX(); + SquareLocalMatrices& localX(proj_matrices->getLocalX()); + + bool dm_distributed = (localX.nmat() > 1); + CAROM_VERIFY(!dm_distributed); + + /* copy density matrix into librom */ + CAROM::Matrix dm(localX.getRawPtr(), localX.m(), localX.n(), dm_distributed, true); + + /* compute ROM rho using hyper reduction */ + CAROM::Vector sample_rho(1, true); // this will be resized in computeRhoOnSamplePts + computeRhoOnSamplePts(dm, psi, rom_psi, sampled_row, sample_rho); + sample_rho.gather(); + CAROM::Vector *rom_rho = pot_rhs_rom.mult(sample_rho); + + /* rescale the electron density */ + const double nel = ct.getNel(); + // volume element is already multiplied in pot_rhs_rescaler. + const double tcharge = rom_rho->inner_product(pot_rhs_rescaler); + *rom_rho *= nel / tcharge; + + /* compute ROM ion density using hyper reduction */ + std::shared_ptr ions = mgmol->getIons(); + CAROM::Vector sampled_rhoc(1, true); // this will be resized in evalIonDensityOnSamplePts + pot.evalIonDensityOnSamplePts(*ions, sampled_row, sampled_rhoc); + sampled_rhoc.gather(); + CAROM::Vector *rom_rhoc = pot_rhs_rom.mult(sampled_rhoc); + + /* rescale the ion density */ + const double ionic_charge = pot.getIonicCharge(); + // volume element is already multiplied in pot_rhs_rescaler. + const double comp_rho = rom_rhoc->inner_product(pot_rhs_rescaler); + *rom_rhoc *= ionic_charge / comp_rho; + + /* right-hand side */ + *rom_rho -= (*rom_rhoc); + *rom_rho *= (4.0 * M_PI); + + /* solve Poisson ROM */ + CAROM::Vector *rom_pot = pot_rom_inv.mult(*rom_rho); + + /* clean up */ + delete rom_rho; + delete rom_rhoc; + delete rom_pot; } /* test routines */ @@ -573,7 +709,7 @@ void testROMRhoOperator(MGmolInterface *mgmol_) MGmol *mgmol = static_cast *>(mgmol_); Poisson *poisson = mgmol->electrostat_->getPoissonSolver(); Potentials& pot = mgmol->getHamiltonian()->potential(); - std::shared_ptr> rho = NULL; // mgmol->getRho(); + std::shared_ptr> rho = mgmol->getRho(); const OrthoType ortho_type = rho->getOrthoType(); assert(ortho_type == OrthoType::Nonorthogonal); @@ -919,13 +1055,13 @@ void testROMIonDensity(MGmolInterface *mgmol_) CAROM_VERIFY(abs(fom_overlap_ions[test_idx][k][d] - ions->overlappingVL_ions()[k]->position(d)) < 1.0e-12); /* eval ion density on sample grid points */ - std::vector sampled_rhoc(sampled_row.size()); + CAROM::Vector sampled_rhoc(1, true); // this will be resized in evalIonDensityOnSamplePts pot.evalIonDensityOnSamplePts(*ions, sampled_row, sampled_rhoc); for (int d = 0; d < sampled_row.size(); d++) { - printf("rank %d, fom rhoc[%d]: %.3e, rom rhoc: %.3e\n", rank, sampled_row[d], fom_rhoc[test_idx][sampled_row[d]], sampled_rhoc[d]); - CAROM_VERIFY(abs(fom_rhoc[test_idx][sampled_row[d]] - sampled_rhoc[d]) < 1.0e-12); + printf("rank %d, fom rhoc[%d]: %.3e, rom rhoc: %.3e\n", rank, sampled_row[d], fom_rhoc[test_idx][sampled_row[d]], sampled_rhoc(d)); + CAROM_VERIFY(abs(fom_rhoc[test_idx][sampled_row[d]] - sampled_rhoc(d)) < 1.0e-12); } } From 4a00f5db9a9f01f9ef0d022ea244da36fc59e813 Mon Sep 17 00:00:00 2001 From: Seung Whan Chung Date: Wed, 30 Oct 2024 11:01:59 -0700 Subject: [PATCH 2/6] minor fix --- src/read_config.cc | 2 +- src/rom_workflows.cc | 4 ++-- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/src/read_config.cc b/src/read_config.cc index 02c0aa23..891c06c7 100644 --- a/src/read_config.cc +++ b/src/read_config.cc @@ -437,7 +437,7 @@ void setupROMConfigOption(po::options_description &rom_cfg) ("ROM.basis.number_of_potential_basis", po::value()->default_value(-1), "Number of potential POD basis to build Hartree potential ROM operator.") ("ROM.potential_rom_file", po::value()->default_value(""), - "File name to save/load potential ROM operators."); + "File name to save/load potential ROM operators.") ("ROM.online.test_restart_file", po::value()->default_value(""), "File name to test online ROM."); } diff --git a/src/rom_workflows.cc b/src/rom_workflows.cc index 49aea1d4..87bf97f7 100644 --- a/src/rom_workflows.cc +++ b/src/rom_workflows.cc @@ -327,9 +327,9 @@ void runPoissonROM(MGmolInterface *mgmol_) num_pot_basis * num_pot_basis, false); /* load right-hand side hyper-reduction sample index */ - h5_helper.getIntegerArray("potential_rhs_sample_idx", global_sampled_row.data(), + h5_helper.getIntegerArray("potential_rhs_hr_idx", global_sampled_row.data(), global_sampled_row.size(), false); - h5_helper.getIntegerArray("potential_rhs_sampled_rows_per_proc", sampled_rows_per_proc.data(), + h5_helper.getIntegerArray("potential_rhs_hr_idcs_per_proc", sampled_rows_per_proc.data(), sampled_rows_per_proc.size(), false); /* load right-hand side hyper-reduction operator */ From 26feee31fd1b7f7814e1b22a76cde5734a15711f Mon Sep 17 00:00:00 2001 From: Seung Whan Chung Date: Fri, 8 Nov 2024 14:39:29 -0800 Subject: [PATCH 3/6] wip with debugging. --- src/rom_Control.h | 1 + src/rom_workflows.cc | 143 ++++++++++++++++++++++++++++++++++++++++--- src/rom_workflows.h | 1 + 3 files changed, 135 insertions(+), 10 deletions(-) diff --git a/src/rom_Control.h b/src/rom_Control.h index 5653267f..990d3f1c 100644 --- a/src/rom_Control.h +++ b/src/rom_Control.h @@ -33,6 +33,7 @@ enum class ROMVariable { ORBITALS, POTENTIAL, + DENSITY, NONE }; diff --git a/src/rom_workflows.cc b/src/rom_workflows.cc index 87bf97f7..64c962c3 100644 --- a/src/rom_workflows.cc +++ b/src/rom_workflows.cc @@ -66,6 +66,9 @@ void readRestartFiles(MGmolInterface *mgmol_) Control& ct = *(Control::instance()); Mesh* mymesh = Mesh::instance(); const pb::PEenv& myPEenv = mymesh->peenv(); + MGmol_MPI& mmpi = *(MGmol_MPI::instance()); + const int rank = mmpi.mypeGlobal(); + const int nprocs = mmpi.size(); ROMPrivateOptions rom_options = ct.getROMOptions(); /* type of variable we intend to run POD */ @@ -131,6 +134,9 @@ void readRestartFiles(MGmolInterface *mgmol_) /* we save hartree potential */ basis_generator.takeSample(pot.vh_rho()); break; + + case ROMVariable::DENSITY: + basis_prefix += "_density"; } } basis_generator.writeSnapshot(); @@ -202,11 +208,12 @@ void buildROMPoissonOperator(MGmolInterface *mgmol_) delete col; } // for (int c = 0; c < num_pot_basis; c++) - /* DEIM hyperreduction */ + /* hyperreduction */ CAROM::Matrix pot_rhs_rom(num_pot_basis, num_pot_basis, false); std::vector global_sampled_row(num_pot_basis), sampled_rows_per_proc(nprocs); - DEIM(pot_basis, num_pot_basis, global_sampled_row, sampled_rows_per_proc, - pot_rhs_rom, rank, nprocs); + CAROM::Hyperreduction HR("deim"); + HR.ComputeSamples(pot_basis, num_pot_basis, global_sampled_row, sampled_rows_per_proc, + pot_rhs_rom, rank, nprocs, num_pot_basis); if (rank == 0) { int num_sample_rows = 0; @@ -276,6 +283,19 @@ void runPoissonROM(MGmolInterface *mgmol_) MPI_Abort(MPI_COMM_WORLD, 0); } + /* Load MGmol pointer and Potentials */ + MGmol *mgmol = static_cast *>(mgmol_); + Poisson *poisson = mgmol->electrostat_->getPoissonSolver(); + Potentials& pot = mgmol->getHamiltonian()->potential(); + std::shared_ptr> rho = mgmol->getRho(); + const int dim = pot.size(); + + /* GridFunc initialization inputs */ + const pb::Grid &grid(poisson->vh().grid()); + short bc[3]; + for (int d = 0; d < 3; d++) + bc[d] = poisson->vh().bc(d); + /* Load Hartree potential basis matrix */ std::string basis_file = rom_options.basis_file; const int num_pot_basis = rom_options.num_potbasis; @@ -355,13 +375,6 @@ void runPoissonROM(MGmolInterface *mgmol_) int num_global_sample = CAROM::get_global_offsets(sampled_rows_per_proc[rank], offsets); for (int s = 0, gs = offsets[rank]; gs < offsets[rank+1]; gs++, s++) sampled_row[s] = global_sampled_row[gs]; - - /* Load MGmol pointer and Potentials */ - MGmol *mgmol = static_cast *>(mgmol_); - Poisson *poisson = mgmol->electrostat_->getPoissonSolver(); - Potentials& pot = mgmol->getHamiltonian()->potential(); - std::shared_ptr> rho = mgmol->getRho(); - const int dim = pot.size(); /* ROM currently support only nspin=1 */ CAROM_VERIFY(rho->rho_.size() == 1); @@ -404,19 +417,70 @@ void runPoissonROM(MGmolInterface *mgmol_) /* compute ROM rho using hyper reduction */ CAROM::Vector sample_rho(1, true); // this will be resized in computeRhoOnSamplePts computeRhoOnSamplePts(dm, psi, rom_psi, sampled_row, sample_rho); + + /* check sampled rho */ + for (int s = 0; s < sampled_row.size(); s++) + { + printf("rank %d, rho[%d]: %.5e, sample_rho: %.5e\n", + rank, sampled_row[s], rho->rho_[0][sampled_row[s]], sample_rho(s)); + } + sample_rho.gather(); CAROM::Vector *rom_rho = pot_rhs_rom.mult(sample_rho); + if (rank == 0) + { + printf("rom rho before projection\n"); + for (int d = 0; d < num_pot_basis; d++) + printf("%.5e\t", rom_rho->item(d)); + printf("\n"); + + printf("pot_rhs_rescaler\n"); + for (int d = 0; d < num_pot_basis; d++) + printf("%.5e\t", pot_rhs_rescaler.item(d)); + printf("\n"); + } /* rescale the electron density */ const double nel = ct.getNel(); // volume element is already multiplied in pot_rhs_rescaler. const double tcharge = rom_rho->inner_product(pot_rhs_rescaler); *rom_rho *= nel / tcharge; + if (rank == 0) + printf("rank %d, rho scaler: %.3e\n", rank, nel / tcharge); + + /* check rho projection */ + CAROM::Vector fom_rho_vec(&rho->rho_[0][0], dim, true, false); + CAROM::Vector *rho_proj = pot_basis->transposeMult(fom_rho_vec); + CAROM::Vector *rho_proj2 = pot_basis->mult(*rho_proj); + (*rho_proj2) -= fom_rho_vec; + double rho_proj_error = rho_proj2->norm() / fom_rho_vec.norm(); + if (rank == 0) + { + printf("rom rho\n"); + for (int d = 0; d < num_pot_basis; d++) + printf("%.5e\t", rom_rho->item(d)); + printf("\n"); + printf("fom rho projection\n"); + for (int d = 0; d < num_pot_basis; d++) + printf("%.5e\t", rho_proj->item(d)); + printf("\n"); + + printf("rho proj error: %.5e\n", rho_proj_error); + } /* compute ROM ion density using hyper reduction */ std::shared_ptr ions = mgmol->getIons(); CAROM::Vector sampled_rhoc(1, true); // this will be resized in evalIonDensityOnSamplePts pot.evalIonDensityOnSamplePts(*ions, sampled_row, sampled_rhoc); + + /* check sampled rhoc */ + RHODTYPE *rho_comp = pot.rho_comp(); + for (int s = 0; s < sampled_row.size(); s++) + { + printf("rank %d, rhoc[%d]: %.5e, sampled_rhoc: %.5e\n", + rank, sampled_row[s], rho_comp[sampled_row[s]], sampled_rhoc(s)); + } + sampled_rhoc.gather(); CAROM::Vector *rom_rhoc = pot_rhs_rom.mult(sampled_rhoc); @@ -425,6 +489,8 @@ void runPoissonROM(MGmolInterface *mgmol_) // volume element is already multiplied in pot_rhs_rescaler. const double comp_rho = rom_rhoc->inner_product(pot_rhs_rescaler); *rom_rhoc *= ionic_charge / comp_rho; + if (rank == 0) + printf("rank %d, rhoc scaler: %.3e\n", rank, ionic_charge / comp_rho); /* right-hand side */ *rom_rho -= (*rom_rhoc); @@ -433,10 +499,48 @@ void runPoissonROM(MGmolInterface *mgmol_) /* solve Poisson ROM */ CAROM::Vector *rom_pot = pot_rom_inv.mult(*rom_rho); + /* data array to lift up rom solution */ + std::vector test_sol(dim); + /* get librom view-vector of test_sol[s] */ + CAROM::Vector test_sol_vec(test_sol.data(), dim, true, false); + pot_basis->mult(*rom_pot, test_sol_vec); + + /* mgmol grid function for lifted-up fom solution */ + pb::GridFunc testsol_gf(grid, bc[0], bc[1], bc[2]); + pb::GridFunc fomsol_gf(grid, bc[0], bc[1], bc[2]); + testsol_gf.assign(test_sol.data(), 'd'); + // POTDTYPE *vh_rho = pot.vh_rho(); + // for (int d = 0; d < dim; d++) + // (vh_rho[d]) += 1.0e-3; + fomsol_gf.assign(pot.vh_rho(), 'd'); + + testsol_gf -= fomsol_gf; + double rel_error = testsol_gf.norm2() / fomsol_gf.norm2(); + if (rank == 0) + printf("relative error: %.3e\n", rel_error); + + /* librom view vector for fom solution */ + CAROM::Vector fom_sol_vec(pot.vh_rho(), dim, true, false); + CAROM::Vector *fom_proj = pot_basis->transposeMult(fom_sol_vec); + if (rank == 0) + { + printf("rom vector\n"); + for (int d = 0; d < num_pot_basis; d++) + printf("%.5e\t", rom_pot->item(d)); + printf("\n"); + printf("fom projection\n"); + for (int d = 0; d < num_pot_basis; d++) + printf("%.5e\t", fom_proj->item(d)); + printf("\n"); + } + /* clean up */ delete rom_rho; delete rom_rhoc; delete rom_pot; + delete fom_proj; + delete rho_proj; + delete rho_proj2; } /* test routines */ @@ -713,6 +817,12 @@ void testROMRhoOperator(MGmolInterface *mgmol_) const OrthoType ortho_type = rho->getOrthoType(); assert(ortho_type == OrthoType::Nonorthogonal); + /* GridFunc initialization inputs */ + const pb::Grid &grid(poisson->vh().grid()); + short bc[3]; + for (int d = 0; d < 3; d++) + bc[d] = poisson->vh().bc(d); + /* potential should have the same size as rho */ const int dim = pot.size(); @@ -734,11 +844,24 @@ void testROMRhoOperator(MGmolInterface *mgmol_) /* Collect the restart files */ std::string filename; + pb::GridFunc rhogf(grid, bc[0], bc[1], bc[2]); + std::vector rho_outvec(rho->rho_[0].size()); for (int k = minidx; k <= maxidx; k++) { filename = string_format(rom_options.restart_file_fmt, k); mgmol->loadRestartFile(filename); basis_generator.takeSample(&rho->rho_[0][0]); + + rhogf.assign(&rho->rho_[0][0]); + rhogf.init_vect(rho_outvec.data(), 'd'); + for (int d = 0; d < rho_outvec.size(); d++) + { + const double error = abs(rho->rho_[0][d] - rho_outvec[d]); + if (error > 0.0) + printf("rank %d, rho[%d]: %.15e, sample_rho: %.15e\n", + rank, d, rho->rho_[0][d], rho_outvec[d]); + CAROM_VERIFY(error == 0.0); + } } // basis_generator.writeSnapshot(); const CAROM::Matrix rho_snapshots(*basis_generator.getSnapshotMatrix()); diff --git a/src/rom_workflows.h b/src/rom_workflows.h index 9222249f..313f6f59 100644 --- a/src/rom_workflows.h +++ b/src/rom_workflows.h @@ -36,6 +36,7 @@ namespace po = boost::program_options; #include "librom.h" +#include "hyperreduction/Hyperreduction.h" #include "utils/HDFDatabase.h" #include "utils/mpi_utils.h" From 3e31f71df827a834f524a9b88887c540b7d7b030 Mon Sep 17 00:00:00 2001 From: Seung Whan Chung Date: Tue, 7 Jan 2025 12:21:01 -0800 Subject: [PATCH 4/6] rho hyperreduction test passed. --- src/rom_main.cc | 2 +- src/rom_workflows.cc | 18 ++++++++++++------ 2 files changed, 13 insertions(+), 7 deletions(-) diff --git a/src/rom_main.cc b/src/rom_main.cc index 7eb323a7..6e56dd26 100644 --- a/src/rom_main.cc +++ b/src/rom_main.cc @@ -137,13 +137,13 @@ int main(int argc, char** argv) testROMRhoOperator(mgmol); else testROMRhoOperator(mgmol); + break; case (ROMStage::TEST_ION): if (ct.isLocMode()) testROMIonDensity(mgmol); else testROMIonDensity(mgmol); - break; default: diff --git a/src/rom_workflows.cc b/src/rom_workflows.cc index 64c962c3..6380d473 100644 --- a/src/rom_workflows.cc +++ b/src/rom_workflows.cc @@ -802,6 +802,9 @@ void testROMRhoOperator(MGmolInterface *mgmol_) const int rank = mmpi.mypeGlobal(); const int nprocs = mmpi.size(); + static std::random_device rd; // Will be used to obtain a seed for the random number engine + static std::mt19937 gen(rd()); // Standard mersenne_twister_engine seeded with rd(){} + // if (ct.isLocMode()) // printf("LocMode is On!\n"); // else @@ -890,7 +893,10 @@ void testROMRhoOperator(MGmolInterface *mgmol_) sampled_row[s] = global_sampled_row[gs]; /* load only the first restart file for now */ - const int test_idx = 2; + std::uniform_int_distribution<> distrib(0, num_restart-1); + int test_idx = distrib(gen); + mmpi.bcastGlobal(&test_idx); + if (rank == 0) printf("test index: %d\n", test_idx); filename = string_format(rom_options.restart_file_fmt, test_idx + minidx); /* @@ -967,10 +973,10 @@ void testROMRhoOperator(MGmolInterface *mgmol_) for (int s = 0; s < sampled_row.size(); s++) { const double error = abs(rho->rho_[0][sampled_row[s]] - sample_rho(s)); - if (error > 1.0e-4) + if (error > 1.0e-10) printf("rank %d, rho[%d]: %.5e, sample_rho: %.5e, librom_snapshot: %.5e\n", rank, sampled_row[s], rho->rho_[0][sampled_row[s]], sample_rho(s), rho_snapshots(sampled_row[s], test_idx)); - CAROM_VERIFY(error < 1.0e-4); + CAROM_VERIFY(error < 1.0e-10); } sample_rho.gather(); @@ -978,16 +984,16 @@ void testROMRhoOperator(MGmolInterface *mgmol_) CAROM::Vector *rom_rho = rho_basis_inv.mult(sample_rho); for (int d = 0; d < rom_rho->dim(); d++) { - if ((rank == 0) && (abs(proj_rho->item(d, test_idx) - rom_rho->item(d)) > 1.0e-3)) + if ((rank == 0)) printf("rom_rho error: %.3e\n", abs(proj_rho->item(d, test_idx) - rom_rho->item(d))); - CAROM_VERIFY(abs(proj_rho->item(d, test_idx) - rom_rho->item(d)) < 1.0e-3); + CAROM_VERIFY(abs(proj_rho->item(d, test_idx) - rom_rho->item(d)) < 1.0e-10); } CAROM::Vector *fom_rho = rho_basis->mult(*rom_rho); CAROM_VERIFY(fom_rho->dim() == rho->rho_[0].size()); for (int d = 0; d < fom_rho->dim(); d++) - CAROM_VERIFY(abs(fom_rho->item(d) - rho->rho_[0][d]) < 1.0e-4); + CAROM_VERIFY(abs(fom_rho->item(d) - rho->rho_[0][d]) < 1.0e-10); delete rom_rho; delete fom_rho; From de75383255d6a00a4ee4fbe848604c48c41c7ce4 Mon Sep 17 00:00:00 2001 From: Seung Whan Chung Date: Tue, 7 Jan 2025 17:48:43 -0800 Subject: [PATCH 5/6] testPoisson includes hyperreduction. hyperreduction must be done with right-handside pod basis. --- src/rom_workflows.cc | 214 +++++++++++++++++++++++++++++-------------- 1 file changed, 143 insertions(+), 71 deletions(-) diff --git a/src/rom_workflows.cc b/src/rom_workflows.cc index 6380d473..7112860e 100644 --- a/src/rom_workflows.cc +++ b/src/rom_workflows.cc @@ -130,11 +130,20 @@ void readRestartFiles(MGmolInterface *mgmol_) break; case ROMVariable::POTENTIAL: + { basis_prefix += "_potential"; + /* + potential in restart file is not consistent with + density/density matrix in the same file. + here we recompute potential. + */ + std::shared_ptr ions = mgmol->getIons(); + mgmol->update_pot(*ions); + /* we save hartree potential */ basis_generator.takeSample(pot.vh_rho()); break; - + } case ROMVariable::DENSITY: basis_prefix += "_density"; } @@ -164,6 +173,7 @@ void buildROMPoissonOperator(MGmolInterface *mgmol_) /* Load Hartree potential basis matrix */ std::string basis_file = rom_options.basis_file; + basis_file += "_potential"; const int num_pot_basis = rom_options.num_potbasis; CAROM::BasisReader basis_reader(basis_file); CAROM::Matrix *pot_basis = basis_reader.getSpatialBasis(num_pot_basis); @@ -298,6 +308,7 @@ void runPoissonROM(MGmolInterface *mgmol_) /* Load Hartree potential basis matrix */ std::string basis_file = rom_options.basis_file; + basis_file += "_potential"; const int num_pot_basis = rom_options.num_potbasis; CAROM::BasisReader basis_reader(basis_file); CAROM::Matrix *pot_basis = basis_reader.getSpatialBasis(num_pot_basis); @@ -421,52 +432,52 @@ void runPoissonROM(MGmolInterface *mgmol_) /* check sampled rho */ for (int s = 0; s < sampled_row.size(); s++) { - printf("rank %d, rho[%d]: %.5e, sample_rho: %.5e\n", - rank, sampled_row[s], rho->rho_[0][sampled_row[s]], sample_rho(s)); + printf("rank %d, rho[%d]: %.5e, sample_rho: %.5e, diff: %.5e\n", + rank, sampled_row[s], rho->rho_[0][sampled_row[s]], sample_rho(s), abs(rho->rho_[0][sampled_row[s]] - sample_rho(s))); } - sample_rho.gather(); - CAROM::Vector *rom_rho = pot_rhs_rom.mult(sample_rho); - if (rank == 0) - { - printf("rom rho before projection\n"); - for (int d = 0; d < num_pot_basis; d++) - printf("%.5e\t", rom_rho->item(d)); - printf("\n"); - - printf("pot_rhs_rescaler\n"); - for (int d = 0; d < num_pot_basis; d++) - printf("%.5e\t", pot_rhs_rescaler.item(d)); - printf("\n"); - } - - /* rescale the electron density */ - const double nel = ct.getNel(); - // volume element is already multiplied in pot_rhs_rescaler. - const double tcharge = rom_rho->inner_product(pot_rhs_rescaler); - *rom_rho *= nel / tcharge; - if (rank == 0) - printf("rank %d, rho scaler: %.3e\n", rank, nel / tcharge); - - /* check rho projection */ - CAROM::Vector fom_rho_vec(&rho->rho_[0][0], dim, true, false); - CAROM::Vector *rho_proj = pot_basis->transposeMult(fom_rho_vec); - CAROM::Vector *rho_proj2 = pot_basis->mult(*rho_proj); - (*rho_proj2) -= fom_rho_vec; - double rho_proj_error = rho_proj2->norm() / fom_rho_vec.norm(); - if (rank == 0) - { - printf("rom rho\n"); - for (int d = 0; d < num_pot_basis; d++) - printf("%.5e\t", rom_rho->item(d)); - printf("\n"); - printf("fom rho projection\n"); - for (int d = 0; d < num_pot_basis; d++) - printf("%.5e\t", rho_proj->item(d)); - printf("\n"); + // sample_rho.gather(); + // CAROM::Vector *rom_rho = pot_rhs_rom.mult(sample_rho); + // if (rank == 0) + // { + // printf("rom rho before projection\n"); + // for (int d = 0; d < num_pot_basis; d++) + // printf("%.5e\t", rom_rho->item(d)); + // printf("\n"); + + // printf("pot_rhs_rescaler\n"); + // for (int d = 0; d < num_pot_basis; d++) + // printf("%.5e\t", pot_rhs_rescaler.item(d)); + // printf("\n"); + // } - printf("rho proj error: %.5e\n", rho_proj_error); - } + // /* rescale the electron density */ + // const double nel = ct.getNel(); + // // volume element is already multiplied in pot_rhs_rescaler. + // const double tcharge = rom_rho->inner_product(pot_rhs_rescaler); + // *rom_rho *= nel / tcharge; + // if (rank == 0) + // printf("rank %d, rho scaler: %.3e\n", rank, nel / tcharge); + + // /* check rho projection */ + // CAROM::Vector fom_rho_vec(&rho->rho_[0][0], dim, true, false); + // CAROM::Vector *rho_proj = pot_basis->transposeMult(fom_rho_vec); + // CAROM::Vector *rho_proj2 = pot_basis->mult(*rho_proj); + // (*rho_proj2) -= fom_rho_vec; + // double rho_proj_error = rho_proj2->norm() / fom_rho_vec.norm(); + // if (rank == 0) + // { + // printf("rom rho\n"); + // for (int d = 0; d < num_pot_basis; d++) + // printf("%.5e\t", rom_rho->item(d)); + // printf("\n"); + // printf("fom rho projection\n"); + // for (int d = 0; d < num_pot_basis; d++) + // printf("%.5e\t", rho_proj->item(d)); + // printf("\n"); + + // printf("rho proj error: %.5e\n", rho_proj_error); + // } /* compute ROM ion density using hyper reduction */ std::shared_ptr ions = mgmol->getIons(); @@ -477,27 +488,49 @@ void runPoissonROM(MGmolInterface *mgmol_) RHODTYPE *rho_comp = pot.rho_comp(); for (int s = 0; s < sampled_row.size(); s++) { - printf("rank %d, rhoc[%d]: %.5e, sampled_rhoc: %.5e\n", - rank, sampled_row[s], rho_comp[sampled_row[s]], sampled_rhoc(s)); + printf("rank %d, rhoc[%d]: %.5e, sampled_rhoc: %.5e, diff: %.5e\n", + rank, sampled_row[s], rho_comp[sampled_row[s]], sampled_rhoc(s), abs(rho_comp[sampled_row[s]] - sampled_rhoc(s))); } - sampled_rhoc.gather(); - CAROM::Vector *rom_rhoc = pot_rhs_rom.mult(sampled_rhoc); + // sampled_rhoc.gather(); + // CAROM::Vector *rom_rhoc = pot_rhs_rom.mult(sampled_rhoc); - /* rescale the ion density */ - const double ionic_charge = pot.getIonicCharge(); - // volume element is already multiplied in pot_rhs_rescaler. - const double comp_rho = rom_rhoc->inner_product(pot_rhs_rescaler); - *rom_rhoc *= ionic_charge / comp_rho; - if (rank == 0) - printf("rank %d, rhoc scaler: %.3e\n", rank, ionic_charge / comp_rho); + // /* rescale the ion density */ + // const double ionic_charge = pot.getIonicCharge(); + // // volume element is already multiplied in pot_rhs_rescaler. + // const double comp_rho = rom_rhoc->inner_product(pot_rhs_rescaler); + // *rom_rhoc *= ionic_charge / comp_rho; + // if (rank == 0) + // printf("rank %d, rhoc scaler: %.3e\n", rank, ionic_charge / comp_rho); /* right-hand side */ - *rom_rho -= (*rom_rhoc); - *rom_rho *= (4.0 * M_PI); + CAROM::Vector sampled_rhs(sample_rho); + sampled_rhs -= sampled_rhoc; + sampled_rhs *= (4.0 * M_PI); + + /* get poisson right-hand side by applying Laplace operator to potential */ + pb::GridFunc fomsol_gf(grid, bc[0], bc[1], bc[2]); + fomsol_gf.assign(pot.vh_rho(), 'd'); + /* apply Laplace operator */ + pb::GridFunc fomrhs_gf(grid, bc[0], bc[1], bc[2]); + poisson->applyOperator(fomsol_gf, fomrhs_gf); + CAROM::Vector fomrhs(pot_basis->numRows(), true); + fomrhs_gf.init_vect(fomrhs.getData(), 'd'); + /* check sampled rhs */ + for (int s = 0; s < sampled_row.size(); s++) + { + printf("rank %d, rhs[%d]: %.5e, sampled_rhs: %.5e, diff: %.5e\n", + rank, sampled_row[s], fomrhs(sampled_row[s]), sampled_rhs(s), abs(fomrhs(sampled_row[s]) - sampled_rhs(s))); + } + + /* ROM right-hand side */ + sampled_rhs.gather(); + CAROM::Vector *rom_rhs = pot_rhs_rom.mult(sampled_rhs); + // *rom_rho -= (*rom_rhoc); + // *rom_rho *= (4.0 * M_PI); /* solve Poisson ROM */ - CAROM::Vector *rom_pot = pot_rom_inv.mult(*rom_rho); + CAROM::Vector *rom_pot = pot_rom_inv.mult(*rom_rhs); /* data array to lift up rom solution */ std::vector test_sol(dim); @@ -507,12 +540,7 @@ void runPoissonROM(MGmolInterface *mgmol_) /* mgmol grid function for lifted-up fom solution */ pb::GridFunc testsol_gf(grid, bc[0], bc[1], bc[2]); - pb::GridFunc fomsol_gf(grid, bc[0], bc[1], bc[2]); testsol_gf.assign(test_sol.data(), 'd'); - // POTDTYPE *vh_rho = pot.vh_rho(); - // for (int d = 0; d < dim; d++) - // (vh_rho[d]) += 1.0e-3; - fomsol_gf.assign(pot.vh_rho(), 'd'); testsol_gf -= fomsol_gf; double rel_error = testsol_gf.norm2() / fomsol_gf.norm2(); @@ -532,15 +560,20 @@ void runPoissonROM(MGmolInterface *mgmol_) for (int d = 0; d < num_pot_basis; d++) printf("%.5e\t", fom_proj->item(d)); printf("\n"); + printf("ratio\n"); + for (int d = 0; d < num_pot_basis; d++) + printf("%.5e\t", rom_pot->item(d) / fom_proj->item(d)); + printf("\n"); } /* clean up */ - delete rom_rho; - delete rom_rhoc; - delete rom_pot; - delete fom_proj; - delete rho_proj; - delete rho_proj2; + delete rom_rhs; + // delete rom_rho; + // delete rom_rhoc; + // delete rom_pot; + // delete fom_proj; + // delete rho_proj; + // delete rho_proj2; } /* test routines */ @@ -551,6 +584,9 @@ void testROMPoissonOperator(MGmolInterface *mgmol_) Control& ct = *(Control::instance()); Mesh* mymesh = Mesh::instance(); const pb::PEenv& myPEenv = mymesh->peenv(); + MGmol_MPI& mmpi = *(MGmol_MPI::instance()); + const int rank = mmpi.mypeGlobal(); + const int nprocs = mmpi.size(); ROMPrivateOptions rom_options = ct.getROMOptions(); @@ -611,14 +647,21 @@ void testROMPoissonOperator(MGmolInterface *mgmol_) std::string basis_prefix = "test_poisson"; CAROM::Options svd_options(dim, nsnapshot, 1); CAROM::BasisGenerator basis_generator(svd_options, false, basis_prefix); + CAROM::BasisGenerator rhs_basis_generator(svd_options, false, "test_rhs"); /* Collect snapshots and train POD basis */ for (int s = 0; s < nsnapshot; s++) + { basis_generator.takeSample(fom_sol[s].data()); + rhs_basis_generator.takeSample(rhs[s].data()); + } basis_generator.endSamples(); + rhs_basis_generator.endSamples(); /* Load POD basis. We use the maximum number of basis vectors. */ const CAROM::Matrix *pot_basis = basis_generator.getSpatialBasis(); + const CAROM::Matrix *rhs_basis = rhs_basis_generator.getSpatialBasis(); + CAROM::Matrix *rhs2pot = pot_basis->transposeMult(rhs_basis); /* Check if full projection preserves FOM solution */ for (int c = 0; c < nsnapshot; c++) @@ -737,6 +780,27 @@ void testROMPoissonOperator(MGmolInterface *mgmol_) } delete identity; + /* set up hyper-reduction for right-hand side. */ + CAROM::Matrix rhs_basis_inv(nsnapshot, nsnapshot, false); + std::vector global_sampled_row(nsnapshot), sampled_rows_per_proc(nprocs); + CAROM::Hyperreduction HR("deim"); + HR.ComputeSamples(rhs_basis, nsnapshot, global_sampled_row, sampled_rows_per_proc, + rhs_basis_inv, rank, nprocs, nsnapshot); + CAROM::Matrix *pot_rhs_rom = rhs2pot->mult(rhs_basis_inv); + if (rank == 0) + { + int num_sample_rows = 0; + for (int k = 0; k < sampled_rows_per_proc.size(); k++) + num_sample_rows += sampled_rows_per_proc[k]; + printf("number of sampled row: %d\n", num_sample_rows); + } + + /* get local sampled row */ + std::vector offsets, sampled_row(sampled_rows_per_proc[rank]); + int num_global_sample = CAROM::get_global_offsets(sampled_rows_per_proc[rank], offsets); + for (int s = 0, gs = offsets[rank]; gs < offsets[rank+1]; gs++, s++) + sampled_row[s] = global_sampled_row[gs]; + /* Test with sample RHS. ROM must be able to 100% reproduce the FOM solution. */ std::vector rom_sol(0), rom_rhs(0); std::vector> test_sol(nsnapshot); @@ -745,8 +809,16 @@ void testROMPoissonOperator(MGmolInterface *mgmol_) /* get librom view-vector of rhs[s] */ CAROM::Vector fom_rhs(rhs[s].data(), dim, true, false); - /* project onto POD basis */ - rom_rhs.push_back(pot_basis->transposeMult(fom_rhs)); + // /* project onto POD basis */ + // rom_rhs.push_back(pot_basis->transposeMult(fom_rhs)); + /* get sampled rhs */ + CAROM::Vector sample_rhs(1, true); + sample_rhs.setSize(sampled_row.size()); + for (int s = 0; s < sampled_row.size(); s++) + sample_rhs(s) = fom_rhs(sampled_row[s]); + /* get ROM rhs projection */ + sample_rhs.gather(); + rom_rhs.push_back(pot_rhs_rom->mult(sample_rhs)); /* FOM FD operator scales rhs by 4pi */ *rom_rhs.back() *= 4. * M_PI; From 3fd6a236f82efb8880d44799866ecd93d546bd92 Mon Sep 17 00:00:00 2001 From: Seung Whan Chung Date: Wed, 8 Jan 2025 14:24:05 -0800 Subject: [PATCH 6/6] hyperreduction using rho basis. --- src/Control.cc | 6 ++ src/read_config.cc | 2 + src/rom_Control.h | 1 + src/rom_workflows.cc | 171 ++++++++++++++++++++++++------------------- 4 files changed, 105 insertions(+), 75 deletions(-) diff --git a/src/Control.cc b/src/Control.cc index cae4f08b..4f2fb2d1 100644 --- a/src/Control.cc +++ b/src/Control.cc @@ -2055,6 +2055,8 @@ void Control::setROMOptions(const boost::program_options::variables_map& vm) rom_pri_option.variable = ROMVariable::ORBITALS; else if (str.compare("potential") == 0) rom_pri_option.variable = ROMVariable::POTENTIAL; + else if (str.compare("density") == 0) + rom_pri_option.variable = ROMVariable::DENSITY; else rom_pri_option.variable = ROMVariable::NONE; @@ -2064,6 +2066,7 @@ void Control::setROMOptions(const boost::program_options::variables_map& vm) rom_pri_option.compare_md = vm["ROM.basis.compare_md"].as(); rom_pri_option.num_orbbasis = vm["ROM.basis.number_of_orbital_basis"].as(); rom_pri_option.num_potbasis = vm["ROM.basis.number_of_potential_basis"].as(); + rom_pri_option.num_rhobasis = vm["ROM.basis.number_of_density_basis"].as(); rom_pri_option.pot_rom_file = vm["ROM.potential_rom_file"].as(); rom_pri_option.test_restart_file = vm["ROM.online.test_restart_file"].as(); @@ -2125,4 +2128,7 @@ void Control::syncROMOptions() mpirc = MPI_Bcast(&rom_pri_option.num_potbasis, 1, MPI_INT, 0, comm_global_); bcast_check(mpirc); + + mpirc = MPI_Bcast(&rom_pri_option.num_rhobasis, 1, MPI_INT, 0, comm_global_); + bcast_check(mpirc); } diff --git a/src/read_config.cc b/src/read_config.cc index 891c06c7..1d15e62e 100644 --- a/src/read_config.cc +++ b/src/read_config.cc @@ -436,6 +436,8 @@ void setupROMConfigOption(po::options_description &rom_cfg) "Number of orbital POD basis.") ("ROM.basis.number_of_potential_basis", po::value()->default_value(-1), "Number of potential POD basis to build Hartree potential ROM operator.") + ("ROM.basis.number_of_density_basis", po::value()->default_value(-1), + "Number of density POD basis to build Hartree potential ROM operator.") ("ROM.potential_rom_file", po::value()->default_value(""), "File name to save/load potential ROM operators.") ("ROM.online.test_restart_file", po::value()->default_value(""), diff --git a/src/rom_Control.h b/src/rom_Control.h index 990d3f1c..819b2b53 100644 --- a/src/rom_Control.h +++ b/src/rom_Control.h @@ -56,6 +56,7 @@ struct ROMPrivateOptions bool compare_md = false; int num_orbbasis = -1; int num_potbasis = -1; + int num_rhobasis = -1; std::string pot_rom_file = ""; /* options for online Poisson ROM */ diff --git a/src/rom_workflows.cc b/src/rom_workflows.cc index 7112860e..1630781b 100644 --- a/src/rom_workflows.cc +++ b/src/rom_workflows.cc @@ -84,6 +84,7 @@ void readRestartFiles(MGmolInterface *mgmol_) MGmol *mgmol = static_cast *>(mgmol_); OrbitalsType *orbitals = mgmol->getOrbitals(); Potentials& pot = mgmol->getHamiltonian()->potential(); + std::shared_ptr> rho = mgmol->getRho(); std::string filename; /* Determine basis prefix, dimension, and sample size */ @@ -104,6 +105,11 @@ void readRestartFiles(MGmolInterface *mgmol_) basis_prefix += "_potential"; dim = pot.size(); break; + + case ROMVariable::DENSITY: // electronic density + basis_prefix += "_density"; + dim = pot.size(); // same as potential + break; default: ct.global_exit(-1); @@ -131,7 +137,6 @@ void readRestartFiles(MGmolInterface *mgmol_) case ROMVariable::POTENTIAL: { - basis_prefix += "_potential"; /* potential in restart file is not consistent with density/density matrix in the same file. @@ -145,7 +150,8 @@ void readRestartFiles(MGmolInterface *mgmol_) break; } case ROMVariable::DENSITY: - basis_prefix += "_density"; + basis_generator.takeSample(&rho->rho_[0][0]); + break; } } basis_generator.writeSnapshot(); @@ -173,10 +179,12 @@ void buildROMPoissonOperator(MGmolInterface *mgmol_) /* Load Hartree potential basis matrix */ std::string basis_file = rom_options.basis_file; - basis_file += "_potential"; const int num_pot_basis = rom_options.num_potbasis; - CAROM::BasisReader basis_reader(basis_file); - CAROM::Matrix *pot_basis = basis_reader.getSpatialBasis(num_pot_basis); + CAROM::BasisReader pot_basis_reader(basis_file + "_potential"); + CAROM::Matrix *pot_basis = pot_basis_reader.getSpatialBasis(num_pot_basis); + const int num_rho_basis = rom_options.num_rhobasis; + CAROM::BasisReader rho_basis_reader(basis_file + "_density"); + CAROM::Matrix *rho_basis = rho_basis_reader.getSpatialBasis(num_rho_basis); /* Load PoissonSolver pointer */ MGmol *mgmol = static_cast *>(mgmol_); @@ -219,11 +227,12 @@ void buildROMPoissonOperator(MGmolInterface *mgmol_) } // for (int c = 0; c < num_pot_basis; c++) /* hyperreduction */ - CAROM::Matrix pot_rhs_rom(num_pot_basis, num_pot_basis, false); - std::vector global_sampled_row(num_pot_basis), sampled_rows_per_proc(nprocs); + CAROM::Matrix rho_basis_inv(num_rho_basis, num_rho_basis, false); + + std::vector global_sampled_row(num_rho_basis), sampled_rows_per_proc(nprocs); CAROM::Hyperreduction HR("deim"); - HR.ComputeSamples(pot_basis, num_pot_basis, global_sampled_row, sampled_rows_per_proc, - pot_rhs_rom, rank, nprocs, num_pot_basis); + HR.ComputeSamples(rho_basis, num_rho_basis, global_sampled_row, sampled_rows_per_proc, + rho_basis_inv, rank, nprocs, num_rho_basis); if (rank == 0) { int num_sample_rows = 0; @@ -232,11 +241,17 @@ void buildROMPoissonOperator(MGmolInterface *mgmol_) printf("number of sampled row: %d\n", num_sample_rows); } + /* projection of rho onto potential pod basis */ + CAROM::Matrix pot_rho_rom(num_pot_basis, num_rho_basis, false); + CAROM::Matrix tmp(pot_rho_rom); + pot_basis->transposeMult(*rho_basis, tmp); + tmp.mult(rho_basis_inv, pot_rho_rom); + /* ROM rescaleTotalCharge operator */ - CAROM::Vector fom_ones(pot_basis->numRows(), true); - CAROM::Vector rom_ones(num_pot_basis, false); + CAROM::Vector fom_ones(rho_basis->numRows(), true); + CAROM::Vector rom_ones(num_rho_basis, false); fom_ones = mymesh->grid().vel(); // volume element - pot_basis->transposeMult(fom_ones, rom_ones); + rho_basis->transposeMult(fom_ones, rom_ones); /* Save ROM operator */ // write the file from PE0 only @@ -263,12 +278,12 @@ void buildROMPoissonOperator(MGmolInterface *mgmol_) h5_helper.putInteger("potential_rhs_hr_idx_size", global_sampled_row.size()); /* save right-hand side hyper-reduction operator */ - h5_helper.putDoubleArray("potential_rhs_rom_inverse", pot_rhs_rom.getData(), - num_pot_basis * num_pot_basis, false); + h5_helper.putDoubleArray("potential_rho_rom_inverse", pot_rho_rom.getData(), + num_pot_basis * num_rho_basis, false); /* save right-hand side rescaling operator */ - h5_helper.putDoubleArray("potential_rhs_rescaler", rom_ones.getData(), - num_pot_basis, false); + h5_helper.putDoubleArray("potential_rho_rescaler", rom_ones.getData(), + num_rho_basis, false); h5_helper.close(); } @@ -310,14 +325,15 @@ void runPoissonROM(MGmolInterface *mgmol_) std::string basis_file = rom_options.basis_file; basis_file += "_potential"; const int num_pot_basis = rom_options.num_potbasis; + const int num_rho_basis = rom_options.num_rhobasis; CAROM::BasisReader basis_reader(basis_file); CAROM::Matrix *pot_basis = basis_reader.getSpatialBasis(num_pot_basis); /* initialize rom operator variables */ CAROM::Matrix pot_rom(num_pot_basis, num_pot_basis, false); CAROM::Matrix pot_rom_inv(num_pot_basis, num_pot_basis, false); - CAROM::Matrix pot_rhs_rom(num_pot_basis, num_pot_basis, false); - CAROM::Vector pot_rhs_rescaler(num_pot_basis, false); + CAROM::Matrix pot_rho_rom(num_pot_basis, num_rho_basis, false); + CAROM::Vector pot_rho_rescaler(num_rho_basis, false); int nprocs0 = -1, hr_idx_size = -1; if (MPIdata::onpe0) @@ -364,11 +380,11 @@ void runPoissonROM(MGmolInterface *mgmol_) sampled_rows_per_proc.size(), false); /* load right-hand side hyper-reduction operator */ - h5_helper.getDoubleArray("potential_rhs_rom_inverse", pot_rhs_rom.getData(), - num_pot_basis * num_pot_basis, false); + h5_helper.getDoubleArray("potential_rho_rom_inverse", pot_rho_rom.getData(), + num_pot_basis * num_rho_basis, false); /* load right-hand side rescaling operator */ - h5_helper.getDoubleArray("potential_rhs_rescaler", pot_rhs_rescaler.getData(), + h5_helper.getDoubleArray("potential_rho_rescaler", pot_rho_rescaler.getData(), num_pot_basis, false); h5_helper.close(); @@ -376,8 +392,8 @@ void runPoissonROM(MGmolInterface *mgmol_) // broadcast over all processes mmpi.bcastGlobal(pot_rom.getData(), num_pot_basis * num_pot_basis, 0); mmpi.bcastGlobal(pot_rom_inv.getData(), num_pot_basis * num_pot_basis, 0); - mmpi.bcastGlobal(pot_rhs_rom.getData(), num_pot_basis * num_pot_basis, 0); - mmpi.bcastGlobal(pot_rhs_rescaler.getData(), num_pot_basis, 0); + mmpi.bcastGlobal(pot_rho_rom.getData(), num_pot_basis * num_rho_basis, 0); + mmpi.bcastGlobal(pot_rho_rescaler.getData(), num_rho_basis, 0); mmpi.bcastGlobal(global_sampled_row.data(), hr_idx_size, 0); mmpi.bcastGlobal(sampled_rows_per_proc.data(), nprocs, 0); @@ -436,8 +452,9 @@ void runPoissonROM(MGmolInterface *mgmol_) rank, sampled_row[s], rho->rho_[0][sampled_row[s]], sample_rho(s), abs(rho->rho_[0][sampled_row[s]] - sample_rho(s))); } - // sample_rho.gather(); - // CAROM::Vector *rom_rho = pot_rhs_rom.mult(sample_rho); + sample_rho.gather(); + CAROM::Vector rom_rho(num_pot_basis, false); + pot_rho_rom.mult(sample_rho, rom_rho); // if (rank == 0) // { // printf("rom rho before projection\n"); @@ -459,38 +476,44 @@ void runPoissonROM(MGmolInterface *mgmol_) // if (rank == 0) // printf("rank %d, rho scaler: %.3e\n", rank, nel / tcharge); - // /* check rho projection */ - // CAROM::Vector fom_rho_vec(&rho->rho_[0][0], dim, true, false); - // CAROM::Vector *rho_proj = pot_basis->transposeMult(fom_rho_vec); + /* check rho projection */ + CAROM::Vector fom_rho_vec(&rho->rho_[0][0], dim, true, false); + CAROM::Vector rho_proj(num_pot_basis, false); + pot_basis->transposeMult(fom_rho_vec, rho_proj); // CAROM::Vector *rho_proj2 = pot_basis->mult(*rho_proj); // (*rho_proj2) -= fom_rho_vec; // double rho_proj_error = rho_proj2->norm() / fom_rho_vec.norm(); - // if (rank == 0) - // { - // printf("rom rho\n"); - // for (int d = 0; d < num_pot_basis; d++) - // printf("%.5e\t", rom_rho->item(d)); - // printf("\n"); - // printf("fom rho projection\n"); - // for (int d = 0; d < num_pot_basis; d++) - // printf("%.5e\t", rho_proj->item(d)); - // printf("\n"); + if (rank == 0) + { + printf("rom rho\n"); + for (int d = 0; d < num_pot_basis; d++) + printf("%.5e\t", rom_rho.item(d)); + printf("\n"); + printf("fom rho projection\n"); + for (int d = 0; d < num_pot_basis; d++) + printf("%.5e\t", rho_proj.item(d)); + printf("\n"); - // printf("rho proj error: %.5e\n", rho_proj_error); - // } + // printf("rho proj error: %.5e\n", rho_proj_error); + } - /* compute ROM ion density using hyper reduction */ - std::shared_ptr ions = mgmol->getIons(); - CAROM::Vector sampled_rhoc(1, true); // this will be resized in evalIonDensityOnSamplePts - pot.evalIonDensityOnSamplePts(*ions, sampled_row, sampled_rhoc); + /* project FOM ion density onto potential pod basis */ + CAROM::Vector fom_rhoc_vec(pot.rho_comp(), dim, true, false); + CAROM::Vector rom_rhoc(num_pot_basis, false); + pot_basis->transposeMult(fom_rhoc_vec, rom_rhoc); - /* check sampled rhoc */ - RHODTYPE *rho_comp = pot.rho_comp(); - for (int s = 0; s < sampled_row.size(); s++) - { - printf("rank %d, rhoc[%d]: %.5e, sampled_rhoc: %.5e, diff: %.5e\n", - rank, sampled_row[s], rho_comp[sampled_row[s]], sampled_rhoc(s), abs(rho_comp[sampled_row[s]] - sampled_rhoc(s))); - } + // /* compute ROM ion density using hyper reduction */ + // std::shared_ptr ions = mgmol->getIons(); + // CAROM::Vector sampled_rhoc(1, true); // this will be resized in evalIonDensityOnSamplePts + // pot.evalIonDensityOnSamplePts(*ions, sampled_row, sampled_rhoc); + + // /* check sampled rhoc */ + // RHODTYPE *rho_comp = pot.rho_comp(); + // for (int s = 0; s < sampled_row.size(); s++) + // { + // printf("rank %d, rhoc[%d]: %.5e, sampled_rhoc: %.5e, diff: %.5e\n", + // rank, sampled_row[s], rho_comp[sampled_row[s]], sampled_rhoc(s), abs(rho_comp[sampled_row[s]] - sampled_rhoc(s))); + // } // sampled_rhoc.gather(); // CAROM::Vector *rom_rhoc = pot_rhs_rom.mult(sampled_rhoc); @@ -503,10 +526,10 @@ void runPoissonROM(MGmolInterface *mgmol_) // if (rank == 0) // printf("rank %d, rhoc scaler: %.3e\n", rank, ionic_charge / comp_rho); - /* right-hand side */ - CAROM::Vector sampled_rhs(sample_rho); - sampled_rhs -= sampled_rhoc; - sampled_rhs *= (4.0 * M_PI); + // /* right-hand side */ + // CAROM::Vector sampled_rhs(sample_rho); + // sampled_rhs -= sampled_rhoc; + // sampled_rhs *= (4.0 * M_PI); /* get poisson right-hand side by applying Laplace operator to potential */ pb::GridFunc fomsol_gf(grid, bc[0], bc[1], bc[2]); @@ -516,27 +539,27 @@ void runPoissonROM(MGmolInterface *mgmol_) poisson->applyOperator(fomsol_gf, fomrhs_gf); CAROM::Vector fomrhs(pot_basis->numRows(), true); fomrhs_gf.init_vect(fomrhs.getData(), 'd'); - /* check sampled rhs */ - for (int s = 0; s < sampled_row.size(); s++) - { - printf("rank %d, rhs[%d]: %.5e, sampled_rhs: %.5e, diff: %.5e\n", - rank, sampled_row[s], fomrhs(sampled_row[s]), sampled_rhs(s), abs(fomrhs(sampled_row[s]) - sampled_rhs(s))); - } + // /* check sampled rhs */ + // for (int s = 0; s < sampled_row.size(); s++) + // { + // printf("rank %d, rhs[%d]: %.5e, sampled_rhs: %.5e, diff: %.5e\n", + // rank, sampled_row[s], fomrhs(sampled_row[s]), sampled_rhs(s), abs(fomrhs(sampled_row[s]) - sampled_rhs(s))); + // } /* ROM right-hand side */ - sampled_rhs.gather(); - CAROM::Vector *rom_rhs = pot_rhs_rom.mult(sampled_rhs); - // *rom_rho -= (*rom_rhoc); - // *rom_rho *= (4.0 * M_PI); + CAROM::Vector rom_rhs(rom_rho); + rom_rhs -= rom_rhoc; + rom_rhs *= 4.0 * M_PI; /* solve Poisson ROM */ - CAROM::Vector *rom_pot = pot_rom_inv.mult(*rom_rhs); + CAROM::Vector rom_pot(num_pot_basis, false); + pot_rom_inv.mult(rom_rhs, rom_pot); /* data array to lift up rom solution */ std::vector test_sol(dim); /* get librom view-vector of test_sol[s] */ CAROM::Vector test_sol_vec(test_sol.data(), dim, true, false); - pot_basis->mult(*rom_pot, test_sol_vec); + pot_basis->mult(rom_pot, test_sol_vec); /* mgmol grid function for lifted-up fom solution */ pb::GridFunc testsol_gf(grid, bc[0], bc[1], bc[2]); @@ -549,30 +572,28 @@ void runPoissonROM(MGmolInterface *mgmol_) /* librom view vector for fom solution */ CAROM::Vector fom_sol_vec(pot.vh_rho(), dim, true, false); - CAROM::Vector *fom_proj = pot_basis->transposeMult(fom_sol_vec); + CAROM::Vector fom_proj(num_pot_basis, false); + pot_basis->transposeMult(fom_sol_vec, fom_proj); if (rank == 0) { printf("rom vector\n"); for (int d = 0; d < num_pot_basis; d++) - printf("%.5e\t", rom_pot->item(d)); + printf("%.5e\t", rom_pot.item(d)); printf("\n"); printf("fom projection\n"); for (int d = 0; d < num_pot_basis; d++) - printf("%.5e\t", fom_proj->item(d)); + printf("%.5e\t", fom_proj.item(d)); printf("\n"); printf("ratio\n"); for (int d = 0; d < num_pot_basis; d++) - printf("%.5e\t", rom_pot->item(d) / fom_proj->item(d)); + printf("%.5e\t", rom_pot.item(d) / fom_proj.item(d)); printf("\n"); } /* clean up */ - delete rom_rhs; + // delete rom_rhs; // delete rom_rho; // delete rom_rhoc; - // delete rom_pot; - // delete fom_proj; - // delete rho_proj; // delete rho_proj2; }