Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Modifications to elliptic eigenproblem example #280

Merged
merged 28 commits into from
Jun 24, 2024
Merged

Modifications to elliptic eigenproblem example #280

Changes from 10 commits
Commits
Show all changes
28 commits
Select commit Hold shift + click to select a range
ff004a4
Add boundary condition switch
siuwuncheung Apr 4, 2024
af1f38d
Output eigenvalue absolute error
siuwuncheung Apr 4, 2024
861bcdb
Make eigenfunction sign in FOM and ROM consistent
siuwuncheung Apr 4, 2024
e60ca5d
Clean up problems
siuwuncheung Apr 23, 2024
c523237
Change Gaussian radius
siuwuncheung Apr 23, 2024
333103d
Change calculation of Gaussian width
siuwuncheung Apr 24, 2024
0e1f903
Slight changes in example command lines
siuwuncheung May 22, 2024
8239b14
Change filenames
siuwuncheung May 22, 2024
157a2e9
Remove TODO
siuwuncheung May 22, 2024
6c7a0f8
Astyle
siuwuncheung May 22, 2024
8ab3481
Change filenames. Change calculation of Gaussian width
siuwuncheung May 30, 2024
e2add1b
Slight changes to filenames
siuwuncheung May 30, 2024
8e17778
Consistent sign over parameters and eigenvector ordering
siuwuncheung May 30, 2024
69812cf
Astyle
siuwuncheung May 30, 2024
958bd25
Fix visualization signs
siuwuncheung May 30, 2024
ec0d4c1
Fix stream names
siuwuncheung Jun 13, 2024
9e78052
Reset stream name
siuwuncheung Jun 13, 2024
1e864d6
Reorganize reference mode sign logic
siuwuncheung Jun 13, 2024
95d9533
Clean up
siuwuncheung Jun 13, 2024
8760daf
Change variable names
siuwuncheung Jun 13, 2024
1ebd1f2
Warning meesages
siuwuncheung Jun 13, 2024
91fba19
Astyle
siuwuncheung Jun 13, 2024
11d63b0
Normalize
siuwuncheung Jun 13, 2024
a819f82
Adjust Gaussian
siuwuncheung Jun 19, 2024
b4fc1fd
Rename variables
siuwuncheung Jun 19, 2024
88a9368
Remove verbose
siuwuncheung Jun 20, 2024
60f0013
Update example in header
siuwuncheung Jun 20, 2024
4d0abab
Fix failing CI
siuwuncheung Jun 20, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
192 changes: 54 additions & 138 deletions examples/prom/elliptic_eigenproblem_global_rom.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -32,7 +32,7 @@
// Merge phase: elliptic_eigenproblem_global_rom -merge -p 2 -ns 4
//
// FOM run (for error calculation):
// elliptic_eigenproblem_global_rom -fom -p 2 -a 0.5
// elliptic_eigenproblem_global_rom -fom -p 2 -a 0.50
// Example output:
// Number of unknowns: 289
// Eigenvalue 0: 0.04533314
Expand All @@ -43,7 +43,7 @@
// Elapsed time for assembling FOM: 1.471708e-03 second
// Elapsed time for solving FOM: 3.416310e-01 second
//
// Online phase: elliptic_eigenproblem_global_rom -online -p 2 -a 0.5
// Online phase: elliptic_eigenproblem_global_rom -online -p 2 -a 0.50
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I am getting the following dramatically different outputs in Ruby for ROMs:

