src/femaxxdriver.cpp

Go to the documentation of this file.
00001 //
00002 // C++ Implementation: femaxxdriver
00003 //
00004 // Description:
00005 //
00006 //
00007 // Author: Roman Geus <geus@maxwell>, (C) 2004
00008 //
00009 // Copyright: See COPYING file that comes with this distribution
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     // Init PBE
00072     //
00073     // Timers used by FemaxxDriver and sub-components:
00074 
00075     // Main program parts:
00076     // 50   Life of FemaxxDriver
00077     // 52   Generation of tetrahedral mesh
00078     // 53   DOF mappings (boundary conditions)
00079     // 54   Eigenvalue solver
00080 
00081     // Matrix construction:
00082     // 60   Matrix assembly A and M
00083     // 61   Matrix construction: Y
00084     // 62   Matrix multiplication: C and C^T
00085     // 63   Matrix multiplication: H
00086 
00087     // Preconditioner/Projector construction:
00088     // 70   Preconditioner construction of A-sigma*M
00089     // 71   Preconditioner construction of H
00090     // 72   Solver construction of H
00091     // 75   Nullspace projector construction
00092 
00093     // Operator application:
00094     // 80   Application of A
00095     // 81   Application of M
00096     // 82   Application of Preconditioner for A-sigma*M
00097     // 83   Application of Preconditioner for H
00098     // 84   Application of Nullspace Projector
00099 
00100     // Hierarchical precinditioner:
00101     // 110  Construction
00102     // 111  Operator
00103     // 112  Apply
00104     // 113  (1,1)-solve
00105     // 114  (2,2)-solve
00106     // 115  off-diagonal multiply
00107 
00108     // Utility routines:
00109     // 130  MatrixMarket export
00110     // 131  Octree construction
00111     // 132  Octree traversal
00112     // 133  Quality factors
00113     // 134  VTK export
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     // Finalize PBE
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     // Defaults for cmd line parser
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     // default values for ML preconditioner of H for smoothed aggregation
00160     ML_Epetra::SetDefaults("SA", params.sublist("h_precon").sublist("ml_params"));
00161     // replace Gauss-Seidel smoother by MLS
00162     params.sublist("h_precon").sublist("ml_params").set("smoother: type", "MLS");
00163 
00164     // Amesos params used if h_solver.type == "lu"
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     // Amesos params used if asigma_precon.type == "lu"
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     // ML params for (1,1)-block of A-sigma*M (entire matrix if element_order == 1)
00175     ML_Epetra::SetDefaults("maxwell", params.sublist("asigma_precon").sublist("ml_params"));
00176 
00177     // ML 4.0
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     // FIXME: The repartitioning feature in ML crashes as of Trilinos 6.0.11 -> commented out.
00182     params.sublist("asigma_precon").sublist("ml_params").set("repartition: enable",true);
00183     //how strict the nodal load-balancing is
00184     params.sublist("asigma_precon").sublist("ml_params").set("repartition: node max min ratio",1.1);
00185     //minimum number of nodal unknowns per processor, you might have to change
00186     params.sublist("asigma_precon").sublist("ml_params").set("repartition: node min per proc",200);
00187     //how strict the edge load-balancing is
00188     params.sublist("asigma_precon").sublist("ml_params").set("repartition: edge max min ratio",1.1);
00189     //minimum number of edge unknowns per processor, you might have to change 
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     //params.sublist("asigma_precon").sublist("ml_params").set("aggregation: type", "MIS");
00194     //params.sublist("asigma_precon").sublist("ml_params").set("coarse: max size", 30);
00195     //params.sublist("asigma_precon").sublist("ml_params").set("aggregation: threshold", 0.0);    
00196 
00197     // JDBSYM params
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     //ifstream istr(params_.get<string>("mesh_file_name").c_str());
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     // Generate mesh data from file input on all processors
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 // #ifdef HAVE_MPI
00221 //     BroadcastMeshBuilder* broadcast = new BroadcastMeshBuilder(&comm_, builder);
00222 //     NetgenReader* reader(0);
00223 //     if (comm_.MyPID() == 0) {
00224 //         reader = new NetgenReader(&istr, broadcast);
00225 //         reader->read();
00226 //     } else
00227 //         broadcast->receiver();
00228 //     delete broadcast;
00229 //     if (comm_.MyPID() == 0)
00230 //         delete reader;
00231     HDF5ParallelReader* reader = new HDF5ParallelReader(meshFileName, builder, dynamic_cast<const Epetra_MpiComm&>(comm_).GetMpiComm());
00232     reader->read();
00233     delete reader;
00234 
00235 // #else
00236 //     NetgenReader* reader = new NetgenReader(&istr, builder);
00237 //     reader->read();
00238 //     delete reader;
00239 // #endif
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     // Load material properties (if any)
00249     tmesh_->load_materials(materialFileName);
00250 
00251     // create FE mesh data structures
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     // create operators
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     // Free mesh data structures
00280     // FIXME: these data structures could freed even earlier (after construction of A, M, Y),
00281     // but it would make the program ugly.
00282     log_mem_footprint(comm_, "before mesh data structure dealloc");
00283 //    delete femax_mesh_; /** Do not delete femax_msh here because it will be used in the evaluation of the solution */
00284     
00289 //    delete tmesh_;
00290     
00291 //    femax_mesh_ = 0;
00292 //    tmesh_ = 0;
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     // Here the preconditioner is actually constructed
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     // Eigenvalue solver
00317     //
00318     int myPid = comm_.MyPID();
00319     int numProc = comm_.NumProc();
00320 
00321     // Initialize the random number generator
00322     srand((unsigned) time(0));
00323 
00324     // Initialize memory counters
00325     initMemCounters();
00326 
00327     double highMem = currentSize();
00328     highMem = (highMem < currentSize()) ? currentSize() : highMem;
00329 
00330     int globalSize = (A->OperatorRangeMap()).NumGlobalElements();
00331 
00332     // Get the first eigenvalue for the mass matrix
00333     int kmax = params_.get<int>("kmax");
00334     kmax = (kmax > globalSize) ? globalSize : kmax;
00335 
00336     // Define arrays for the eigensolver algorithm
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     // Projector construction (including construction of solver for H)
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         // tester.set_operator(Prec);
00357         // tester.dump_stats(cout, "A-sigma*M operator");
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     // Define the GeneralEigenSolver object
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_); //applying projector step on blocked set of initial random vectors.
00376         mySolver = new KnyazevLOBPCG(comm_, A, M, Prec, &YHC, eigtol, max_iter, clvl); 
00377     } else if (params_.get<string>("eigsolver") == "lobpcg") {
00378         // FIXME: workspace mgmt of LOBPCG must be corrected to allow for any kmax
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         //hack: 2 eigenvectors are computed. The three other Qs are used as working space
00385         //hack works only if computed eigenvectors are at least 5.
00386         //number of eigenvectors * 2 + 1 = needed size of parameter Q
00387         //if you want to compute x eigenvectors for this version of LOBPCG
00388         //change kmax=x above accordingly and run script
00389         //mpirun ... -kmax=2*x+1 ...
00390     } else {
00391         rError("Unknown eigensolver \"%s\".", params_.get<string>("eigsolver").c_str());
00392     }
00393 
00394     // Solve the eigenproblem
00395     nof_converged_ = mySolver->solve(kmax, *Q_, lambda_);
00396     pbe_stop(54);
00397 
00398     // Converged eigenpairs
00399     Epetra_MultiVector Q = get_eigenvectors();
00400     Epetra_SerialDenseVector lambda = get_eigenvalues();
00401 
00402     // Output information on simulation
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             // Check the computed eigenmodes
00410             CheckingTools myVerify(comm_);
00411             // Check orthonormality of eigenvectors
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             // Print out norm of residuals
00416             myVerify.errorEigenResiduals(Q, lambda_, A, M);
00417 #if 0
00418             // Check C^T Q = 0
00419             Epetra_MultiVector CtQ(C->OperatorDomainMap(), nof_converged_);
00420             (const_cast<Epetra_Operator*>(C))->SetUseTranspose(true);
00421             C->Apply(Q, CtQ);
00422 
00423             // FIXME: Calling the Epetra_MultiVector::NormInf() leads to SIGSEGV for
00424             // small problems. This is probably a bug in Trilinos
00425 
00426             // Get the largest M-norm of Y from H
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         // Memory
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     } // if (nof_converged_ < 0)
00462 
00463     // Release all objects
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     // Export to hdf5
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         // tmesh_->get_octree()->export_vtk();
00493     }
00494 
00495     pbe_stop(131);
00496 
00497     //
00498     // Q-factors
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         // Bring k-th eigenvector to processor 0
00505         //
00506         Epetra_MultiVector q_dist(View, Q, k, 1);
00507 
00508         // create a target map, for which all the elements are on proc 0
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         // redistribute vector to proc 0
00515         Epetra_Vector q(TargetMap);
00516         q.Export(q_dist, Epetra_Export(q_dist.Map(), TargetMap), Add);
00517 
00518         // Compute q-factor on processor 0
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     // Print eigenvalue summary
00527     //
00528     if (comm_.MyPID() == 0) {
00529         // pi
00530         const double pi = ::acos(-1.0);
00531         // speed of light [m / s]
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                 //                << setw(17) << setprecision(8) << scientific << 0.0 << " "
00545                 << setw(12) << setprecision(3) << fixed << q_factors[k];
00546         }
00547         rInfo(buf.str().c_str());
00548     }
00549 
00550  
00551         // Export eigensolutions to VTK file
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 //       h5_write_eigenmodes(params_.get<string>("vtk_file_name"), get_eigenvectors(), get_eigenvalues());
00578 
00579          // log unused parameters
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 }

Generated on Fri Oct 26 13:35:12 2007 for FEMAXX (Finite Element Maxwell Eigensolver) by  doxygen 1.4.7