examples/comet_stems.cpp

Go to the documentation of this file.
00001 //
00002 // C++ Implementation: comet_stems
00003 //
00004 // Description: Main program for eigenmode calculations for the COMET cyclotron.
00005 //
00006 // TODO: Voltage normalisation: 80kV
00007 //
00008 // Author: Roman Geus <roman.geus@psi.ch>, (C) 2004, 2005
00009 //
00010 // Copyright: See COPYING file that comes with this distribution
00011 //
00012 #include <sys/types.h>
00013 #include <sys/stat.h>
00014 #include <fcntl.h>
00015 #include <stdlib.h>
00016 
00017 #include "rlog/rlog.h"
00018 #include "rlog/rloglocation.h"
00019 #include "rlog/Error.h"
00020 #include "rlog/RLogChannel.h"
00021 #include "rlog/StdioNode.h"
00022 #include "rlog/RLogTime.h"
00023 #include "myrlog.h"
00024 
00025 #include <algorithm>
00026 #include <iostream>
00027 #include <fstream>
00028 #include <cmath>
00029 
00030 //#include "Epetra_ConfigDefs.h"
00031 #include <ml_include.h>
00032 
00033 #ifdef HAVE_MPI
00034 #include <mpi.h>
00035 #include <Epetra_MpiComm.h>
00036 #else
00037 #include <Epetra_SerialComm.h>
00038 #endif
00039 
00040 #include "optparse/TypedOptionParser.h"
00041 #include "pbe.h"
00042 
00043 #include <Trilinos_Util_CommandLineParser.h>
00044 #include <Teuchos_ParameterList.hpp>
00045 #include <ml_epetra_preconditioner.h>
00046 
00047 #include "adaptsim.h"
00048 #include "spline.h"
00049 #include "romberg.h"
00050 #include "CircleZ.h"
00051 #include "EfieldIntegrand.h"
00052 #include "LoggingFunction.h"
00053 #include "clp_utils.h"
00054 #include "vector.h"
00055 #include "tetmesh.h"
00056 #include "nedelecmesh.h"
00057 #include "femaxmesh.h"
00058 #include "vtkexport.h"
00059 #include "femaxxdriver.h"
00060 
00061 // Disable Romberg integration and use only adaptsim.
00062 #define DISABLE_QROMB 1
00063 
00064 using namespace rlog;
00065 using namespace std;
00066 
00068 void finalize() {
00069 #ifdef HAVE_MPI
00070     MPI_Finalize();
00071 #endif
00072 }
00073 
00082 static int parse_cmd_args(int argc, char *argv[], const Epetra_Comm& comm, Teuchos::ParameterList& params) {
00083     using optparse::STORE;
00084     using optparse::STORE_TRUE;
00085     using optparse::STORE_FALSE;
00086     using optparse::STRING;
00087     using optparse::BOOL;
00088     using optparse::INT;
00089     using optparse::DOUBLE;
00090 
00091     // Initialise parameter list with driver defaults.
00092     FemaxxDriver::set_defaults(params);
00093     
00094     // Add application-specific parameters
00095     params.set("stem_z_threshold", 0.25400);
00096 
00097     // Setup option parser using defaults from parameter list    
00098     string usage("comet_stems\n\nEigenmode solver application for the COMET cyclotron\n\nUsage: comet_stems [options]\n\noptions:");
00099     optparse::TypedOptionParser parser(usage);
00100 
00101     // Application-specific default values
00102     char* home_dir = getenv("HOME");        // Default mesh is expected to be in $HOME/meshes
00103     string mesh_path;
00104     if (home_dir != 0)
00105         mesh_path += string(home_dir) + "/meshes/";
00106     mesh_path += "P_1499_1641_1641_1612.ng";
00107     parser.add_option("-m", "--mesh", "mesh_file_name", "Mesh file in HDF5 format", STORE, STRING, mesh_path);
00108     parser.add_option("", "--sym-plane-config", "sym_plane_config", "Boundary conditions on symmetry planes", STORE, INT, "1");
00109     parser.add_option("-s", "--sigma", "sigma", "Shift for preconditioner of A - sigma*M", STORE, DOUBLE, "1.1");
00110     
00111     // Application-specific command line parameters
00112     parser.add_option("", "--stem_height_1", "stem_height_1", "Hight of stem 1 in meters", STORE, DOUBLE, "0.218854");
00113     parser.add_option("", "--stem_height_2", "stem_height_2", "Hight of stem 2 in meters", STORE, DOUBLE, "0.239586");
00114     parser.add_option("", "--stem_height_3", "stem_height_3", "Hight of stem 3 in meters", STORE, DOUBLE, "0.239586");
00115     parser.add_option("", "--stem_height_4", "stem_height_4", "Hight of stem 4 in meters", STORE, DOUBLE, "0.235352");
00116    
00117     // Options taken over from femaxx_driver
00118     parser.add_option("-h", "--help", "help", "Show this message", STORE_TRUE, BOOL, "0");
00119     //parser.add_option("-m", "--mesh", "mesh_file_name", "Mesh file in NETGEN format", STORE, STRING);
00120     //parser.add_option("-s", "--sigma", "sigma", "Shift for preconditioner of A - sigma*M", STORE, DOUBLE);
00121     parser.add_option(params, "-n", "--no-logfile", "disable_log", "Disable logfile", STORE_TRUE);
00122     parser.add_option(params, "-o", "--order", "element_order", "Finite element order", STORE, "1,2");
00123     parser.add_option(params, "-k", "--kmax", "kmax", "Number of eigensolutions to be computed", STORE);
00124     parser.add_option(params, "", "--tol", "eigtol", "Error tolerance for eigenpairs", STORE);
00125     parser.add_option(params, "", "--export-mtx", "export_mtx", "Export matrices in MatrixMarket format", STORE_TRUE);
00126     parser.add_option(params, "-v", "--vtk", "vtk_file_name", "Export eigenfields to VTK file", STORE);
00127     //parser.add_option(params, "", "--sym-plane-config", "sym_plane_config", "Boundary conditions on symmetry planes", STORE);
00128     parser.add_option(params, "-d", "--debug", "debug", "Enable debugging code", STORE_TRUE);
00129     parser.add_option(params, "", "--eigsolver", "eigsolver", "Eigensolver", STORE, "jdsym,lobpcg,knyazev");
00130     parser.add_option(params, "", "--aprec", "asigma_precon.type", "Preconditioner type for A - sigma*M", STORE, "neumann,ml,lu,2level,diag,if");
00131     parser.add_option(params, "", "--a11solver", "asigma_precon.solver11", "Solver for (1,1)-block of A - sigma*M", STORE, "lu,ml,if");
00132     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");
00133     parser.add_option(params, "", "--hsolver", "h_solver.type", "Solver type for H", STORE, "lu,pcg,prec");
00134     parser.add_option(params, "", "--hprec", "h_precon.type", "Preconditioner type for H", STORE, "ml,2level,diag,if");
00135     parser.add_option(params, "", "--h11solver", "h_precon.solver11", "Solver for (1,1)-block of H", STORE, "lu,ml,pcg_ml,if");
00136     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");
00137     parser.add_option(params, "", "--jtau", "jd_params.tau", "Target value", STORE);
00138     parser.add_option(params, "", "--jmin", "jd_params.jmin", "Min. size of search space", STORE);
00139     parser.add_option(params, "", "--jmax", "jd_params.jmax", "Max. size of search space", STORE);
00140     parser.add_option(params, "", "--jitmax", "jd_params.max_it", "Max. #JD iterations", STORE);
00141     parser.add_option(params, "", "--jlinitmax", "jd_params.max_inner_it", "Max. #of inner iterations per solve", STORE);
00142     parser.add_option(params, "", "--jtoldecay", "jd_params.tol_decay", "Factor by which inner error tol is decreasing", STORE);
00143     parser.add_option(params, "", "--eps_tr", "jd_params.eps_tr", "Threshold for shift in correction equation", STORE);
00144     try {
00145         parser.parse_args(argc, argv);
00146            
00147         // Test if all required options have been set
00148         if (parser.get_option("mesh_file_name") == "")
00149             throw optparse::OptionError("mesh", "Missing required option");
00150         if (parser.get_option("sigma") == "")
00151             throw optparse::OptionError("sigma", "Missing required option");
00152     
00153         if (parser.get_option("help") == "1") {
00154             if (comm.MyPID() == 0)
00155                 parser.help(cerr);
00156             return 1;
00157         }
00158     }
00159     catch(optparse::OptionError e) {
00160         if (comm.MyPID() == 0) {
00161             cerr << e.what() << endl << endl;
00162             parser.help(cerr);
00163         }
00164         return 1;
00165     }
00166     
00167     if (parser.arguments.size() != 0) {
00168         if (comm.MyPID() == 0) {
00169             cerr << "Unexpected argument \"" << parser.arguments[0] << "\" found." << endl << endl;
00170             parser.help(cerr);
00171         }
00172         return 1;
00173     }
00174     
00175     // Update parameter list with values from parsed options.
00176     parser.set_parameter_list(params);
00177     
00178     return 0;
00179 }
00180 
00182 class CometGapFinder {
00183 public:
00184     CometGapFinder(std::string mesh_type, const mesh::TetMesh& tmesh)
00185         : tmesh_(tmesh)
00186     {
00187         if (mesh_type == "P_1499_1641_1641_1612") {
00188             // Tabulated functions (used for spline interpolation
00189             double rad_entry_in[] = {0.10243483, 0.24759676, 0.37891511, 0.48161345, 0.58906747, 0.67218756, 0.74142964, 0.82499992};
00190             double phi_entry_in[] = {4.06171642, 4.37300720, 4.62760356, 4.77981775, 4.94923057, 5.10438120, 5.27123020, 5.32050234};
00191             double rad_entry_out[] = {0.08249764, 0.15229278, 0.39877397, 0.60499174, 0.81999979};
00192             double phi_entry_out[] = {4.59501326, 4.67390096, 4.91424909, 5.31862873, 5.59166061};
00193             double rad_exit_in[] = {0.13510732, 0.18577222, 0.24290023, 0.27314372, 0.33987641, 0.41934085, 0.46370478, 0.53079440, 0.61135473, 0.68823301, 0.73127536, 0.75081509, 0.76798017};
00194             double phi_exit_in[] = {3.74278291, 3.58438653, 3.62984776, 3.71049239, 3.76697089, 3.87803507, 3.95716874, 4.06445940, 4.22136555, 4.41000000, 4.52095706, 4.58974001, 4.66243071};
00195             double rad_exit_out[] = {0.06894548, 0.15229278, 0.27732126, 0.40163350, 0.60412626, 0.70856021, 0.81999968};
00196             double phi_exit_out[] = {2.96489808, 3.10310464, 3.22001175, 3.46130553, 3.77422794, 4.08845302, 4.24000150};
00197 
00198             // Copy arrays into vectors (insert rather than assign)
00199             std::copy(rad_entry_in, rad_entry_in + sizeof(rad_entry_in)/sizeof(rad_entry_in[0]),
00200                       std::inserter(rad_entry_in_, rad_entry_in_.begin()));
00201             std::copy(phi_entry_in, phi_entry_in + sizeof(phi_entry_in)/sizeof(phi_entry_in[0]),
00202                       std::inserter(phi_entry_in_, phi_entry_in_.begin()));
00203             std::copy(rad_entry_out, rad_entry_out + sizeof(rad_entry_out)/sizeof(rad_entry_out[0]),
00204                       insert_iterator< vector<double> > (rad_entry_out_, rad_entry_out_.begin()));
00205             std::copy(phi_entry_out, phi_entry_out + sizeof(phi_entry_out)/sizeof(phi_entry_out[0]),
00206                       insert_iterator< vector<double> > (phi_entry_out_, phi_entry_out_.begin()));
00207             std::copy(rad_exit_in, rad_exit_in + sizeof(rad_exit_in)/sizeof(rad_exit_in[0]),
00208                       insert_iterator< vector<double> > (rad_exit_in_, rad_exit_in_.begin()));
00209             std::copy(phi_exit_in, phi_exit_in + sizeof(phi_exit_in)/sizeof(phi_exit_in[0]),
00210                       insert_iterator< vector<double> > (phi_exit_in_, phi_exit_in_.begin()));
00211             std::copy(rad_exit_out, rad_exit_out + sizeof(rad_exit_out)/sizeof(rad_exit_out[0]),
00212                       insert_iterator< vector<double> > (rad_exit_out_, rad_exit_out_.begin()));
00213             std::copy(phi_exit_out, phi_exit_out + sizeof(phi_exit_out)/sizeof(phi_exit_out[0]),
00214                       insert_iterator< vector<double> > (phi_exit_out_, phi_exit_out_.begin()));
00215 
00216             // Spline preprocessing (compute coefficients)
00217             NR::spline(rad_entry_in_, phi_entry_in_, 1e40, 1e40, phi_entry_in2_);
00218             NR::spline(rad_entry_out_, phi_entry_out_, 1e40, 1e40, phi_entry_out2_);
00219             NR::spline(rad_exit_in_, phi_exit_in_, 1e40, 1e40, phi_exit_in2_);
00220             NR::spline(rad_exit_out_, phi_exit_out_, 1e40, 1e40, phi_exit_out2_);
00221         } else {
00222             // Mesh type not supported
00223             rExit("Mesh not supported by CometGapFinder.");
00224         }
00225 
00226         // Compute min and max radius
00227         min_radius_ = rad_entry_in_[0];
00228         min_radius_ = std::max<double>(min_radius_, rad_entry_out_[0]);
00229         min_radius_ = std::max<double>(min_radius_, rad_exit_in_[0]);
00230         min_radius_ = std::max<double>(min_radius_, rad_exit_out_[0]);
00231 
00232         max_radius_ = rad_entry_in_.end()[-1];
00233         max_radius_ = std::min<double>(max_radius_, rad_entry_out_.end()[-1]);
00234         max_radius_ = std::min<double>(max_radius_, rad_exit_in_.end()[-1]);
00235         max_radius_ = std::min<double>(max_radius_, rad_exit_out_.end()[-1]);
00236     }
00237 
00242     double get_min_radius() const {
00243         return min_radius_;
00244     }
00245 
00250     double get_max_radius() const {
00251         return max_radius_;
00252     }
00253 
00265     void get_gap(int gap_id, double radius, double& phi_entry, double& phi_exit) {
00266         const double pi = ::acos(-1.0);
00267         const double z_gap_find = 0.2;                 // Height at which the wing walls are detected
00268         const double phi_range = 0.6;                  // Gap width. Angle across which the voltage is computed
00269         double phase[4] = {0.0, 0.5*pi, pi, 1.5*pi};   // Phase of the four wings wrt to wing 0
00270         int wing_id = gap_id/2;
00271         
00272         assert(0 <= gap_id && gap_id < 8);
00273         assert(min_radius_ <= radius && radius <= max_radius_);
00274 
00275         // Find start values for bisection using spline functions
00276         double phi_entry_in, phi_entry_out, phi_exit_in, phi_exit_out;
00277         NR::splint(rad_entry_out_, phi_entry_out_, phi_entry_out2_, radius, phi_entry_out);
00278         NR::splint(rad_entry_in_, phi_entry_in_, phi_entry_in2_, radius, phi_entry_in);
00279         NR::splint(rad_exit_out_, phi_exit_out_, phi_exit_out2_, radius, phi_exit_out);
00280         NR::splint(rad_exit_in_, phi_exit_in_, phi_exit_in2_, radius, phi_exit_in);
00281 
00282         // Find mesh boundaries using bisection
00283         CircleZ circle_f(radius, z_gap_find);
00284 
00285         rDebug("radius=%lf, phi_entry_out=%lf [%s], phi_entry_in=%lf [%s], phi_exit_out=%lf [%s], phi_exit_in=%lf [%s]",
00286                radius,
00287                phi_entry_out, tmesh_.is_inside(circle_f(phi_entry_out)) ? "T" : "F",
00288                phi_entry_in, tmesh_.is_inside(circle_f(phi_entry_in)) ? "T" : "F",
00289                phi_exit_out, tmesh_.is_inside(circle_f(phi_exit_out)) ? "T" : "F",
00290                phi_exit_in, tmesh_.is_inside(circle_f(phi_exit_in)) ? "T" : "F");
00291 
00292         double phi_begin = tmesh_.find_boundary<CircleZ>(circle_f, phi_exit_in, phi_exit_out);
00293         double phi_end = tmesh_.find_boundary<CircleZ>(circle_f, phi_entry_in, phi_entry_out);
00294         rDebug("wing_id=%d, radius=%lf, phi_begin=%lf, phi_end=%lf", wing_id, radius, phi_begin, phi_end);
00295 
00296         phi_begin += phase[wing_id];
00297         phi_end += phase[wing_id];
00298         if (gap_id % 2 == 0) {
00299             // gap at beginning of the wing
00300             phi_entry = phi_begin - 0.5*phi_range;
00301             phi_exit = phi_begin + 0.5*phi_range;
00302         } else {
00303             // gap at end of wing
00304             phi_entry = phi_end - 0.5*phi_range;
00305             phi_exit = phi_end + 0.5*phi_range;
00306         }            
00307     }
00308 
00309 private:
00310     const mesh::TetMesh& tmesh_;
00311     double min_radius_;
00312     double max_radius_;
00313     std::vector<double> rad_entry_in_, rad_entry_out_, phi_entry_in_, phi_entry_out_, phi_entry_in2_, phi_entry_out2_;
00314     std::vector<double> rad_exit_in_, rad_exit_out_, phi_exit_in_, phi_exit_out_, phi_exit_in2_, phi_exit_out2_;
00315 };
00316 
00323 void eval_tangetial_efield(FemaxxDriver& driver, const colarray::Vector<double>& q, string gnuplot_filename) {
00324     const double pi = ::acos(-1.0);
00325     const int nof_radii = 10;
00326     const int nof_angles = 1440;
00327     ofstream of(gnuplot_filename.c_str());
00328     mesh::TetMesh& tmesh(driver.get_tet_mesh());
00329     NedelecMesh& nedelec_mesh = driver.get_femax_mesh().get_nedelec_mesh();
00330     
00331     rAssert(tmesh.get_octree() != 0);
00332 
00333     CometGapFinder gap_finder("P_1499_1641_1641_1612", tmesh);
00334     const double radius_min = gap_finder.get_min_radius();
00335     const double radius_max = gap_finder.get_max_radius();
00336 
00337     for (int j = 0; j < nof_angles; ++ j) {
00338         double phi = 2.0*pi*j/nof_angles;
00339         of << phi;
00340         
00341         for (int i = 0; i <= nof_radii; ++ i) {
00342             double radius = radius_min + i * (radius_max - radius_min) / nof_radii;
00343             CircleZ circle_f(radius, 0.0);
00344             EfieldDs<CircleZ> integrand(circle_f, nedelec_mesh, q);
00345                     
00346             double e_tangential = integrand(phi);
00347             of << " " << e_tangential;
00348         }
00349         of << endl;
00350     }
00351     of.close();
00352 }
00353 
00354 
00361 void eval_gap_voltage(FemaxxDriver& driver, const colarray::Vector<double>& q, string gnuplot_filename) {
00362     const double tol = 1e-10;
00363     int info;
00364     const int nof_radii = 20;
00365     mesh::TetMesh& tmesh(driver.get_tet_mesh());
00366     NedelecMesh& nedelec_mesh = driver.get_femax_mesh().get_nedelec_mesh();
00367     ofstream of(gnuplot_filename.c_str());
00368     
00369     rAssert(tmesh.get_octree() != 0);
00370 
00371     CometGapFinder gap_finder("P_1499_1641_1641_1612", tmesh);
00372     const double radius_min = gap_finder.get_min_radius();
00373     const double radius_max = gap_finder.get_max_radius();
00374 
00375     rDebug("radius_min=%lf, radius_max=%lf, steps=%d", radius_min, radius_max, nof_radii);
00376 
00377     for (int i = 0; i <= nof_radii; ++ i) {
00378         double radius = radius_min + i * (radius_max - radius_min) / nof_radii;
00379         double phi_entry, phi_exit;
00380         CircleZ circle_f(radius, 0.0);
00381         EfieldDs<CircleZ> integrand(circle_f, nedelec_mesh, q);
00382         LoggingFunction<EfieldDs<CircleZ>, double> log_integrand(integrand);
00383         
00384         of << radius;
00385 
00386         for (int gap_id = 0; gap_id < 8; ++ gap_id) {
00387             // Get entry and exit angle for gap gap_id
00388             gap_finder.get_gap(gap_id, radius, phi_entry, phi_exit);
00389             
00390             // Compute voltage using romberg
00391 #ifndef DISABLE_QROMB
00392             double gap_voltage_qromb = NR::qromb(integrand, phi_entry, phi_exit, tol, info);
00393             rInfo("gap %d: radius=%f, phi=%f->%f, voltage=%e", gap_id, radius, phi_entry, phi_exit, gap_voltage_qromb);
00394 #endif
00395             // Open new log file (for use with gnuplot)    
00396             if (driver.get_comm().MyPID() == 0) {
00397                 stringstream os;
00398                 os << "comet_stems-" << i << "-" << gap_id << ".tab";
00399                 log_integrand.enable_plot(os.str());
00400             }
00401             
00402             // Compute voltage using adaptsim
00403             double gap_voltage_adaptsim = NR::adaptsim(log_integrand, phi_entry, phi_exit, tol, info);
00404             rInfo("gap %d: radius=%f, phi=%f->%f, voltage=%e", gap_id, radius, phi_entry, phi_exit, gap_voltage_adaptsim);
00405             rInfo("num_eval=%d, y_min=%e, y_max=%e", log_integrand.get_num_eval(), log_integrand.get_y_min(), log_integrand.get_y_max());
00406 
00407 #ifndef DISABLE_QROMB
00408             if (fabs(gap_voltage_adaptsim - gap_voltage_qromb) > 0.001 * fabs(gap_voltage_adaptsim)) {
00409                 rWarning("Voltages computed by qromb and adaptsim differ by more than 1/1000.");
00410             }
00411 #endif
00412             
00413             // Reset statistics and close logfile
00414             log_integrand.reset();
00415             
00416             // write voltage to output file
00417             of << " " << gap_voltage_adaptsim;
00418         } // for gap_id
00419         of << endl;
00420     } // for i
00421     
00422     of.close();
00423 }
00424 
00426 void export_mesh(string file, const mesh::TetMesh* tmesh) {
00427     int nofPoints = tmesh->get_nof_points();
00428     int nofTets = tmesh->get_nof_tets();
00429 
00430     mesh::Vector3 coord;
00431     mesh::Tet* curTet;
00432     std::ofstream of;
00433 
00434     of.open(file.c_str());
00435     assert(of.is_open());
00436 
00437     of.precision(5);
00438     of << "# vtk DataFile Version 2.0" << endl;                      // first line of a vtk file
00439     of << "generated using VtkExport::export_eigenfields" << endl;   // header
00440     of << "ASCII" << endl << endl;                                   // file format
00441     of << "DATASET UNSTRUCTURED_GRID" << endl;
00442     of << "POINTS " << nofPoints << " float" << endl;
00443     of << scientific;
00444     for (int i = 0; i < nofPoints  ; i++) {
00445         coord   = tmesh->get_point(i)->get_coord();
00446         of << coord.x << " " << coord.y << " " << coord.z << endl;
00447     }
00448     of << endl;
00449     of << "CELLS " <<  nofTets << " " << 5*nofTets << endl;
00450 
00451     // for each cell and write IDs of its points (preceded by number of points
00452     // in that cell)
00453     for (int i = 0; i < nofTets; i ++) {
00454         curTet = tmesh->get_tet(i);
00455         of << 4;    //assume that it is always a tetrahedron, i.e. with 4 corners
00456         for (int j = 0; j < 4; j ++) {
00457             of << " " << curTet->get_corner(j);
00458         }
00459         of << endl;
00460     }
00461     of << endl;
00462 
00463     // cell-type always = 10 as it is a tetrahedron
00464     of << "CELL_TYPES " << nofTets << endl;
00465     for (int i = 0; i < nofTets; i++) {
00466         of << 10 << endl;
00467     }
00468     of << endl;
00469 
00470     // write dummy point data (zero vectors)
00471     of << "POINT_DATA " << nofPoints << endl << endl;
00472 
00473     of << "VECTORS dummy" << " float" << endl;
00474     for (int j =0; j < nofPoints; j++) {
00475         of << 0.0 << " " << 0.0 << " " << 0.0 << endl;
00476     }
00477     of << endl;
00478 }
00479 
00482 int main(int argc, char *argv[]) {
00483 
00484 #ifdef HAVE_MPI
00485     MPI_Init(&argc, &argv);
00486     atexit(finalize);
00487     Epetra_MpiComm comm(MPI_COMM_WORLD);
00488 #else
00489     Epetra_SerialComm comm;
00490 #endif
00491 
00492 #ifdef __GNUC__
00493     std::set_terminate (__gnu_cxx::__verbose_terminate_handler);
00494 #endif
00495 
00496     Teuchos::ParameterList params;
00497     // Parse command line using optparse and initialise parameter list
00498     if (parse_cmd_args(argc, argv, comm, params) != 0)
00499         // optparse discovered errors: exit program
00500         return 1;
00501     
00502     RLogInit(argc, argv);
00503 
00504     MyStdioNode stdLog;     // Logs to stderr
00505     StdioNode* dbgLog(0);   // Logs to file
00506 
00507     DEF_CHANNEL("warning/all", Log_Warning);
00508     DEF_CHANNEL("error/all", Log_Error);
00509     DEF_CHANNEL("info/all", Log_Info);
00510     DEF_CHANNEL("debug/all", Log_Debug);
00511 
00512     if (comm.MyPID() == 0) {
00513         stdLog.subscribeTo(RLOG_CHANNEL("warning"));
00514         stdLog.subscribeTo(RLOG_CHANNEL("error"));
00515         stdLog.subscribeTo(RLOG_CHANNEL("info"));
00516     } else {
00517         stdLog.subscribeTo(RLOG_CHANNEL("error"));
00518         stdLog.subscribeTo(RLOG_CHANNEL("warning/all"));
00519         stdLog.subscribeTo(RLOG_CHANNEL("info/all"));
00520     }
00521     string logname(get_log_file_name("comet_stems_"));
00522     if (!params.get<bool>("disable_log") && comm.MyPID() == 0) {
00523         int fid = ::open(logname.c_str(), O_WRONLY|O_CREAT|O_TRUNC, 0666);
00524         rAssert(fid != -1);
00525         dbgLog = new StdioNode(fid);
00526         dbgLog->subscribeTo(RLOG_CHANNEL("debug"));
00527         dbgLog->subscribeTo(RLOG_CHANNEL("info"));
00528         dbgLog->subscribeTo(RLOG_CHANNEL("warning"));
00529         dbgLog->subscribeTo(RLOG_CHANNEL("error"));
00530     }
00531 
00532     rInfo("Comet_stems running on %d processor(s)...", comm.NumProc());
00533     rLog(RLOG_CHANNEL("info/all"), "Hi from processor %d.", comm.MyPID());
00534     if (!params.get<bool>("disable_log"))
00535         rInfo("Logging debug output to %s", logname.c_str());
00536 
00537     // Log command-line options and parameter list
00538     log_command_line_args(argc, argv);
00539     if (comm.MyPID() == 0) {
00540         rDebug("Input parameters:");
00541         ostringstream buf;
00542         buf << endl;
00543         params.print(buf, 8);
00544         rDebug(buf.str().c_str());
00545     }
00546 
00547     FemaxxDriver driver(comm, params);
00548     driver.load_mesh();
00549 
00550     if (params.get<string>("vtk_file_name") != "") {
00551         string file_name = params.get<string>("vtk_file_name") + "_mesh_orig.vtk";
00552         if (comm.MyPID() == 0) {
00553             export_mesh(file_name, &driver.get_tet_mesh());
00554         }
00555     }
00556 
00557     rInfo("Rescale COMET stems...");
00558     // Extract mesh points belonging to stems
00559     const double z_threshold = params.get<double>("stem_z_threshold");
00560 
00561     int stem_count[4] = {0, 0, 0, 0};
00562     double stem_zmax[4] = {0.0, 0.0, 0.0, 0.0};
00563     double stem_factor[4] = {0.0, 0.0, 0.0, 0.0};
00564     double stem_newheight[4] = { params.get<double>("stem_height_1"),
00565                                  params.get<double>("stem_height_2"),
00566                                  params.get<double>("stem_height_3"),
00567                                  params.get<double>("stem_height_4") };
00568 
00569     pbe_start(10, "Stem node adjustment")
00570     mesh::TetMesh& tmesh(driver.get_tet_mesh());
00571     for (mesh::TetMesh::point_iterator it(tmesh.point_begin()); it < tmesh.point_end(); ++ it) {
00572         mesh::Vector3 coord((*it).get_coord());
00573         if (coord.z > z_threshold) {
00574             int stem_id = 0;
00575             if (coord.x > 0.0)
00576                 stem_id += 1;
00577             if (coord.y > 0.0)
00578                 stem_id += 2;
00579             stem_count[stem_id] ++;
00580             if (coord.z > stem_zmax[stem_id])
00581                 stem_zmax[stem_id] = coord.z;
00582         }
00583     }
00584 
00585     if (comm.MyPID() == 0) {
00586         for (int i = 0; i < 4; i ++)
00587             rDebug("Stem %d: nof_points=%d, zmax=%lf", i, stem_count[i], stem_zmax[i]);
00588     }
00589 
00590     // reposition stem mesh points according to input parameters
00591     for (int i = 0; i < 4; i ++)
00592         stem_factor[i] = stem_newheight[i] / (stem_zmax[i] - z_threshold);
00593 
00594     for (mesh::TetMesh::point_iterator it(tmesh.point_begin()); it < tmesh.point_end(); ++ it) {
00595         mesh::Vector3& coord((*it).get_coord());
00596         if (coord.z > z_threshold) {
00597             int stem_id = 0;
00598             if (coord.x > 0.0)
00599                 stem_id += 1;
00600             if (coord.y > 0.0)
00601                 stem_id += 2;
00602             stem_count[stem_id] ++;
00603             if (coord.z > stem_zmax[stem_id])
00604                 stem_zmax[stem_id] = coord.z;
00605             coord.z = (coord.z - z_threshold) * stem_factor[stem_id] + z_threshold;
00606         }
00607     }
00608 
00609     if (comm.MyPID() == 0) {
00610         for (int i = 0; i < 4; i ++)
00611             rInfo("Stem %d: old height=%lf, new height=%lf", i, stem_zmax[i] - z_threshold, stem_newheight[i]);
00612     }
00613     pbe_stop(10);
00614 
00615     if (params.get<string>("vtk_file_name") != "") {
00616         string file_name = params.get<string>("vtk_file_name") + "_mesh_modded.vtk";
00617         if (comm.MyPID() == 0) {
00618             export_mesh(file_name, &driver.get_tet_mesh());
00619         }
00620     }
00621 
00622     // calculate eigenfields
00623     driver.calculate_eigenfields();
00624 
00625     // reload mesh data structures
00626     driver.load_mesh();
00627 
00628     // print eigenfield info and export eigenfields
00629     driver.postprocess();
00630 
00631     // evaluate gap voltages
00632     rInfo("Evaluate gap voltages...");
00633     postprocess::VtkExport post(driver.get_femax_mesh().get_nedelec_mesh(),
00634                                 driver.get_eigenvalues(),
00635                                 driver.get_eigenvectors());
00636     colarray::Vector<double>* q(0);
00637     double lambda;
00638     if (comm.MyPID() == 0)
00639         q = new colarray::Vector<double>(post.get_nof_dofs());
00640     // get 1st eigenpair
00641     post.get_eigenpair(0, &lambda, q);
00642     if (comm.MyPID() == 0) {
00643         eval_tangetial_efield(driver, *q, "tangential.tab");
00644         eval_gap_voltage(driver, *q, "gap_voltages.tab");
00645         delete q;
00646     }
00647 
00648     return EXIT_SUCCESS;
00649 }

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