00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012 #include <iostream>
00013 #include <fstream>
00014 #include <sstream>
00015 #include <cstdlib>
00016 #include <string>
00017 #include <cassert>
00018 #include <unistd.h>
00019
00020 #include "myrlog.h"
00021
00022 #include <ml_include.h>
00023 #include <Epetra_Operator.h>
00024 #include <Epetra_MultiVector.h>
00025 #include <Epetra_SerialDenseVector.h>
00026 #include <Epetra_Vector.h>
00027 #include <ml_epetra_preconditioner.h>
00028
00029 #include "OutputCapturer.h"
00030 #include "pbedefs.h"
00031
00032 #include "looseoctree.h"
00033 #include "debugmeshbuilder.h"
00034 #include "tetmeshbuilder.h"
00035 #include "broadcastmeshbuilder.h"
00036 #include "netgenreader.h"
00037 #include "femaxmesh.h"
00038 #include "nedelecmesh.h"
00039 #include "lagrangemesh.h"
00040 #include "export_Epetra_CrsMatrix.h"
00041 #include "operatortest.h"
00042 #include "materials.h"
00043
00044 #include "AbstractEigsolvOperators.h"
00045 #include "LinearEigsolvOperators.h"
00046 #include "QuadraticEigsolvOperators.h"
00047 #include "Projector.h"
00048
00049 #include "ModalAnalysisSolver.h"
00050 #include "femaxxdriver.h"
00051
00052 #include "JDBSYM.h"
00053 #include "KnyazevLOBPCG.h"
00054 #include "LOBPCG.h"
00055
00056 #include "CheckingTools.h"
00057 #include "vtkexport.h"
00058 #include "h5_tools.h"
00059
00060 using namespace mesh;
00061
00062 FemaxxDriver::FemaxxDriver(const Epetra_Comm& comm, Teuchos::ParameterList& params_)
00063 : comm_(comm),
00064 params_(params_),
00065 Q_(0),
00066 lambda_(0),
00067 nof_converged_(0),
00068 tmesh_(0),
00069 femax_mesh_(0)
00070 {
00071
00072
00073
00074
00075
00076
00077
00078
00079
00080
00081
00082
00083
00084
00085
00086
00087
00088
00089
00090
00091
00092
00093
00094
00095
00096
00097
00098
00099
00100
00101
00102
00103
00104
00105
00106
00107
00108
00109
00110
00111
00112
00113
00114
00115 const int nof_pbe_timers = 150;
00116 pbe_init("femaxx", nof_pbe_timers, 0, NULL);
00117 pbe_start(50, "Whole program");
00118 }
00119
00120 FemaxxDriver::~FemaxxDriver() {
00121 delete Q_;
00122 delete[] lambda_;
00123 delete femax_mesh_;
00124 delete tmesh_;
00125
00126
00127 pbe_stop(50);
00128
00129 OutputCapturer cap;
00130 cap.begin_capture();
00131 pbe_dump();
00132 pbe_finalize(0);
00133 cap.end_capture();
00134 if (comm_.MyPID() == 0) {
00135 cap.log("PBE output", RLOG_CHANNEL("info"));
00136 }
00137 }
00138
00139 void FemaxxDriver::set_defaults(Teuchos::ParameterList& params) {
00140
00141 params.set("element_order", 1);
00142 params.set("kmax", 5);
00143 params.set("eigtol", 1e-8);
00144 params.set("tau", 0.0);
00145 params.set("sym_plane_config", 0);
00146 params.set("eigsolver", "jdsym");
00147 params.sublist("asigma_precon").set("type", "ml");
00148 params.sublist("asigma_precon").set("solver11", "lu");
00149 params.sublist("asigma_precon").set("solver22", "diag");
00150 params.sublist("h_solver").set("type", "pcg");
00151 params.sublist("h_precon").set("type", "ml");
00152 params.sublist("h_precon").set("solver11", "lu");
00153 params.sublist("h_precon").set("solver22", "diag");
00154 params.set("vtk_file_name", "");
00155 params.set("export_mtx", false);
00156 params.set("debug", false);
00157 params.set("disable_log", false);
00158
00159
00160 ML_Epetra::SetDefaults("SA", params.sublist("h_precon").sublist("ml_params"));
00161
00162 params.sublist("h_precon").sublist("ml_params").set("smoother: type", "MLS");
00163
00164
00165 params.sublist("h_solver").sublist("amesos_params").set("solver_type", "Superludist");
00166 params.sublist("h_solver").sublist("amesos_params").set("MatrixType", "SPD");
00167 params.sublist("h_solver").sublist("amesos_params").set("OutputLevel", 2);
00168
00169
00170 params.sublist("asigma_precon").sublist("amesos_params").set("solver_type", "Superludist");
00171 params.sublist("asigma_precon").sublist("amesos_params").set("MatrixType", "symmetric");
00172 params.sublist("asigma_precon").sublist("amesos_params").set("OutputLevel", 2);
00173
00174
00175 ML_Epetra::SetDefaults("maxwell", params.sublist("asigma_precon").sublist("ml_params"));
00176
00177
00178 params.sublist("asigma_precon").sublist("ml_params").set("negative conductivity", true);
00179 params.sublist("asigma_precon").sublist("ml_params").set("subsmoother: type", "symmetric Gauss-Seidel");
00180 # if 0
00181
00182 params.sublist("asigma_precon").sublist("ml_params").set("repartition: enable",true);
00183
00184 params.sublist("asigma_precon").sublist("ml_params").set("repartition: node max min ratio",1.1);
00185
00186 params.sublist("asigma_precon").sublist("ml_params").set("repartition: node min per proc",200);
00187
00188 params.sublist("asigma_precon").sublist("ml_params").set("repartition: edge max min ratio",1.1);
00189
00190 params.sublist("asigma_precon").sublist("ml_params").set("repartition: edge min per proc",600);
00191 params.sublist("asigma_precon").sublist("ml_params").set("repartition: partitioner","ParMETIS");
00192 # endif
00193
00194
00195
00196
00197
00198 JDBSYM::set_defaults(params.sublist("jd_params"), params.get<int>("kmax"));
00199 params.sublist("jd_params").set("clvl", 2);
00200 }
00201
00202 void FemaxxDriver::load_mesh() {
00203 log_mem_footprint(comm_, "initial");
00204 rInfo("Read mesh file...");
00205 if (access(params_.get<string>("mesh_file_name").c_str(), R_OK)) {
00206 rError("Unable to access mesh file \"%s\".", params_.get<string>("mesh_file_name").c_str());
00207 exit(1);
00208 }
00209
00210
00211 const char* meshFileName = params_.get<string>("mesh_file_name").c_str();
00212 const char* materialFileName = params_.get<string>("material_file_name").c_str();
00213
00214
00215
00216
00217 log_mem_footprint(comm_, "before loading mesh");
00218 pbe_start(52, "Mesh generation");
00219 TetMeshBuilder* builder = new TetMeshBuilder(dynamic_cast<const Epetra_MpiComm&>(comm_).GetMpiComm());
00220
00221
00222
00223
00224
00225
00226
00227
00228
00229
00230
00231 HDF5ParallelReader* reader = new HDF5ParallelReader(meshFileName, builder, dynamic_cast<const Epetra_MpiComm&>(comm_).GetMpiComm());
00232 reader->read();
00233 delete reader;
00234
00235
00236
00237
00238
00239
00240
00241 tmesh_ = builder->get_mesh();
00242 delete builder;
00243 pbe_stop(52);
00244 if (comm_.MyPID() == 0)
00245 tmesh_->log_mesh_info();
00246 log_mem_footprint(comm_, "after loading tetmesh");
00247
00248
00249 tmesh_->load_materials(materialFileName);
00250
00251
00252 int sym_plane_config = params_.get<int>("sym_plane_config");
00253 rInfo("Calculate DOF mappings (sym_plane_config=%d)...", sym_plane_config);
00254 pbe_start(53, "Dof mappings");
00255 femax_mesh_ = new FemaxMesh(*tmesh_, params_.get<int>("element_order"), sym_plane_config);
00256 pbe_stop(53);
00257 log_mem_footprint(comm_, "after calc dof mappings");
00258 }
00259
00260 void FemaxxDriver::calculate_eigenfields() {
00261
00262
00263 AbstractEigsolvOperators* lop;
00264 if (params_.get<int>("element_order") == 1)
00265 lop = new LinearEigsolvOperators(*femax_mesh_, params_, comm_);
00266 else
00267 lop = new QuadraticEigsolvOperators(*femax_mesh_, params_, comm_);
00268 log_mem_footprint(comm_);
00269
00270 if (params_.get<bool>("debug")) {
00271 OperatorTest tester(lop->getA());
00272 tester.dump_stats(cout, "A original");
00273 tester.set_operator(lop->getM());
00274 tester.dump_stats(cout, "M original");
00275 tester.set_operator(lop->getH());
00276 tester.dump_stats(cout, "H original");
00277 }
00278
00279
00280
00281
00282 log_mem_footprint(comm_, "before mesh data structure dealloc");
00283
00284
00289
00290
00291
00292
00293 log_mem_footprint(comm_, "after mesh data structure dealloc");
00294
00295 const Epetra_Operator* A = lop->getA();
00296 const Epetra_Operator* M = lop->getM();
00297
00298 if (params_.get<double>("sigma") == 0.0)
00299 rWarning("Shift sigma is zero: Calculating preconditioner for a singular matrix.");
00300 const Epetra_Operator* Prec = lop->getAsigmaPrec(params_.get<double>("sigma"));
00301 Epetra_Operator* H = lop->getH();
00302 Epetra_Operator* Y = lop->getY();
00303 Epetra_Operator* C = lop->getC();
00304 log_mem_footprint(comm_, "after constr of precon for A-sigma*M");
00305
00306 if (params_.get<bool>("debug")) {
00307 OperatorTest tester(A);
00308 tester.dump_stats(cout, "A reordered");
00309 tester.set_operator(M);
00310 tester.dump_stats(cout, "M reordered");
00311 tester.set_operator(H);
00312 tester.dump_stats(cout, "H reordered");
00313 }
00314
00315
00316
00317
00318 int myPid = comm_.MyPID();
00319 int numProc = comm_.NumProc();
00320
00321
00322 srand((unsigned) time(0));
00323
00324
00325 initMemCounters();
00326
00327 double highMem = currentSize();
00328 highMem = (highMem < currentSize()) ? currentSize() : highMem;
00329
00330 int globalSize = (A->OperatorRangeMap()).NumGlobalElements();
00331
00332
00333 int kmax = params_.get<int>("kmax");
00334 kmax = (kmax > globalSize) ? globalSize : kmax;
00335
00336
00337 lambda_ = new (nothrow) double[kmax];
00338 assert(lambda_ != 0);
00339
00340 Q_ = new Epetra_MultiVector(A->OperatorRangeMap(), kmax, false);
00341 Q_->Random();
00342 log_mem_footprint(comm_, "after alloc of Q and lambda");
00343
00344 highMem = (highMem > currentSize()) ? highMem : currentSize();
00345
00346
00347 pbe_start(75, "Nullspace projector construction");
00348 rInfo("Calculate nullspace projector...");
00349 Projector YHC(Y, lop->getHSolver(), C);
00350 pbe_stop(75);
00351 log_mem_footprint(comm_, "after constr of projector");
00352
00353 if (params_.get<bool>("debug")) {
00354 OperatorTest tester(&YHC);
00355 tester.dump_stats(cout, "Nullspace projector");
00356
00357
00358 tester.set_operator(lop->getHSolver(), true);
00359 tester.dump_stats(cout, "H solver");
00360 tester.set_operator(Prec, true);
00361 tester.dump_stats(cout, "A-sigma*M precond");
00362 }
00363
00364
00365 pbe_start(54, "Eigenvalue solver");
00366 double eigtol = params_.get<double>("eigtol");
00367 ModalAnalysisSolver *mySolver = 0;
00368
00369 if (params_.get<string>("eigsolver") == "jdsym") {
00370 params_.sublist("jd_params").set("tol", eigtol);
00371 mySolver = new JDBSYM(comm_, A, M, Prec, &YHC, params_.sublist("jd_params"));
00372 } else if (params_.get<string>("eigsolver") == "knyazev") {
00373 const int max_iter = 300;
00374 const int clvl = 3;
00375 YHC.Apply(*Q_, *Q_);
00376 mySolver = new KnyazevLOBPCG(comm_, A, M, Prec, &YHC, eigtol, max_iter, clvl);
00377 } else if (params_.get<string>("eigsolver") == "lobpcg") {
00378
00379 const int max_iter = 300;
00380 const int clvl = 3;
00381 YHC.Apply(*Q_, *Q_);
00382 kmax = 2;
00383 mySolver = new LOBPCG(comm_, A, M, Prec, &YHC, kmax, eigtol, max_iter, clvl);
00384
00385
00386
00387
00388
00389
00390 } else {
00391 rError("Unknown eigensolver \"%s\".", params_.get<string>("eigsolver").c_str());
00392 }
00393
00394
00395 nof_converged_ = mySolver->solve(kmax, *Q_, lambda_);
00396 pbe_stop(54);
00397
00398
00399 Epetra_MultiVector Q = get_eigenvectors();
00400 Epetra_SerialDenseVector lambda = get_eigenvalues();
00401
00402
00403 if (nof_converged_ < 0) {
00404 rError("Eigensolver exited with error %d!!!", nof_converged_);
00405 } else {
00406 rInfo("Number of computed eigenmodes: %d", nof_converged_);
00407
00408 if (nof_converged_ > 0) {
00409
00410 CheckingTools myVerify(comm_);
00411
00412 double tmp = myVerify.errorOrthonormality(&Q, M);
00413 if (myPid == 0)
00414 rInfo("Maximum coefficient in matrix Q^T M Q - I = %le", tmp);
00415
00416 myVerify.errorEigenResiduals(Q, lambda_, A, M);
00417 #if 0
00418
00419 Epetra_MultiVector CtQ(C->OperatorDomainMap(), nof_converged_);
00420 (const_cast<Epetra_Operator*>(C))->SetUseTranspose(true);
00421 C->Apply(Q, CtQ);
00422
00423
00424
00425
00426
00427 double maxNormY = 0.0;
00428 Epetra_Vector Hdiag(dynamic_cast<Epetra_RowMatrix*>(H)->RowMatrixRowMap());
00429 dynamic_cast<Epetra_RowMatrix*>(H)->ExtractDiagonalCopy(Hdiag);
00430 Hdiag.NormInf(&maxNormY);
00431 double maxError = 0.0;
00432 for (int i = 0; i < nof_converged_; ++i) {
00433 double infNorm = 0.0;
00434 CtQ(i)->NormInf(&infNorm);
00435 maxError = (infNorm > maxError) ? infNorm : maxError;
00436 }
00437 if (myPid == 0) {
00438 cout << endl;
00439 cout << " Maximum M-norm of vectors Y = " << maxNormY << endl;
00440 cout << " Maximum coefficient in matrix C^T Q = " << maxError << endl;
00441 cout << endl;
00442 }
00443 #endif
00444 }
00445
00446
00447 Teuchos::ParameterList& jd_params(params_.sublist("jd_params"));
00448 if (params_.get<string>("eigsolver") == "jdsym") {
00449 double maxHighMem = 0.0;
00450 comm_.MaxAll(&highMem, &maxHighMem, 1);
00451 rDebug("Memory: High water mark in set-up (EST) : %.2lf MB", maxHighMem);
00452 rDebug("Memory requested per processor for eigenvectors (EST): %.2lf MB",
00453 Q_->GlobalLength()*kmax*sizeof(double)/1024.0/1024.0/numProc);
00454 rDebug("Memory requested per processor for working space (EST): %.2lf MB",
00455 Q_->GlobalLength()*(2*kmax + jd_params.get<int>("jmax"))*sizeof(double)/1024.0/1024.0/numProc);
00456 }
00457 mySolver->memoryInfo();
00458 mySolver->operationInfo();
00459 mySolver->timeInfo();
00460
00461 }
00462
00463
00464 log_mem_footprint(comm_, "before dealloc eigenvalue solver");
00465 delete mySolver; mySolver = 0;
00466 log_mem_footprint(comm_, "after dealloc eigenvalue solver");
00467
00468 delete lop; lop = 0;
00469 log_mem_footprint(comm_, "after dealloc operators");
00470
00471 #if 0
00472
00473 if (myPid == 0) {
00474 h5_create_empty_file("eigenmodes.h5");
00475 h5_write_param_list("eigenmodes.h5", params_);
00476 }
00477 comm_.Barrier();
00478 h5_write_eigenmodes("eigenmodes.h5", Q, lambda);
00479 #endif
00480 }
00481
00482 void FemaxxDriver::postprocess() {
00483 rAssert(tmesh_);
00484 rAssert(femax_mesh_);
00485
00486 Epetra_MultiVector Q = get_eigenvectors();
00487 Epetra_SerialDenseVector lambda = get_eigenvalues();
00488
00489 pbe_start(131, "Construction of octree");
00490 if (comm_.MyPID() == 0) {
00491 tmesh_->construct_octree();
00492
00493 }
00494
00495 pbe_stop(131);
00496
00497
00498
00499
00500 pbe_start(133, "Computation of quality factors");
00501 vector<double> q_factors;
00502 for (int k = 0; k < nof_converged_; ++ k) {
00503
00504
00505
00506 Epetra_MultiVector q_dist(View, Q, k, 1);
00507
00508
00509 int NumMyElements_target = 0;
00510 if (comm_.MyPID() == 0)
00511 NumMyElements_target = Q.GlobalLength();
00512 Epetra_Map TargetMap(-1, NumMyElements_target, 0, comm_);
00513
00514
00515 Epetra_Vector q(TargetMap);
00516 q.Export(q_dist, Epetra_Export(q_dist.Map(), TargetMap), Add);
00517
00518
00519 if (comm_.MyPID() == 0) {
00520 q_factors.push_back(femax_mesh_->get_nedelec_mesh().q_factor(lambda[k], q.Values()));
00521 }
00522 }
00523 pbe_stop(133);
00524
00525
00526
00527
00528 if (comm_.MyPID() == 0) {
00529
00530 const double pi = ::acos(-1.0);
00531
00532 const double cLight = 2.99792458e8;
00533
00534 ostringstream buf;
00535 buf << "Eigenvalue summary:" << endl;
00536 buf << "nr. norm. frequency freq. [MHz] Q";
00537 for (int k = 0; k < nof_converged_; ++ k) {
00538 double omega = ::sqrt(lambda[k]) * cLight;
00539 double f = omega / 2.0 / pi;
00540 buf << endl << right
00541 << setw(3) << k+1 << " "
00542 << setw(17) << setprecision(8) << fixed << lambda[k] << " "
00543 << setw(17) << setprecision(8) << fixed << f / 1e6 << " "
00544
00545 << setw(12) << setprecision(3) << fixed << q_factors[k];
00546 }
00547 rInfo(buf.str().c_str());
00548 }
00549
00550
00551
00552
00553 if (params_.get<string>("vtk_file_name") != "")
00554 {
00555 pbe_start(134, "VTK export");
00556
00557 rInfo("Export eigensolutions to VTK file...");
00558
00559 postprocess::VtkExport exporter(femax_mesh_->get_nedelec_mesh(),
00560 lambda,
00561 Q);
00562
00563 exporter.export_eigenfields(params_.get<string>("vtk_file_name"));
00564 #if 0
00565 int nof_cells = 1000;
00566 exporter.export_eigenfields2("test_rectangular.vtk", nof_cells);
00567 #endif
00568
00569 #if 0
00570 mesh::Vector3 a(1.0, 1.0, 0.0);
00571 mesh::Vector3 b(1.0, 1.0, 0.77);
00572 exporter.export_a_to_b("test_a_b.vtk", a, b, 100);
00573 #endif
00574 pbe_stop(134);
00575 }
00576
00577
00578
00579
00580 ostringstream buf;
00581 buf << "Parameter list:" << endl;
00582 buf << " [unused] flag is not shown properly for ML paramter lists :-(" << endl;
00583 params_.print(buf, 8);
00584 rDebug(buf.str().c_str());
00585 }
00586
00587 void FemaxxDriver::run()
00588 {
00590 load_mesh();
00591
00593 rInfo("starting calculation of eigenmodal solutions");
00594 calculate_eigenfields();
00595
00597 rInfo("preparing storage of computed eigenmodal solution");
00598 h5_create_empty_file(params_.get<string>("eigensolutiondata_filename").c_str());
00599
00601 rInfo("writing eigenmesh into HDF5 output file = %s", params_.get<string>("eigensolutiondata_filename").c_str());
00602 h5_write_eigenmesh(params_.get<string>("eigensolutiondata_filename").c_str(), comm_, tmesh_);
00603
00605 rInfo("writing mesh nodal coordinates into HDF5 output file");
00606 h5_write_eigenpoint(params_.get<string>("eigensolutiondata_filename").c_str(), comm_, tmesh_);
00607
00608
00610 rInfo("writing eigenvalues into HDF5 file [master process only]");
00611 h5_write_eigenvalue(params_.get<string>("eigensolutiondata_filename").c_str(), comm_, tmesh_, femax_mesh_->get_nedelec_mesh(),*Q_,lambda_);
00612 rInfo("wrote eigenvalues into HDF5 file [master process only]");
00613
00615 rInfo("writing eigenmodal fields sampled at centroids into HDF5 output file");
00616 h5_write_eigenfield(params_.get<string>("eigensolutiondata_filename").c_str(), comm_, tmesh_, femax_mesh_->get_nedelec_mesh(),*Q_,lambda_);
00617
00618
00620 rInfo("sampling computed modal solution on a cartesian grid");
00621 rInfo("this is contribution of processor %d", comm_.MyPID());
00622
00624 double xmin=params_.get<double>("xaxis_cartesian_sampling_start_location");
00625 double xmax=params_.get<double>("xaxis_cartesian_sampling_stop_location");
00626 int nx=params_.get<int>("number_samples_x_axis");
00627 rInfo("sampling interval on x-axis = [%12.10e, %12.10e] & number of samples on x-axis = %d",xmin, xmax, nx);
00628
00630 double ymin=params_.get<double>("yaxis_cartesian_sampling_start_location");
00631 double ymax=params_.get<double>("yaxis_cartesian_sampling_stop_location");
00632 int ny=params_.get<int>("number_samples_y_axis");
00633 rInfo("sampling interval on y-axis = [%12.10e, %12.10e] & number of samples on y-axis = %d",ymin, ymax, ny);
00634
00636 double zmin=params_.get<double>("zaxis_cartesian_sampling_start_location");
00637 double zmax=params_.get<double>("zaxis_cartesian_sampling_stop_location");
00638 int nz=params_.get<int>("number_samples_z_axis");
00639 rInfo("sampling interval on z-axis = [%12.10e, %12.10e] & number of samples on z-axis = %d",zmin, zmax, nz);
00640
00641 rInfo("total number of samples[nx*ny*nz] = %d",nx*ny*nz);
00642
00645 h5_cartesian_sampling(params_.get<string>("eigensolutiondata_filename").c_str(), comm_,
00646 tmesh_, femax_mesh_->get_nedelec_mesh(),
00647 *Q_, lambda_,
00648 xmin, xmax, nx,
00649 ymin, ymax, ny,
00650 zmin, zmax, nz);
00651
00652 }
00653
00657 const Epetra_MultiVector FemaxxDriver::get_eigenvectors() {
00658 return Epetra_MultiVector(View, *Q_, 0, nof_converged_);
00659 }
00660
00664 Epetra_SerialDenseVector FemaxxDriver::get_eigenvalues() {
00665 return Epetra_SerialDenseVector(View, lambda_, nof_converged_);
00666 }