Number of unknowns: 289
Opening file: elliptic_eigenproblem_basis.000000
spatial basis dimension is 289 x 14
Eigenvalue 0: = 26.893636
Eigenvalue 1: = 67.093674
Eigenvalue 2: = 69.162828
Eigenvalue 3: = 109.44197
Eigenvalue 4: = 136.84225
norm_ref[i] = 273.80793
norm_ev[i] = 270.02649
ip[i] = 1.8178431
norm_ref[i] = 259.39492
norm_ev[i] = 264.66499
ip[i] = -12.477348
norm_ref[i] = 264.55127
norm_ev[i] = 264.93123
ip[i] = -22.66543
norm_ref[i] = 264.64771
norm_ev[i] = 259.31325
ip[i] = -97.06303
norm_ref[i] = 269.74102
norm_ev[i] = 273.78523
ip[i] = -18.535971
FOM solution for eigenvalue 0 = 26.893632
ROM solution for eigenvalue 0 = 26.893636
Absolute error of ROM solution for eigenvalue 0 = 3.97966e-06
Relative error of ROM solution for eigenvalue 0 = 1.4797778e-07
FOM solution for eigenvalue 1 = 67.093653
ROM solution for eigenvalue 1 = 67.093674
Absolute error of ROM solution for eigenvalue 1 = 2.1396013e-05
Relative error of ROM solution for eigenvalue 1 = 3.1889772e-07
FOM solution for eigenvalue 2 = 69.162823
ROM solution for eigenvalue 2 = 69.162828
Absolute error of ROM solution for eigenvalue 2 = 4.853771e-06
Relative error of ROM solution for eigenvalue 2 = 7.0178902e-08
FOM solution for eigenvalue 3 = 109.44196
ROM solution for eigenvalue 3 = 109.44197
Absolute error of ROM solution for eigenvalue 3 = 4.1389821e-06
Relative error of ROM solution for eigenvalue 3 = 3.7818969e-08
FOM solution for eigenvalue 4 = 136.84157
ROM solution for eigenvalue 4 = 136.84225
Absolute error of ROM solution for eigenvalue 4 = 0.0006786573
Relative error of ROM solution for eigenvalue 4 = 4.9594382e-06
Relative l2 error of ROM eigenvector 0 = 0.51582637
Relative l2 error of ROM eigenvector 1 = 1.9924533
Relative l2 error of ROM eigenvector 2 = 1.8788293
Relative l2 error of ROM eigenvector 3 = 1.9193804
Relative l2 error of ROM eigenvector 4 = 0.18407987
Unable to connect to GLVis server at localhost:19916
GLVis visualization disabled.
Elapsed time for assembling ROM: 1.196660e-02 second
Elapsed time for solving ROM: 5.293224e-02 second

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I am getting the same results as you. I think the sample output comments in the header are old and haven't been updated with new values from the changes made in this PR.

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

If this is correct, then the relative l2 error for eigenvectors are pretty off.

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I am updating some examples and parameters with outputing more information. Eventually, I want to use this example:

./elliptic_eigenproblem_global_rom -offline -p 2 -id 0 -a 0.0
./elliptic_eigenproblem_global_rom -offline -p 2 -id 1 -a 1.0
./elliptic_eigenproblem_global_rom -merge -p 2 -ns 2
./elliptic_eigenproblem_global_rom -fom -p 2 -a 0.5 -visit
./elliptic_eigenproblem_global_rom -online -p 2 -a 0.5 -visit

With the most recent commit 11d63b0, the results are

Warning: eigenvector 1 in FOM and ROM are not directly comparable.
Inner product = 0.78603098
TODO: Projection error of eigenvector.
Warning: eigenvector 2 in FOM and ROM are not directly comparable.
Inner product = 0.78603092
TODO: Projection error of eigenvector.
FOM solution for eigenvalue 0 = 19.921322
ROM solution for eigenvalue 0 = 19.92139
Absolute error of ROM solution for eigenvalue 0 = 6.7559202e-05
Relative error of ROM solution for eigenvalue 0 = 3.391301e-06
FOM solution for eigenvalue 1 = 52.415509
ROM solution for eigenvalue 1 = 52.416907
Absolute error of ROM solution for eigenvalue 1 = 0.0013989485
Relative error of ROM solution for eigenvalue 1 = 2.6689592e-05
FOM solution for eigenvalue 2 = 52.415509
ROM solution for eigenvalue 2 = 52.416908
Absolute error of ROM solution for eigenvalue 2 = 0.0013990685
Relative error of ROM solution for eigenvalue 2 = 2.6691881e-05
FOM solution for eigenvalue 3 = 81.244631
ROM solution for eigenvalue 3 = 81.245048
Absolute error of ROM solution for eigenvalue 3 = 0.00041707441
Relative error of ROM solution for eigenvalue 3 = 5.1335627e-06
Relative l2 error of ROM eigenvector 0 = 0.0010013981
Relative l2 error of ROM eigenvector 1 = 0.65416973
Relative l2 error of ROM eigenvector 2 = 0.65416982
Relative l2 error of ROM eigenvector 3 = 0.0025055664

