Skip to content

Commit c2c7fd5

Browse files
cval26cval26
committed
fix variables offset
1 parent 679c952 commit c2c7fd5

2 files changed

Lines changed: 77 additions & 59 deletions

File tree

‎rom/laghos_rom.cpp‎

Lines changed: 73 additions & 58 deletions
Original file line numberDiff line numberDiff line change
@@ -57,9 +57,10 @@ void ROM_Sampler::SampleSolution(const double t, const double dt, const double p
5757
SetStateVariables(S);
5858
SetStateVariableRates(dt);
5959

60-
const bool sampleX = generator_X->isNextSample(t);
61-
62-
Vector dSdt;
60+
const bool sampleX = (hyperreductionSamplingType == eqp_energy) ? true :
61+
generator_X->isNextSample(t);
62+
63+
Vector dSdt;
6364
if (!sns && rhsBasis)
6465
{
6566
dSdt.SetSize(S.Size());
@@ -98,7 +99,8 @@ void ROM_Sampler::SampleSolution(const double t, const double dt, const double p
9899
}
99100
}
100101

101-
const bool sampleV = generator_V->isNextSample(t);
102+
const bool sampleV = (hyperreductionSamplingType == eqp_energy) ? true :
103+
generator_V->isNextSample(t);
102104

103105
//TODO: use this, plus generator_Fv->computeNextSampleTime? So far, it seems sampleV == true on every step.
104106
//const bool sampleFv = generator_Fv->isNextSample(t);
@@ -152,7 +154,8 @@ void ROM_Sampler::SampleSolution(const double t, const double dt, const double p
152154
}
153155
}
154156

155-
const bool sampleE = generator_E->isNextSample(t);
157+
const bool sampleE = (hyperreductionSamplingType == eqp_energy) ? true :
158+
generator_E->isNextSample(t);
156159

157160
if (sampleE)
158161
{
@@ -672,7 +675,7 @@ void ROM_Sampler::SetupEQP_Force_Eq(const CAROM::Matrix* snapX,
672675
} // j
673676
} // i
674677

675-
CAROM::Vector w(ne * nqe, true);
678+
CAROM::Vector w(ne * nqe, true);
676679

