FEMAXX (Finite Element Maxwell Eigensolver): examples/femaxx

00001 
00050 #include <sys/types.h>
00051 #include <sys/stat.h>
00052 #include <fcntl.h>
00053 
00054 #include "rlog/rlog.h"
00055 #include "rlog/StdioNode.h"
00056 #include "rlog/RLogChannel.h"
00057 #include "myrlog.h"
00058 
00059 //#include "Epetra_ConfigDefs.h"
00060 #include <ml_include.h>
00061 
00062 #ifdef HAVE_MPI
00063 #include <mpi.h>
00064 #include <Epetra_MpiComm.h>
00065 #else
00066 #include <Epetra_SerialComm.h>
00067 #endif
00068 
00069 // include after mpi.h
00070 #include "pbe.h"
00071 
00072 #include <Teuchos_ParameterList.hpp>
00073 #include <ml_epetra_preconditioner.h>
00074 
00075 #include "optparse/TypedOptionParser.h"
00076 #include "femaxxdriver.h"
00077 
00078 using namespace rlog;
00079 
00083 void finalize() {
00084 #ifdef HAVE_MPI
00085     MPI_Finalize();
00086 #endif
00087 }
00088 
00097 static int parse_cmd_args(int argc, char *argv[], const Epetra_Comm& comm, Teuchos::ParameterList& params) {
00098     using optparse::STORE;
00099     using optparse::STORE_TRUE;
00100     using optparse::STORE_FALSE;
00101     using optparse::STRING;
00102     using optparse::BOOL;
00103     using optparse::INT;
00104     using optparse::DOUBLE;
00105 
00106     // Initialise parameter list with driver defaults.
00107     FemaxxDriver::set_defaults(params);
00108 
00109     // Setup option parser using defaults from parameter list    
00110     string usage("femaxx_driver\n\nDriver for parallel eigenmode solver femaXX\n\nUsage: femaxx_driver [options]\n\noptions:");
00111     optparse::TypedOptionParser parser(usage);
00112     parser.add_option("-h", "--help", "help", "Show this message", STORE_TRUE, BOOL, "0");
00113     parser.add_option("-m", "--mesh", "mesh_file_name", "Mesh file in HDF5 format", STORE);
00114     parser.add_option("", "--material", "material_file_name", "Material file in regular format", STORE, STRING);
00115     parser.add_option("-s", "--sigma", "sigma", "Shift for preconditioner of A - sigma*M", STORE, DOUBLE);
00116     parser.add_option(params, "-n", "--no-logfile", "disable_log", "Disable logfile", STORE_TRUE);
00117     parser.add_option(params, "-o", "--order", "element_order", "Finite element order", STORE, "1,2");
00118     parser.add_option(params, "-k", "--kmax", "kmax", "Number of eigensolutions to be computed", STORE);
00119     parser.add_option(params, "", "--tol", "eigtol", "Error tolerance for eigenpairs", STORE);
00120     parser.add_option(params, "", "--export-mtx", "export_mtx", "Export matrices in MatrixMarket format", STORE_TRUE);
00121     parser.add_option(params, "-v", "--vtk", "vtk_file_name", "Export eigenfields to VTK file", STORE);
00122     parser.add_option(params, "", "--sym-plane-config", "sym_plane_config", "Boundary conditions on symmetry planes", STORE);
00123     parser.add_option(params, "-d", "--debug", "debug", "Enable debugging code", STORE_TRUE);
00124     parser.add_option(params, "", "--eigsolver", "eigsolver", "Eigensolver", STORE, "jdsym,lobpcg,knyazev");
00125     parser.add_option(params, "", "--aprec", "asigma_precon.type", "Preconditioner type for A - sigma*M", STORE, "neumann,ml,lu,2level,diag,if");
00126     parser.add_option(params, "", "--a11solver", "asigma_precon.solver11", "Solver for (1,1)-block of A - sigma*M", STORE, "lu,ml,if");
00127     parser.add_option(params, "", "--a22solver", "asigma_precon.solver22", "Solver for (2,2)-block of A - sigma*M", STORE, "diag,gmres_incomplete_block_jacobi,pcg_jacobi,blockdiag,none");
00128     parser.add_option(params, "", "--hsolver", "h_solver.type", "Solver type for H", STORE, "lu,pcg,prec");
00129     parser.add_option(params, "", "--hprec", "h_precon.type", "Preconditioner type for H", STORE, "ml,2level,diag,if");
00130     parser.add_option(params, "", "--h11solver", "h_precon.solver11", "Solver for (1,1)-block of H", STORE, "lu,ml,pcg_ml,if");
00131     parser.add_option(params, "", "--h22solver", "h_precon.solver22", "Solver for (2,2)-block of H", STORE, "diag,gmres_incomplete_block_jacobi,pcg_jacobi,blockdiag,none");
00132     parser.add_option(params, "", "--jtau", "jd_params.tau", "Target value", STORE);
00133     parser.add_option(params, "", "--jmin", "jd_params.jmin", "Min. size of search space", STORE);
00134     parser.add_option(params, "", "--jmax", "jd_params.jmax", "Max. size of search space", STORE);
00135     parser.add_option(params, "", "--jitmax", "jd_params.max_it", "Max. #JD iterations", STORE);
00136     parser.add_option(params, "", "--jlinitmax", "jd_params.max_inner_it", "Max. #of inner iterations per solve", STORE);
00137     parser.add_option(params, "", "--jtoldecay", "jd_params.tol_decay", "Factor by which inner error tol is decreasing", STORE);
00138     parser.add_option(params, "", "--eps_tr", "jd_params.eps_tr", "Threshold for shift in correction equation", STORE);
00139     parser.add_option("-eigdf", "--eigendatafile", "eigensolutiondata_filename","name of file containing computed HDF5 formatted eigendata", STORE, STRING);
00140 
00144     parser.add_option("-csxstart", "--cartesian-sampling-x-start", "xaxis_cartesian_sampling_start_location", "x-axis start location for cartesian field sampling",STORE,DOUBLE);
00145     parser.add_option("-csxstop",  "--cartesian-sampling-x-stop",  "xaxis_cartesian_sampling_stop_location",  "x-axis stop location for cartesian field sampling", STORE,DOUBLE);
00146     parser.add_option("-nosx",     "--number-samples-x-axis",      "number_samples_x_axis",                   "number of samples in x-axis direction", STORE, INT); 
00147 
00149     parser.add_option("-csystart", "--cartesian-sampling-y-start", "yaxis_cartesian_sampling_start_location", "y-axis start location for cartesian field sampling",STORE,DOUBLE);
00150     parser.add_option("-csystop",  "--cartesian-sampling-y-stop",  "yaxis_cartesian_sampling_stop_location",  "y-axis stop location for cartesian field sampling", STORE,DOUBLE);
00151     parser.add_option("-nosy",     "--number-samples-y-axis",      "number_samples_y_axis",                   "number of samples in y-axis direction", STORE, INT); 
00152 
00154     parser.add_option("-cszstart", "--cartesian-sampling-z-start", "zaxis_cartesian_sampling_start_location", "z-axis start location for cartesian field sampling",STORE,DOUBLE);
00155     parser.add_option("-cszstop",  "--cartesian-sampling-z-stop",  "zaxis_cartesian_sampling_stop_location",  "z-axis stop location for cartesian field sampling", STORE,DOUBLE);
00156     parser.add_option("-nosz",     "--number-samples-z-axis",      "number_samples_z_axis",                   "number of samples in z-axis direction", STORE, INT); 
00157 
00158 
00159     try {
00160         parser.parse_args(argc, argv);
00161            
00162         // Test if all required options have been set
00163         if (parser.get_option("mesh_file_name") == "")
00164             throw optparse::OptionError("mesh", "Missing required option");
00165         if (parser.get_option("sigma") == "")
00166             throw optparse::OptionError("sigma", "Missing required option");
00167                 if (parser.get_option("eigensolutiondata_filename") == "")
00168                         throw optparse::OptionError("eigensolutiondata_filename", "Missing required option");
00169 
00170         if (parser.get_option("help") == "1") {
00171             if (comm.MyPID() == 0)
00172                 parser.help(cerr);
00173             return 1;
00174         }
00175     }
00176     catch(optparse::OptionError e) {
00177         if (comm.MyPID() == 0) {
00178             cerr << e.what() << endl << endl;
00179             parser.help(cerr);
00180         }
00181         return 1;
00182     }
00183     
00184     if (parser.arguments.size() != 0) {
00185         if (comm.MyPID() == 0) {
00186             cerr << "Unexpected argument \"" << parser.arguments[0] << "\" found." << endl << endl;
00187             parser.help(cerr);
00188         }
00189         return 1;
00190     }
00191     
00192     // Update parameter list with values from parsed options.
00193     parser.set_parameter_list(params);
00194     
00195     return 0;
00196 }
00197 
00198 
00201 int main(int argc, char *argv[])
00202 {
00203 
00204 #ifdef HAVE_MPI
00205     MPI_Init(&argc, &argv);
00206     atexit(finalize);
00207     Epetra_MpiComm Comm(MPI_COMM_WORLD);
00208 #else
00209     Epetra_SerialComm Comm;
00210 #endif
00211 
00212 #ifdef __GNUC__
00213     std::set_terminate (__gnu_cxx::__verbose_terminate_handler);
00214 #endif
00215 
00216     /* FIXME: On jupiter (IBM SP4 @ CSCS) the program crashes (core dump) 
00217        somewhere in the constructor of Teuchos::ParameterList. The cause is 
00218        unknown. */
00219     Teuchos::ParameterList params;
00220     // Parse command line using optparse and initialise parameter list
00221     if (parse_cmd_args(argc, argv, Comm, params) != 0)
00222         // optparse discovered errors: exit program
00223         return 1;
00224     
00225     // Initialise Rlog
00226     RLogInit(argc, argv);
00227 #ifdef HAVE_MPI
00228     char prefix_buf[16];
00229     sprintf(prefix_buf, "p%02d", Comm.MyPID());
00230     MyStdioNode stdLog(prefix_buf);
00231 #else
00232     MyStdioNode stdLog;
00233 #endif
00234     StdioNode* dbgLog(0);
00235 
00236     DEF_CHANNEL("warning/all", Log_Warning);
00237     DEF_CHANNEL("error/all", Log_Error);
00238     DEF_CHANNEL("info/all", Log_Info);
00239     DEF_CHANNEL("debug/all", Log_Debug);
00240 
00241     // stderr output
00242     if (Comm.MyPID() == 0) {
00243         stdLog.subscribeTo(RLOG_CHANNEL("warning"));
00244         stdLog.subscribeTo(RLOG_CHANNEL("error"));
00245         stdLog.subscribeTo(RLOG_CHANNEL("info"));
00246     } else {
00247         stdLog.subscribeTo(RLOG_CHANNEL("error/all"));
00248         stdLog.subscribeTo(RLOG_CHANNEL("warning/all"));
00249         stdLog.subscribeTo(RLOG_CHANNEL("info/all"));
00250     }
00251 
00252     // logfile output
00253     string logname(get_log_file_name("femaxx_"));
00254     if (!params.get<bool>("disable_log") && Comm.MyPID() == 0) {
00255         int fid = ::open(logname.c_str(), O_WRONLY|O_CREAT|O_TRUNC, 0666);
00256         rAssert(fid != -1);
00257         dbgLog = new StdioNode(fid);
00258         dbgLog->subscribeTo(RLOG_CHANNEL("debug"));
00259         dbgLog->subscribeTo(RLOG_CHANNEL("info"));
00260         dbgLog->subscribeTo(RLOG_CHANNEL("warning"));
00261         dbgLog->subscribeTo(RLOG_CHANNEL("error"));
00262     }
00263 
00264     rInfo("Femaxx running on %d processor(s)...", Comm.NumProc());
00265     rDebugAll("Hi from processor %d.", Comm.MyPID());
00266     if (!params.get<bool>("disable_log"))
00267         rInfo("Logging debug output to %s", logname.c_str());
00268 
00269     // Log command-line options and parameter list
00270     log_command_line_args(argc, argv);
00271     if (Comm.MyPID() == 0) {
00272         ostringstream buf;
00273         buf << "Input parameters:" << endl;
00274         params.print(buf, 8);
00275         rDebug(buf.str().c_str());
00276     }
00277 
00278     //
00279     // Use driver to execute computation
00280     //
00281     FemaxxDriver driver(Comm, params);
00282     driver.run();
00283 
00284     return EXIT_SUCCESS;
00285 }
examples/femaxx_main.cpp