Right now, the L2 error of the eigenvectors with non-simple eigenvalues (eigenvector 1 and 2 above) are large because of the way it is calculated, which is direct comparison of the FOM/ROM eigenvectors of the same index.

I will update the results in the header of the script and on the libROM webpage soon after the next PR.

// Example output:
// Eigenvalue 0: = 0.048430949
// Eigenvalue 1: = 0.12021157
Expand Down Expand Up @@ -74,6 +74,7 @@ double Conductivity(const Vector &x);
double Potential(const Vector &x);
int problem = 1;
double kappa;
double gaussian_width;
Vector bb_min, bb_max;

int main(int argc, char *argv[])
Expand All @@ -88,6 +89,7 @@ int main(int argc, char *argv[])
const char *mesh_file = "";
int ser_ref_levels = 2;
int par_ref_levels = 1;
bool dirichlet = true;
int order = 1;
int nev = 5;
int seed = 75;
Expand Down Expand Up @@ -117,6 +119,8 @@ int main(int argc, char *argv[])
"Mesh file to use.");
args.AddOption(&problem, "-p", "--problem",
"Problem setup to use. See options in Conductivity and Potential functions.");
args.AddOption(&dirichlet, "-dir", "--dirichlet", "-neu", "--neumann",
"BC switch.");
args.AddOption(&ser_ref_levels, "-rs", "--refine-serial",
"Number of times to refine the mesh uniformly in serial.");
args.AddOption(&par_ref_levels, "-rp", "--refine-parallel",
Expand Down Expand Up @@ -194,8 +198,6 @@ int main(int argc, char *argv[])
args.PrintOptions(cout);
}

kappa = alpha * M_PI;

