00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
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
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
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
00092 FemaxxDriver::set_defaults(params);
00093
00094
00095 params.set("stem_z_threshold", 0.25400);
00096
00097
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
00102 char* home_dir = getenv("HOME");
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
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
00118 parser.add_option("-h", "--help", "help", "Show this message", STORE_TRUE, BOOL, "0");
00119
00120
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
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
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
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
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
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
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
00223 rExit("Mesh not supported by CometGapFinder.");
00224 }
00225
00226
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;
00268 const double phi_range = 0.6;
00269 double phase[4] = {0.0, 0.5*pi, pi, 1.5*pi};
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
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
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
00300 phi_entry = phi_begin - 0.5*phi_range;
00301 phi_exit = phi_begin + 0.5*phi_range;
00302 } else {
00303
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
00388 gap_finder.get_gap(gap_id, radius, phi_entry, phi_exit);
00389
00390
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
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
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
00414 log_integrand.reset();
00415
00416
00417 of << " " << gap_voltage_adaptsim;
00418 }
00419 of << endl;
00420 }
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;
00439 of << "generated using VtkExport::export_eigenfields" << endl;
00440 of << "ASCII" << endl << endl;
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
00452
00453 for (int i = 0; i < nofTets; i ++) {
00454 curTet = tmesh->get_tet(i);
00455 of << 4;
00456 for (int j = 0; j < 4; j ++) {
00457 of << " " << curTet->get_corner(j);
00458 }
00459 of << endl;
00460 }
00461 of << endl;
00462
00463
00464 of << "CELL_TYPES " << nofTets << endl;
00465 for (int i = 0; i < nofTets; i++) {
00466 of << 10 << endl;
00467 }
00468 of << endl;
00469
00470
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
00498 if (parse_cmd_args(argc, argv, comm, params) != 0)
00499
00500 return 1;
00501
00502 RLogInit(argc, argv);
00503
00504 MyStdioNode stdLog;
00505 StdioNode* dbgLog(0);
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
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
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
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
00623 driver.calculate_eigenfields();
00624
00625
00626 driver.load_mesh();
00627
00628
00629 driver.postprocess();
00630
00631
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
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 }