677680
for (int i=0; i<ne; ++i)
678681
{
@@ -724,13 +727,13 @@ void ROM_Sampler::SetupEQP_En_Force_Eq(const CAROM::Matrix* snapX,
724727
// G is the matrix of accuracy constraints used to enforce that the
725728
// evaluated quantity remains close to the result of the full quadrature rule.
726729

727-
// Declare G of size ((NBv + NBe) * (nsnap+1)) x NQ; store its transpose Gt.
728-
// The first NBv*(nsnap+1) rows of G hold the velocity constraints;
729-
// the remaining NBe*(nsnap+1) rows hold the energy constraints.
730-
CAROM::Matrix Gt(NQ, (NBv + NBe) * (nsnap+1), true);
730+
// Declare G of size ((NBv + NBe) * nsnap) x NQ; store its transpose Gt.
731+
// The first NBv * nsnap rows of G hold the velocity constraints;
732+
// the remaining NBe * nsnap rows hold the energy constraints.
733+
CAROM::Matrix Gt(NQ, (NBv + NBe) * nsnap, true);
731734

732735
// row index of G where energy constraints start
733-
const int estart = NBv * (nsnap + 1);
736+
const int estart = NBv * nsnap;
734737

735738
Vector v_i(tH1size), x_i(tH1size), e_i(tL2size);
736739
Vector w_j_e, v_i_e, v_j_e;
@@ -763,43 +766,35 @@ void ROM_Sampler::SetupEQP_En_Force_Eq(const CAROM::Matrix* snapX,
763766

764767
ParGridFunction gf2H1(gfH1);
765768

766-
for (int i=0; i<nsnap+1; ++i)
769+
for (int i=0; i<nsnap; ++i)
767770
{
768-
if (i == 0) // Use the initial state as the first snapshot.
769-
{
770-
v_i = 0.0;
771-
x_i = 0.0;
772-
e_i = 0.0;
773-
}
774-
else
771+
if (i == 0)
775772
{
776-
if (i == 1 && numSkipped[0] == 1)
777-
{
773+
if (numSkipped[0] == 1)
778774
x_i = 0.0;
779-
}
780775
else
781-
{
782-
for (int j = 0; j < tH1size; ++j)
783-
x_i[j] = (*snapX)(j, i - 1 - numSkipped[0]);
784-
}
776+
for (int j = 0; j < tH1size; ++j) x_i[j] = (*snapX)(j, 0);
785777

786-
if (i == 1 && numSkipped[1] == 1)
778+
if (numSkipped[1] == 1)
787779
v_i = 0.0;
788780
else
789-
{
790-
for (int j = 0; j < tH1size; ++j)
791-
v_i[j] = (*snapV)(j, i - 1 - numSkipped[1]);
792-
}
781+
for (int j = 0; j < tH1size; ++j) v_i[j] = (*snapV)(j, 0);
793782

794-
if (i == 1 && numSkipped[2] == 1)
783+
if (numSkipped[2] == 1)
795784
e_i = 0.0;
796785
else
797-
{
798-
for (int j = 0; j < tL2size; ++j)
799-
{
800-
e_i[j] = (*snapE)(j, i - 1 - numSkipped[2]);
801-
}
802-
}
786+
for (int j = 0; j < tL2size; ++j) e_i[j] = (*snapE)(j, 0);
787+
}
788+
else
789+
{
790+
for (int j = 0; j < tH1size; ++j)
791+
x_i[j] = (*snapX)(j, i - numSkipped[0]);
792+
793+
for (int j = 0; j < tH1size; ++j)
794+
v_i[j] = (*snapV)(j, i - numSkipped[1]);
795+
796+
for (int j = 0; j < tL2size; ++j)
797+
e_i[j] = (*snapE)(j, i - numSkipped[2]);
803798
}
804799

805800
SetStateFromTrueDOFs(x_i, v_i, e_i, S);
@@ -810,7 +805,7 @@ void ROM_Sampler::SetupEQP_En_Force_Eq(const CAROM::Matrix* snapX,
810805
input.FOMoper->ResetQuadratureData();
811806

812807
// Velocity equation.
813-
// For 0 <= j < NBv, 0 <= i <= nsnap, 0 <= e < ne, 0 <= m < nqe,
808+
// For 0 <= j < NBv, 0 <= i < nsnap, 0 <= e < ne, 0 <= m < nqe,
814809
// entry G(j + (i*NBv), (e*nqe) + m) is the coefficient of
815810
//
816811
// e_j^T * F(v_i,e_i,x_i) * 1E
@@ -821,7 +816,7 @@ void ROM_Sampler::SetupEQP_En_Force_Eq(const CAROM::Matrix* snapX,
821816
// identity vector in the energy space.
822817

823818
// Energy equation.
824-
// For 0 <= j < NBe, 0 <= i <= nsnap, 0 <= e < ne, 0 <= m < nqe,
819+
// For 0 <= j < NBe, 0 <= i < nsnap, 0 <= e < ne, 0 <= m < nqe,
825820
// entry G(estart + j + (i*NBe), (e*nqe) + m) is the coefficient of
826821
//
827822
// e_j^T * F(v_i,e_i,x_i)^T * v_i
@@ -1097,17 +1092,10 @@ void ROM_Sampler::Finalize(Array<int> &cutoff, ROM_Options& input)
10971092
{
10981093
if (rank == 0)
10991094
{
1100-
// For the energy-conserving EQP case, we increase the energy basis
1095+
// For the energy-conserving EQP case, increase the energy basis
11011096
// dimension by 1 to accomodate for the addition of the energy
11021097
// identity.
1103-
// The change is done here, before the window basis parameters are
1104-
// written to their files, and is also carried to the caller's scope.
11051098

1106-
// TODO: if cutoff[2] == windowNumSamples then increasing this by
1107-
// 1 will lead to a function call trying to get more columns than
1108-
// there are available in basisE. In this case, we'd need to take
1109-
// all columns of basisE and then append one.
1110-
// And then do the same during the online stage.
11111099
cutoff[2] += 1;
11121100
}
11131101

@@ -1120,17 +1108,21 @@ void ROM_Sampler::Finalize(Array<int> &cutoff, ROM_Options& input)
11201108
// For the energy basis, the cutoff parameter has already been
11211109
// increased by 1 to accomodate the energy identity.
11221110
CAROM::Matrix *tBasisV = basisV->getFirstNColumns(cutoff[1]);
1123-
CAROM::Matrix *tBasisE = basisE->getFirstNColumns(cutoff[2]);
1124-
1125-
// Form the energy identity and replace the last basis vector by it.
1126-
Vector unitE(tL2size);
1127-
unitE = 1.0;
1111+
1112+
CAROM::Matrix *tBasisE = new CAROM::Matrix(tL2size, cutoff[2], false);
11281113

1114+
// Get the first cutoff[2]-1 columns of basisE
11291115
for (int i = 0; i < tL2size; i++)
11301116
{
1131-
(*tBasisE)(i, cutoff[2]-1) = unitE[i];
1117+
for (int j = 0; j < cutoff[2] - 1; j++)
1118+
(*tBasisE)(i, j) = (*basisE)(i, j);
11321119
}
11331120

1121+
// Form the energy identity and include it as the last basis vector.
1122+
Vector unitE(tL2size);
1123+
unitE = 1.0;
1124+
for (int i = 0; i < tL2size; i++)
1125+
(*tBasisE)(i, cutoff[2]-1) = unitE[i];
11341126
tBasisE->orthogonalize();
11351127

11361128
SetupEQP_Force(generator_X->getSnapshotMatrix(),
@@ -1288,17 +1280,18 @@ ROM_Basis::ROM_Basis(ROM_Options const& input, MPI_Comm comm_, const double sFac
12881280
mfH1.SetSize(tH1size);
12891281
mfL2.SetSize(tL2size);
12901282

1283+
if (hyperreduce && hyperreductionSamplingType == eqp_energy)
1284+
basisE = new CAROM::Matrix(tL2size, rdime, true);
1285+
12911286
ReadSolutionBases(input.window);
12921287

12931288
if (hyperreduce && hyperreductionSamplingType == eqp_energy)
12941289
{
1295-
// Form the energy identity and replace the last E-basis vector by it.
1290+
// Form the energy identity and include it as the last basis vector.
12961291
Vector unitE(tL2size);
12971292
unitE = 1.0;
12981293
for (int i = 0; i < tL2size; i++)
1299-
{
13001294
(*basisE)(i, rdime-1) = unitE[i];
1301-
}
13021295
basisE->orthogonalize();
13031296
}
13041297

@@ -2222,7 +2215,29 @@ void ROM_Basis::ReadSolutionBases(const int window)
22222215
if (!useVX)
22232216
basisV = ReadBasisROM(rank, basename + "/" + ROMBasisName::V + std::to_string(window) + basisIdentifier, tH1size, rdimv);
22242217

2225-
basisE = ReadBasisROM(rank, basename + "/" + ROMBasisName::E + std::to_string(window) + basisIdentifier, tL2size, rdime);
2218+
// In the energy-conserving EQP case we read the first rdime-1 basis
2219+
// vectors, since the rdime parameter has been increased by 1 to
2220+
// accommodate the addition of the energy identity.
2221+
if (hyperreduce && hyperreductionSamplingType == eqp_energy)
2222+
{
2223+
int tmp_rdime = rdime - 1;
2224+
CAROM::Matrix *tmp_basisE = 0;
2225+
2226+
tmp_basisE = ReadBasisROM(rank, basename + "/" + ROMBasisName::E +
2227+
std::to_string(window) + basisIdentifier, tL2size, tmp_rdime);
2228+
2229+
for (int i = 0; i < tL2size; i++)
2230+
{
2231+
for (int j = 0; j < tmp_rdime; j++)
2232+
(*basisE)(i, j) = (*tmp_basisE)(i, j);
2233+
}
2234+
}
2235+
else
2236+
{
2237+
basisE = ReadBasisROM(rank, basename + "/" + ROMBasisName::E +
2238+
std::to_string(window) + basisIdentifier, tL2size, rdime);
2239+
}
2240+
22262241

22272242
if (useXV)
22282243
basisX = basisV;

‎rom/laghos_rom.hpp‎

Lines changed: 4 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -481,7 +481,8 @@ class ROM_Sampler
481481
energyFraction_X(input.energyFraction_X), sns(input.SNS), lhoper(input.FOMoper), writeSnapshots(input.parameterID >= 0),
482482
parameterID(input.parameterID), basename(*input.basename), Voffset(!input.useXV && !input.useVX && !input.mergeXV),
483483
useXV(input.useXV), useVX(input.useVX), VTos(input.VTos), spaceTime(input.spaceTimeMethod != no_space_time),
484-
rhsBasis(input.hyperreductionSamplingType != eqp && input.hyperreductionSamplingType != eqp_energy)
484+
rhsBasis(input.hyperreductionSamplingType != eqp && input.hyperreductionSamplingType != eqp_energy),
485+
hyperreductionSamplingType(input.hyperreductionSamplingType)
485486
{
486487
const int window = input.window;
487488

@@ -685,6 +686,8 @@ class ROM_Sampler
685686

686687
const bool spaceTime;
687688

689+
const bool hyperreductionSamplingType;
690+
688691
hydrodynamics::LagrangianHydroOperator *lhoper;
689692

690693
int VTos = 0; // Velocity temporal index offset, used for V and Fe. This fixes the issue that V and Fe are not sampled at t=0, since they are initially zero. This is valid for the Sedov test but not in general when the initial velocity is nonzero.

0 commit comments

Comments
 (0)