if (fom)
{
MFEM_VERIFY(fom && !offline && !online
Expand All @@ -214,29 +216,16 @@ int main(int argc, char *argv[])
Mesh *mesh;
if (mesh_file == "")
{
// square domain with side length 20 at origin (0,0)
double x_max = 20.0;
double y_max = 20.0;
mesh = new Mesh(Mesh::MakeCartesian2D(2, 2, Element::QUADRILATERAL, false,
x_max, y_max));
mesh->GetBoundingBox(bb_min, bb_max, max(order, 1));

// shift mesh around origin (from [0, bb_max] -> [-bb_max/2, bb_max/2])
auto *mesh_transform = +[](const Vector &v_orig, Vector &v_new) -> void {
v_new = v_orig;
// shift mesh vertices by (bb_max[i] / 2)
v_new(0) -= 0.5*bb_max[0];
v_new(1) -= 0.5*bb_max[1];
};

// performs the transform - bb_min/bb_max are updated again below
mesh->Transform(mesh_transform);
mesh = new Mesh(Mesh::MakeCartesian2D(2, 2, Element::QUADRILATERAL));
}
else
{
mesh = new Mesh(mesh_file, 1, 1);
}
int dim = mesh->Dimension();
mesh->GetBoundingBox(bb_min, bb_max, max(order, 1));
gaussian_width = mesh->GetElementVolume(0) * pow(2.0,
1.0 - ser_ref_levels - par_ref_levels);

// 4. Refine the mesh in serial to increase the resolution. In this example
// we do 'ser_ref_levels' of uniform refinement, where 'ser_ref_levels' is
Expand All @@ -245,7 +234,6 @@ int main(int argc, char *argv[])
{
mesh->UniformRefinement();
}
mesh->GetBoundingBox(bb_min, bb_max, max(order, 1));

// 5. Define a parallel mesh by a partitioning of the serial mesh. Refine
// this mesh further in parallel to increase the resolution. Once the
Expand Down Expand Up @@ -284,8 +272,8 @@ int main(int argc, char *argv[])
int max_num_snapshots = 100;
bool update_right_SV = false;
bool isIncremental = false;
const std::string basisName = "basis";
const std::string basisFileName = basisName + std::to_string(id);
const std::string baseName = "elliptic_eigenproblem_";
const std::string basis_filename = baseName + "basis";
CAROM::Options* options;
CAROM::BasisGenerator *generator;
StopWatch solveTimer, assembleTimer, mergeTimer;
Expand All @@ -295,7 +283,9 @@ int main(int argc, char *argv[])
{
options = new CAROM::Options(fespace->GetTrueVSize(), nev, nev,
update_right_SV);
generator = new CAROM::BasisGenerator(*options, isIncremental, basisFileName);
std::string snapshot_basename = baseName + "par" + std::to_string(id);
generator = new CAROM::BasisGenerator(*options, isIncremental,
snapshot_basename);
}

// 9. The merge phase
Expand All @@ -304,10 +294,10 @@ int main(int argc, char *argv[])
mergeTimer.Start();
options = new CAROM::Options(fespace->GetTrueVSize(), max_num_snapshots, nev,
update_right_SV);
generator = new CAROM::BasisGenerator(*options, isIncremental, basisName);
generator = new CAROM::BasisGenerator(*options, isIncremental, basis_filename);
for (int paramID=0; paramID<nsets; ++paramID)
{
std::string snapshot_filename = basisName + std::to_string(
std::string snapshot_filename = baseName + "par" + std::to_string(
paramID) + "_snapshot";
generator->loadSamples(snapshot_filename, "snapshot");
}
Expand Down Expand Up @@ -335,10 +325,11 @@ int main(int argc, char *argv[])
assembleTimer.Start();
ConstantCoefficient one(1.0);

kappa = alpha * M_PI;
FunctionCoefficient kappa_0(Conductivity);
FunctionCoefficient v_0(Potential);

// Project initial conductivity and potential to visualize initial condition
// Project conductivity and potential to visualize initial condition
ParGridFunction c_gf(fespace);
c_gf.ProjectCoefficient(kappa_0);

Expand All @@ -349,7 +340,7 @@ int main(int argc, char *argv[])
if (pmesh->bdr_attributes.Size())
{
ess_bdr.SetSize(pmesh->bdr_attributes.Max());
ess_bdr = 1;
ess_bdr = (dirichlet) ? 1 : 0;
}

ParBilinearForm *a = new ParBilinearForm(fespace);
Expand Down Expand Up @@ -492,10 +483,11 @@ int main(int argc, char *argv[])
}

// 15. The online phase
Vector sign_ev(nev);
if (online) {
// 16. read the reduced basis
assembleTimer.Start();
CAROM::BasisReader reader(basisName);
CAROM::BasisReader reader(basis_filename);

Vector ev;
const CAROM::Matrix *spatialbasis;
Expand All @@ -514,23 +506,23 @@ int main(int argc, char *argv[])
numColumnRB);

// 17. form ROM operator
CAROM::Matrix invReducedA;
ComputeCtAB(*A, *spatialbasis, *spatialbasis, invReducedA);
DenseMatrix *A_mat = new DenseMatrix(invReducedA.numRows(),
invReducedA.numColumns());
A_mat->Set(1, invReducedA.getData());

CAROM::Matrix invReducedM;
ComputeCtAB(*M, *spatialbasis, *spatialbasis, invReducedM);
DenseMatrix *M_mat = new DenseMatrix(invReducedM.numRows(),
invReducedM.numColumns());
M_mat->Set(1, invReducedM.getData());
CAROM::Matrix ReducedA;
ComputeCtAB(*A, *spatialbasis, *spatialbasis, ReducedA);
DenseMatrix *A_mat = new DenseMatrix(ReducedA.numRows(),
ReducedA.numColumns());
A_mat->Set(1, ReducedA.getData());

CAROM::Matrix ReducedM;
ComputeCtAB(*M, *spatialbasis, *spatialbasis, ReducedM);
DenseMatrix *M_mat = new DenseMatrix(ReducedM.numRows(),
ReducedM.numColumns());
M_mat->Set(1, ReducedM.getData());

assembleTimer.Stop();

// 18. solve ROM
solveTimer.Start();
// (Q^T A Q) c = \lamba (Q^T M Q) c
// (Q^T A Q) c = \lambda (Q^T M Q) c
A_mat->Eigenvalues(*M_mat, ev, evect);
solveTimer.Stop();

Expand Down Expand Up @@ -641,12 +633,14 @@ int main(int argc, char *argv[])
for (int i = 0; i < eigenvalues.Size() && i < nev; i++)
{
diff_ev[i] = ev_fom[i] - eigenvalues[i];
double ev_diff_norm = sqrt(diff_ev[i] * diff_ev[i]);
double ev_fom_norm = sqrt(ev_fom[i] * ev_fom[i]);
if (myid == 0)
{
std::cout << "FOM solution for eigenvalue " << i << " = " <<
ev_fom[i] << std::endl;
std::cout << "Absolute error of ROM solution for eigenvalue " << i << " = " <<
abs(diff_ev[i]) << std::endl;
std::cout << "Relative error of ROM solution for eigenvalue " << i << " = " <<
ev_diff_norm / ev_fom_norm << std::endl;
abs(diff_ev[i]) / abs(ev_fom[i]) << std::endl;
}
}

Expand All @@ -672,7 +666,8 @@ int main(int argc, char *argv[])
mode_fom_ifs.close();

const double fomNorm = sqrt(InnerProduct(mode_fom, mode_fom));
mode_fom -= mode_rom;
sign_ev[i] = (InnerProduct(mode_fom, mode_rom) >= 0) ? 1 : -1;
mode_fom.Add(-sign_ev[i], mode_rom);
const double diffNorm = sqrt(InnerProduct(mode_fom, mode_fom));
if (myid == 0) std::cout << "Relative l2 error of ROM eigenvector " << i <<
" = " << diffNorm /
Expand All @@ -686,8 +681,8 @@ int main(int argc, char *argv[])
VisItDataCollection visit_dc("EllipticEigenproblem", pmesh);
if (visit)
{
visit_dc.RegisterField("InitialConductivity", &c_gf);
visit_dc.RegisterField("InitialPotential", &p_gf);
visit_dc.RegisterField("Conductivity", &c_gf);
visit_dc.RegisterField("Potential", &p_gf);
std::vector<ParGridFunction*> visit_evs;
for (int i = 0; i < nev && i < eigenvalues.Size(); i++)
{
Expand All @@ -699,6 +694,7 @@ int main(int argc, char *argv[])
Vector ev;
evect.GetRow(i, ev);
x = ev;
x *= sign_ev[i];
}
visit_evs.push_back(new ParGridFunction(x));
visit_dc.RegisterField("x" + std::to_string(i), visit_evs.back());
Expand Down Expand Up @@ -750,6 +746,7 @@ int main(int argc, char *argv[])
Vector ev;
evect.GetRow(i, ev);
x = ev;
x *= sign_ev[i];
}

sout << "parallel " << num_procs << " " << myid << "\n"
Expand Down Expand Up @@ -827,113 +824,32 @@ double Conductivity(const Vector &x)
return 1.0 + cos(kappa * x(1)) * sin(kappa * x(0));
case 3:
case 4:
case 5:
case 6:
case 7:
case 8:
case 9:
return 1.0;
}
return 0.0;
}

double Potential(const Vector &x)
{
int dim = x.Size();

Vector X(dim);
Vector center(dim); // center of gaussian for problem 4 and 5
center = kappa / M_PI;

Vector neg_center(center);
neg_center.Neg();

// amplitude of gaussians for problems 4-6
const double D = 100.0;
// width of gaussians for problems 4-6
const double c = 0.05;

const double min_d = 2.0 * sqrt(2.0);

X = x;
for (int i = 0; i < dim; i++)
{
// map alpha parameter from [-1,1] to the mesh bounding box (controls the center for problems 4 and 5)
center(i) = CAROM::map_to_ref_mesh(bb_min[i], bb_max[i], center(i));
neg_center(i) = CAROM::map_to_ref_mesh(bb_min[i], bb_max[i], neg_center(i));
}

double D = 100.0;
Vector center(x.Size());
switch (problem)
{
case 1:
case 2:
return 0.0;
case 3:
return 1.0;
case 4:
return D * std::exp(-X.DistanceSquaredTo(center) / c);
case 5:
return -D * std::exp(-X.DistanceSquaredTo(center) / c);
case 6:
return -D * (std::exp(-X.DistanceSquaredTo(center) / c) + std::exp(
-X.DistanceSquaredTo(neg_center) / c));
case 7:
center = 0.25;
center(0) += 0.05 * cos(2.0*kappa);
center(1) += 0.05 * sin(2.0*kappa);

neg_center = center;
neg_center.Neg();

for (int i = 0; i < dim; i++)
for (int i = 0; i < x.Size(); i++)
{
// map alpha parameter from [-1,1] to the mesh bounding box (controls the center for problems 4 and 5)
center(i) = CAROM::map_to_ref_mesh(bb_min[i], bb_max[i], center(i));
neg_center(i) = CAROM::map_to_ref_mesh(bb_min[i], bb_max[i], neg_center(i));
center(i) = 0.5 * (bb_min[i] + bb_max[i]);
}
return -D * (std::exp(-X.DistanceSquaredTo(center) / c) + std::exp(
-X.DistanceSquaredTo(neg_center) / c));
case 8:
// Similar to case 6, but t is restricted to inner 20% of the domain
// in this case, the radius of the gaussian is (0.05*min_d)^2 where
// min_d is the lower bound of the atom distance over time

// verify t is within inner 20% of domain (centered around mesh origin)
CAROM::verify_within_portion(bb_min, bb_max, center, 0.2);

return -D * (std::exp(-X.DistanceSquaredTo(center) / std::pow(c * min_d,
2)) + std::exp(
-X.DistanceSquaredTo(neg_center) / std::pow(c * min_d, 2)));
case 9:
// Similar to case 7, but t is restricted to inner 20% of the domain and t is defined as:
// t = (1.5 + 0.5*cos(2*k), 1.5 + 0.5*sin(2*k)) where k = alpha * PI, alpha is the input parameter given with the -a option.
// The radius of the gaussian follows case 8: (0.05*min_d)^2

// Sets the center to vary around +/- 1.5 (absolute location on the mesh)
center = 1.5;
center(0) += 0.5 * cos(2.0*kappa);
center(1) += 0.5 * sin(2.0*kappa);

// map the absolute location back to a fraction of the mesh domain
center(0) = CAROM::map_from_ref_mesh(bb_min[0], bb_max[0], center(0));
center(1) = CAROM::map_from_ref_mesh(bb_min[1], bb_max[1], center(1));

neg_center = center;
neg_center.Neg();

for (int i = 0; i < dim; i++)
return D * std::exp(-x.DistanceSquaredTo(center) / pow(gaussian_width, 2.0));
case 4:
for (int i = 0; i < x.Size(); i++)
{
// map alpha parameter from [-1,1] to the mesh bounding box (controls the center for problems 4 and 5)
center(i) = CAROM::map_to_ref_mesh(bb_min[i], bb_max[i], center(i));
neg_center(i) = CAROM::map_to_ref_mesh(bb_min[i], bb_max[i], neg_center(i));
center(i) = 0.5 * (bb_min[i] + bb_max[i]);
}

// verify t is within inner 20% of domain (centered around mesh origin)
CAROM::verify_within_portion(bb_min, bb_max, center, 0.2);

return -D * (std::exp(-X.DistanceSquaredTo(center) / std::pow(c * min_d,
2)) + std::exp(
-X.DistanceSquaredTo(neg_center) / std::pow(c * min_d, 2)));
return D / x.DistanceSquaredTo(center);
}
return 0.0;
}
Loading