From ef962f6008f25ab7cbd4ca21bcc72b97a1e2d76f Mon Sep 17 00:00:00 2001 From: Ruben Undheim Date: Tue, 5 Jul 2016 18:02:38 +0200 Subject: Imported Upstream version 0.0.34 --- openEMS/openems.cpp | 1231 +++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 1231 insertions(+) create mode 100644 openEMS/openems.cpp (limited to 'openEMS/openems.cpp') diff --git a/openEMS/openems.cpp b/openEMS/openems.cpp new file mode 100644 index 0000000..72476d3 --- /dev/null +++ b/openEMS/openems.cpp @@ -0,0 +1,1231 @@ +/* +* Copyright (C) 2010-2015 Thorsten Liebig (Thorsten.Liebig@gmx.de) +* +* This program is free software: you can redistribute it and/or modify +* it under the terms of the GNU General Public License as published by +* the Free Software Foundation, either version 3 of the License, or +* (at your option) any later version. +* +* This program is distributed in the hope that it will be useful, +* but WITHOUT ANY WARRANTY; without even the implied warranty of +* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +* GNU General Public License for more details. +* +* You should have received a copy of the GNU General Public License +* along with this program. If not, see . +*/ + +#include "openems.h" +#include +#include +#include +#include "tools/array_ops.h" +#include "tools/useful.h" +#include "FDTD/operator_cylinder.h" +#include "FDTD/operator_cylindermultigrid.h" +#include "FDTD/engine_multithread.h" +#include "FDTD/operator_multithread.h" +#include "FDTD/extensions/operator_ext_excitation.h" +#include "FDTD/extensions/operator_ext_tfsf.h" +#include "FDTD/extensions/operator_ext_mur_abc.h" +#include "FDTD/extensions/operator_ext_upml.h" +#include "FDTD/extensions/operator_ext_lorentzmaterial.h" +#include "FDTD/extensions/operator_ext_conductingsheet.h" +#include "FDTD/extensions/operator_ext_steadystate.h" +#include "FDTD/extensions/engine_ext_steadystate.h" +#include "FDTD/engine_interface_fdtd.h" +#include "FDTD/engine_interface_cylindrical_fdtd.h" +#include "Common/processvoltage.h" +#include "Common/processcurrent.h" +#include "Common/processfieldprobe.h" +#include "Common/processmodematch.h" +#include "Common/processfields_td.h" +#include "Common/processfields_fd.h" +#include "Common/processfields_sar.h" +#include // only for H5get_libversion() +#include // only for BOOST_LIB_VERSION +#include + +//external libs +#include "tinyxml.h" +#include "ContinuousStructure.h" +#include "CSPropProbeBox.h" +#include "CSPropDumpBox.h" + +using namespace std; + +double CalcDiffTime(timeval t1, timeval t2) +{ + double s_diff = t1.tv_sec - t2.tv_sec; + s_diff += (t1.tv_usec-t2.tv_usec)*1e-6; + return s_diff; +} + +openEMS::openEMS() +{ + FDTD_Op=NULL; + FDTD_Eng=NULL; + Eng_Ext_SSD=NULL; + m_CSX=NULL; + PA=NULL; + CylinderCoords = false; + Enable_Dumps = true; + DebugMat = false; + DebugOp = false; + m_debugCSX = false; + m_debugBox = m_debugPEC = m_no_simulation = false; + m_DumpStats = false; + endCrit = 1e-6; + m_OverSampling = 4; + m_CellConstantMaterial=false; + + m_engine = EngineType_Multithreaded; //default engine type + m_engine_numThreads = 0; + + m_Abort = false; + m_Exc = 0; + + m_TS_method=3; + m_TS=0; + m_TS_fac=1.0; + m_maxTime=0.0; + + for (int n=0;n<6;++n) + { + m_BC_type[n] = 0; + m_PML_size[n] = 8; + m_Mur_v_ph[n] = 0; + } +} + +openEMS::~openEMS() +{ + Reset(); +} + +void openEMS::Reset() +{ + if (PA) PA->DeleteAll(); + delete PA; + PA=0; + delete FDTD_Eng; + FDTD_Eng=0; + delete FDTD_Op; + FDTD_Op=0; + delete m_CSX; + m_CSX=0; + delete m_Exc; + m_Exc=0; +} + +void openEMS::showUsage() +{ + cout << " Usage: openEMS [...]" << endl << endl; + cout << " " << endl; + cout << "\t--disable-dumps\t\tDisable all field dumps for faster simulation" << endl; + cout << "\t--debug-material\tDump material distribution to a vtk file for debugging" << endl; + cout << "\t--debug-PEC\t\tDump metal distribution to a vtk file for debugging" << endl; + cout << "\t--debug-operator\tDump operator to vtk file for debugging" << endl; + cout << "\t--debug-boxes\t\tDump e.g. probe boxes to vtk file for debugging" << endl; + cout << "\t--debug-CSX\t\tWrite CSX geometry file to debugCSX.xml" << endl; + cout << "\t--engine=\t\tChoose engine type" << endl; + cout << "\t\t--engine=fastest\t\tfastest available engine (default)" << endl; + cout << "\t\t--engine=basic\t\t\tbasic FDTD engine" << endl; + cout << "\t\t--engine=sse\t\t\tengine using sse vector extensions" << endl; + cout << "\t\t--engine=sse-compressed\t\tengine using compressed operator + sse vector extensions" << endl; +#ifdef MPI_SUPPORT + cout << "\t\t--engine=MPI\t\t\tengine using compressed operator + sse vector extensions + MPI parallel processing" << endl; + cout << "\t\t--engine=multithreaded\t\tengine using compressed operator + sse vector extensions + MPI + multithreading" << endl; +#else + cout << "\t\t--engine=multithreaded\t\tengine using compressed operator + sse vector extensions + multithreading" << endl; +#endif + cout << "\t--numThreads=\tForce use n threads for multithreaded engine (needs: --engine=multithreaded)" << endl; + cout << "\t--no-simulation\t\tonly run preprocessing; do not simulate" << endl; + cout << "\t--dump-statistics\tdump simulation statistics to '" << __OPENEMS_RUN_STAT_FILE__ << "' and '" << __OPENEMS_STAT_FILE__ << "'" << endl; + cout << "\n\t Additional global arguments " << endl; + g_settings.ShowArguments(cout,"\t"); + cout << endl; +} + +//! \brief processes a command line argument +//! \return true if argument is known +//! \return false if argument is unknown +bool openEMS::parseCommandLineArgument( const char *argv ) +{ + if (!argv) + return false; + + if (strcmp(argv,"--disable-dumps")==0) + { + cout << "openEMS - disabling all field dumps" << endl; + SetEnableDumps(false); + return true; + } + else if (strcmp(argv,"--debug-material")==0) + { + cout << "openEMS - dumping material to 'material_dump.vtk'" << endl; + DebugMaterial(); + return true; + } + else if (strcmp(argv,"--debug-operator")==0) + { + cout << "openEMS - dumping operator to 'operator_dump.vtk'" << endl; + DebugOperator(); + return true; + } + else if (strcmp(argv,"--debug-boxes")==0) + { + cout << "openEMS - dumping boxes to 'box_dump*.vtk'" << endl; + DebugBox(); + return true; + } + else if (strcmp(argv,"--debug-PEC")==0) + { + cout << "openEMS - dumping PEC info to 'PEC_dump.vtk'" << endl; + m_debugPEC = true; + return true; + } + else if (strcmp(argv,"--debug-CSX")==0) + { + cout << "openEMS - dumping CSX geometry to 'debugCSX.xml'" << endl; + m_debugCSX = true; + return true; + } + else if (strcmp(argv,"--engine=basic")==0) + { + cout << "openEMS - enabled basic engine" << endl; + m_engine = EngineType_Basic; + return true; + } + else if (strcmp(argv,"--engine=sse")==0) + { + cout << "openEMS - enabled sse engine" << endl; + m_engine = EngineType_SSE; + return true; + } + else if (strcmp(argv,"--engine=sse-compressed")==0) + { + cout << "openEMS - enabled compressed sse engine" << endl; + m_engine = EngineType_SSE_Compressed; + return true; + } + else if (strcmp(argv,"--engine=multithreaded")==0) + { + cout << "openEMS - enabled multithreading" << endl; + m_engine = EngineType_Multithreaded; + return true; + } + else if (strncmp(argv,"--numThreads=",13)==0) + { + m_engine_numThreads = atoi(argv+13); + cout << "openEMS - fixed number of threads: " << m_engine_numThreads << endl; + return true; + } + else if (strcmp(argv,"--engine=fastest")==0) + { + cout << "openEMS - enabled multithreading engine" << endl; + m_engine = EngineType_Multithreaded; + return true; + } + else if (strcmp(argv,"--no-simulation")==0) + { + cout << "openEMS - disabling simulation => preprocessing only" << endl; + m_no_simulation = true; + return true; + } + else if (strcmp(argv,"--dump-statistics")==0) + { + cout << "openEMS - dump simulation statistics to '" << __OPENEMS_RUN_STAT_FILE__ << "' and '" << __OPENEMS_STAT_FILE__ << "'" << endl; + m_DumpStats = true; + return true; + } + + return false; +} + +string openEMS::GetExtLibsInfo(string prefix) +{ + stringstream str; + + str << prefix << "Used external libraries:" << endl; + str << prefix << "\t" << ContinuousStructure::GetInfoLine(true) << endl; + + // libhdf5 + unsigned int major, minor, release; + if (H5get_libversion( &major, &minor, &release ) >= 0) + { + str << prefix << "\t" << "hdf5 -- Version: " << major << '.' << minor << '.' << release << endl; + str << prefix << "\t" << " compiled against: " H5_VERS_INFO << endl; + } + + // tinyxml + str << prefix << "\t" << "tinyxml -- compiled against: " << TIXML_MAJOR_VERSION << '.' << TIXML_MINOR_VERSION << '.' << TIXML_PATCH_VERSION << endl; + + // fparser + str << prefix << "\t" << "fparser" << endl; + + // boost + str << prefix << "\t" << "boost -- compiled against: " << BOOST_LIB_VERSION << endl; + + //vtk + str << prefix << "\t" << "vtk -- Version: " << vtkVersion::GetVTKMajorVersion() << "." << vtkVersion::GetVTKMinorVersion() << "." << vtkVersion::GetVTKBuildVersion() << endl; + str << prefix << "\t" << " compiled against: " << VTK_VERSION << endl; + + return str.str(); +} + +void openEMS::WelcomeScreen() +{ +#if defined(_LP64) || defined(_WIN64) + string bits = "64bit"; +#else + string bits = "32bit"; +#endif + + cout << " ---------------------------------------------------------------------- " << endl; + cout << " | openEMS " << bits << " -- version " GIT_VERSION << endl; + cout << " | (C) 2010-2016 Thorsten Liebig GPL license" << endl; + cout << " ---------------------------------------------------------------------- " << endl; + cout << openEMS::GetExtLibsInfo("\t") << endl; +} + +bool openEMS::SetupBoundaryConditions() +{ + FDTD_Op->SetBoundaryCondition(m_BC_type); //operator only knows about PEC and PMC, everything else is defined by extensions (see below) + + /**************************** create all operator/engine extensions here !!!! **********************************/ + for (int n=0; n<6; ++n) + { + FDTD_Op->SetBCSize(n, 0); + if (m_BC_type[n]==2) //Mur-ABC + { + FDTD_Op->SetBCSize(n, 1); + Operator_Ext_Mur_ABC* op_ext_mur = new Operator_Ext_Mur_ABC(FDTD_Op); + op_ext_mur->SetDirection(n/2,n%2); + if (m_Mur_v_ph[n]>0) + op_ext_mur->SetPhaseVelocity(m_Mur_v_ph[n]); + FDTD_Op->AddExtension(op_ext_mur); + } + if (m_BC_type[n]==3) + FDTD_Op->SetBCSize(n, m_PML_size[n]); + } + + + //create the upml + Operator_Ext_UPML::Create_UPML(FDTD_Op, m_BC_type, m_PML_size, string()); + + return true; +} + +Engine_Interface_FDTD* openEMS::NewEngineInterface(int multigridlevel) +{ + Operator_CylinderMultiGrid* op_cyl_mg = dynamic_cast(FDTD_Op); + while (op_cyl_mg && multigridlevel>0) + { + int mgl = op_cyl_mg->GetMultiGridLevel(); + if (mgl==multigridlevel) + { + if (g_settings.GetVerboseLevel()>0) + cout << __func__ << ": Operator with requested multi-grid level found." << endl; + return new Engine_Interface_Cylindrical_FDTD(op_cyl_mg); + } + Operator_Cylinder* op_cyl_inner = op_cyl_mg->GetInnerOperator(); + op_cyl_mg = dynamic_cast(op_cyl_inner); + if (op_cyl_mg==NULL) //inner most operator reached + { + if (g_settings.GetVerboseLevel()>0) + cout << __func__ << ": Operator with highest multi-grid level chosen." << endl; + return new Engine_Interface_Cylindrical_FDTD(op_cyl_inner); + } + // try next level + } + Operator_Cylinder* op_cyl = dynamic_cast(FDTD_Op); + if (op_cyl) + return new Engine_Interface_Cylindrical_FDTD(op_cyl); + Operator_sse* op_sse = dynamic_cast(FDTD_Op); + if (op_sse) + return new Engine_Interface_SSE_FDTD(op_sse); + return new Engine_Interface_FDTD(FDTD_Op); +} + +void openEMS::SetVerboseLevel(int level) +{ + g_settings.SetVerboseLevel(level); +} + +bool openEMS::SetupProcessing() +{ + //*************** setup processing ************// + if (g_settings.GetVerboseLevel()>0) + cout << "Setting up processing..." << endl; + + unsigned int Nyquist = FDTD_Op->GetExcitationSignal()->GetNyquistNum(); + PA = new ProcessingArray(Nyquist); + + double start[3]; + double stop[3]; + bool l_MultiBox = false; + vector Probes = m_CSX->GetPropertyByType(CSProperties::PROBEBOX); + for (size_t i=0; iGetQtyPrimitives()>1); + + for (size_t nb=0; nbGetQtyPrimitives(); ++nb) + { + CSPrimitives* prim = Probes.at(i)->GetPrimitive(nb); + if (prim!=NULL) + { + double bnd[6] = {0,0,0,0,0,0}; + prim->GetBoundBox(bnd,true); + start[0]= bnd[0]; + start[1]=bnd[2]; + start[2]=bnd[4]; + stop[0] = bnd[1]; + stop[1] =bnd[3]; + stop[2] =bnd[5]; + CSPropProbeBox* pb = Probes.at(i)->ToProbeBox(); + ProcessIntegral* proc = NULL; + if (pb) + { + if (pb->GetProbeType()==0) + { + ProcessVoltage* procVolt = new ProcessVoltage(NewEngineInterface()); + proc=procVolt; + } + else if (pb->GetProbeType()==1) + { + ProcessCurrent* procCurr = new ProcessCurrent(NewEngineInterface()); + proc=procCurr; + } + else if (pb->GetProbeType()==2) + proc = new ProcessFieldProbe(NewEngineInterface(),0); + else if (pb->GetProbeType()==3) + proc = new ProcessFieldProbe(NewEngineInterface(),1); + else if ((pb->GetProbeType()==10) || (pb->GetProbeType()==11)) + { + ProcessModeMatch* pmm = new ProcessModeMatch(NewEngineInterface()); + pmm->SetFieldType(pb->GetProbeType()-10); + pmm->SetModeFunction(0,pb->GetAttributeValue("ModeFunctionX")); + pmm->SetModeFunction(1,pb->GetAttributeValue("ModeFunctionY")); + pmm->SetModeFunction(2,pb->GetAttributeValue("ModeFunctionZ")); + proc = pmm; + } + else + { + cerr << "openEMS::SetupFDTD: Warning: Probe type " << pb->GetProbeType() << " of property '" << pb->GetName() << "' is unknown..." << endl; + continue; + } + if (CylinderCoords) + proc->SetMeshType(Processing::CYLINDRICAL_MESH); + if ((pb->GetProbeType()==1) || (pb->GetProbeType()==3) || (pb->GetProbeType()==11)) + { + proc->SetDualTime(true); + proc->SetDualMesh(true); + } + proc->SetProcessInterval(Nyquist/m_OverSampling); + if (pb->GetStartTime()>0 || pb->GetStopTime()>0) + proc->SetProcessStartStopTime(pb->GetStartTime(), pb->GetStopTime()); + proc->AddFrequency(pb->GetFDSamples()); + proc->GetNormalDir(pb->GetNormalDir()); + if (l_MultiBox==false) + proc->SetName(pb->GetName()); + else + proc->SetName(pb->GetName(),nb); + proc->DefineStartStopCoord(start,stop); + if (g_settings.showProbeDiscretization()) + proc->ShowSnappedCoords(); + proc->SetWeight(pb->GetWeighting()); + PA->AddProcessing(proc); + prim->SetPrimitiveUsed(true); + } + else + delete proc; + } + } + } + + vector DumpProps = m_CSX->GetPropertyByType(CSProperties::DUMPBOX); + for (size_t i=0; iGetQtyPrimitives()>1); + + for (size_t nb=0; nbGetQtyPrimitives(); ++nb) + { + + CSPrimitives* prim = DumpProps.at(i)->GetPrimitive(nb); + if (prim!=NULL) + { + double bnd[6] = {0,0,0,0,0,0}; + prim->GetBoundBox(bnd,true); + start[0]= bnd[0]; + start[1]=bnd[2]; + start[2]=bnd[4]; + stop[0] = bnd[1]; + stop[1] =bnd[3]; + stop[2] =bnd[5]; + CSPropDumpBox* db = DumpProps.at(i)->ToDumpBox(); + if (db) + { + if ((db->GetDumpType()>=0) && (db->GetDumpType()<=3)) + ProcField = new ProcessFieldsTD(NewEngineInterface(db->GetMultiGridLevel())); + else if ((db->GetDumpType()>=10) && (db->GetDumpType()<=13)) + ProcField = new ProcessFieldsFD(NewEngineInterface(db->GetMultiGridLevel())); + else if ( ((db->GetDumpType()>=20) && (db->GetDumpType()<=22)) || (db->GetDumpType()==29) ) + { + ProcessFieldsSAR* procSAR = new ProcessFieldsSAR(NewEngineInterface(db->GetMultiGridLevel())); + ProcField = procSAR; + string method = db->GetAttributeValue("SAR_Method"); + if (!method.empty()) + procSAR->SetSARAveragingMethod(method); + // use (center)-cell based conductivity only + procSAR->SetUseCellConductivity(true); + } + else + cerr << "openEMS::SetupFDTD: unknown dump box type... skipping!" << endl; + if (ProcField) + { + ProcField->SetEnable(Enable_Dumps); + ProcField->SetProcessInterval(Nyquist/m_OverSampling); + if (db->GetStopTime()>0 || db->GetStartTime()>0) + ProcField->SetProcessStartStopTime(db->GetStartTime(), db->GetStopTime()); + if ((db->GetDumpType()==1) || (db->GetDumpType()==11)) + { + ProcField->SetDualTime(true); + //make dualMesh the default mesh for h-field dumps, maybe overwritten by interpolation type (node-interpolation) + ProcField->SetDualMesh(true); + } + if (db->GetDumpType()>=10) + { + ProcField->AddFrequency(db->GetFDSamples()); + ProcField->SetDumpType((ProcessFields::DumpType)(db->GetDumpType()-10)); + } + else + ProcField->SetDumpType((ProcessFields::DumpType)db->GetDumpType()); + + if (db->GetDumpType()==20) + ProcField->SetDumpType(ProcessFields::SAR_LOCAL_DUMP); + if (db->GetDumpType()==21) + ProcField->SetDumpType(ProcessFields::SAR_1G_DUMP); + if (db->GetDumpType()==22) + ProcField->SetDumpType(ProcessFields::SAR_10G_DUMP); + if (db->GetDumpType()==29) + ProcField->SetDumpType(ProcessFields::SAR_RAW_DATA); + + //SetupMaterialStorages() has previewed storage needs... refresh here to prevent cleanup!!! + if ( ProcField->NeedConductivity() && Enable_Dumps ) + FDTD_Op->SetMaterialStoreFlags(1,true); + + ProcField->SetDumpMode((Engine_Interface_Base::InterpolationType)db->GetDumpMode()); + ProcField->SetFileType((ProcessFields::FileType)db->GetFileType()); + if (CylinderCoords) + ProcField->SetMeshType(Processing::CYLINDRICAL_MESH); + if (db->GetSubSampling()) + for (int n=0; n<3; ++n) + ProcField->SetSubSampling(db->GetSubSampling(n),n); + if (db->GetOptResolution()) + for (int n=0; n<3; ++n) + ProcField->SetOptResolution(db->GetOptResolution(n),n); + + if (l_MultiBox==false) + ProcField->SetName(db->GetName()); + else + ProcField->SetName(db->GetName(),nb); + + ProcField->SetFileName(ProcField->GetName()); + ProcField->DefineStartStopCoord(start,stop); + if (g_settings.showProbeDiscretization()) + ProcField->ShowSnappedCoords(); + PA->AddProcessing(ProcField); + prim->SetPrimitiveUsed(true); + } + } + } + } + } + + return true; +} + +bool openEMS::SetupMaterialStorages() +{ + vector DumpProps = m_CSX->GetPropertyByType(CSProperties::DUMPBOX); + for (size_t i=0; iToDumpBox(); + if (!db) + continue; + if (db->GetQtyPrimitives()==0) + continue; + //check for current density dump types + if ( ((db->GetDumpType()==2) || (db->GetDumpType()==12) || // current density storage + (db->GetDumpType()==20) || (db->GetDumpType()==21) || (db->GetDumpType()==22)) && // SAR dump types + Enable_Dumps ) + FDTD_Op->SetMaterialStoreFlags(1,true); //tell operator to store kappa material data + } + return true; +} + +void openEMS::SetupCylinderMultiGrid(std::string val) +{ + m_CC_MultiGrid.clear(); + m_CC_MultiGrid = SplitString2Double(val,','); +} + +bool openEMS::SetupOperator() +{ + if (CylinderCoords) + { + if (m_CC_MultiGrid.size()>0) + { + FDTD_Op = Operator_CylinderMultiGrid::New(m_CC_MultiGrid, m_engine_numThreads); + if (FDTD_Op==NULL) + FDTD_Op = Operator_Cylinder::New(m_engine_numThreads); + } + else + FDTD_Op = Operator_Cylinder::New(m_engine_numThreads); + } + else if (m_engine == EngineType_SSE) + { + FDTD_Op = Operator_sse::New(); + } + else if (m_engine == EngineType_SSE_Compressed) + { + FDTD_Op = Operator_SSE_Compressed::New(); + } + else if (m_engine == EngineType_Multithreaded) + { + FDTD_Op = Operator_Multithread::New(m_engine_numThreads); + } + else + { + FDTD_Op = Operator::New(); + } + return true; +} + +void openEMS::Set_BC_Type(int idx, int type) +{ + if ((idx<0) || (idx>5)) + return; + m_BC_type[idx] = type; +} + +int openEMS::Get_BC_Type(int idx) +{ + if ((idx<0) || (idx>5)) + return -1; + return m_BC_type[idx]; +} + +void openEMS::Set_BC_PML(int idx, unsigned int size) +{ + if ((idx<0) || (idx>5)) + return; + m_BC_type[idx] = 3; + m_PML_size[idx] = size; +} + +int openEMS::Get_PML_Size(int idx) +{ + if ((idx<0) || (idx>5)) + return -1; + if (m_BC_type[idx]!=3) + return -1; // return -1 if BC was *not* a PML + return m_PML_size[idx]; +} + +void openEMS::Set_Mur_PhaseVel(int idx, double val) +{ + if ((idx<0) || (idx>5)) + return; + m_Mur_v_ph[idx] = val; +} + +bool openEMS::ParseFDTDSetup(std::string file) +{ + Reset(); + + if (g_settings.GetVerboseLevel()>0) + cout << "Read openEMS xml file: " << file << " ..." << endl; + + TiXmlDocument doc(file); + if (!doc.LoadFile()) + { + cerr << "openEMS: Error File-Loading failed!!! File: " << file << endl; + exit(-1); + } + + if (g_settings.GetVerboseLevel()>0) + cout << "Read openEMS Settings..." << endl; + TiXmlElement* openEMSxml = doc.FirstChildElement("openEMS"); + if (openEMSxml==NULL) + { + cerr << "Can't read openEMS ... " << endl; + exit(-1); + } + TiXmlElement* FDTD_Opts = openEMSxml->FirstChildElement("FDTD"); + + if (FDTD_Opts==NULL) + { + cerr << "Can't read openEMS FDTD Settings... " << endl; + exit(-1); + } + + if (g_settings.GetVerboseLevel()>0) + cout << "Read Geometry..." << endl; + ContinuousStructure* csx = new ContinuousStructure(); + string EC(csx->ReadFromXML(openEMSxml)); + if (EC.empty()==false) + cerr << EC << endl; + this->SetCSX(csx); + + + return this->Parse_XML_FDTDSetup(FDTD_Opts); +} + +bool openEMS::Parse_XML_FDTDSetup(TiXmlElement* FDTD_Opts) +{ + double dhelp=0; + FDTD_Opts->QueryDoubleAttribute("NumberOfTimesteps",&dhelp); + if (dhelp<0) + this->SetNumberOfTimeSteps(0); + else + this->SetNumberOfTimeSteps((unsigned int)dhelp); + + int ihelp = 0; + FDTD_Opts->QueryIntAttribute("CylinderCoords",&ihelp); + if (ihelp==1) + { + this->SetCylinderCoords(true); + const char* cchelp = FDTD_Opts->Attribute("MultiGrid"); + if (cchelp!=NULL) + this->SetupCylinderMultiGrid(string(cchelp)); + } + + dhelp = 0; + FDTD_Opts->QueryDoubleAttribute("endCriteria",&dhelp); + if (dhelp==0) + this->SetEndCriteria(1e-6); + else + this->SetEndCriteria(dhelp); + + ihelp = 0; + FDTD_Opts->QueryIntAttribute("OverSampling",&ihelp); + if (ihelp<2) + this->SetOverSampling(2); + else + this->SetOverSampling(ihelp); + + // check for cell constant material averaging + if (FDTD_Opts->QueryIntAttribute("CellConstantMaterial",&ihelp)==TIXML_SUCCESS) + this->SetCellConstantMaterial(ihelp==1); + + TiXmlElement* BC = FDTD_Opts->FirstChildElement("BoundaryCond"); + if (BC==NULL) + { + cerr << "Can't read openEMS boundary cond Settings... " << endl; + exit(-3); + } + +// const char* tmp = BC->Attribute("PML_Grading"); +// string pml_gradFunc; +// if (tmp) +// pml_gradFunc = string(tmp); + + string bound_names[] = {"xmin","xmax","ymin","ymax","zmin","zmax"}; + string s_bc; + for (int n=0; n<6; ++n) + { + int EC = BC->QueryIntAttribute(bound_names[n].c_str(),&ihelp); + if (EC==TIXML_SUCCESS) + { + this->Set_BC_Type(n, ihelp); + continue; + } + if (EC==TIXML_WRONG_TYPE) + { + const char* tmp = BC->Attribute(bound_names[n].c_str()); + if (tmp) + s_bc = string(tmp); + else + cerr << "openEMS::SetupBoundaryConditions: Warning, boundary condition for \"" << bound_names[n] << "\" unknown... set to PEC " << endl; + if (s_bc=="PEC") + this->Set_BC_Type(n, 0); + else if (s_bc=="PMC") + this->Set_BC_Type(n, 1); + else if (s_bc=="MUR") + this->Set_BC_Type(n, 2); + else if (strncmp(s_bc.c_str(),"PML_=",4)==0) + this->Set_BC_PML(n, atoi(s_bc.c_str()+4)); + else + cerr << "openEMS::SetupBoundaryConditions: Warning, boundary condition for \"" << bound_names[n] << "\" unknown... set to PEC " << endl; + } + else + cerr << "openEMS::SetupBoundaryConditions: Warning, boundary condition for \"" << bound_names[n] << "\" not found... set to PEC " << endl; + } + + //read general mur phase velocity + if (BC->QueryDoubleAttribute("MUR_PhaseVelocity",&dhelp) == TIXML_SUCCESS) + for (int n=0;n<6;++n) + this->Set_Mur_PhaseVel(n, dhelp); + + string mur_v_ph_names[6] = {"MUR_PhaseVelocity_xmin", "MUR_PhaseVelocity_xmax", "MUR_PhaseVelocity_ymin", "MUR_PhaseVelocity_ymax", "MUR_PhaseVelocity_zmin", "MUR_PhaseVelocity_zmax"}; + for (int n=0; n<6; ++n) + if (BC->QueryDoubleAttribute(mur_v_ph_names[n].c_str(),&dhelp) == TIXML_SUCCESS) + this->Set_Mur_PhaseVel(n, dhelp); + + TiXmlElement* m_Excite_Elem = FDTD_Opts->FirstChildElement("Excitation"); + if (!m_Excite_Elem) + { + cerr << "Excitation::setupExcitation: Error, can't read openEMS excitation settings... " << endl; + return false; + } + + Excitation* exc = this->InitExcitation(); + double f0=0, fc=0, f_max=0; + ihelp = -1; + m_Excite_Elem->QueryIntAttribute("Type",&ihelp); + switch (ihelp) + { + case Excitation::GaissianPulse: + m_Excite_Elem->QueryDoubleAttribute("f0",&f0); + m_Excite_Elem->QueryDoubleAttribute("fc",&fc); + exc->SetupGaussianPulse(f0, fc); + break; + case Excitation::Sinusoidal: // sinusoidal excite + m_Excite_Elem->QueryDoubleAttribute("f0",&f0); + exc->SetupSinusoidal(f0); + break; + case Excitation::DiracPulse: + FDTD_Opts->QueryDoubleAttribute("f_max",&f_max); + exc->SetupDiracPulse(f_max); + break; + case Excitation::Step: + FDTD_Opts->QueryDoubleAttribute("f_max",&f_max); + exc->SetupStepExcite(f_max); + break; + case Excitation::CustomExcite: + m_Excite_Elem->QueryDoubleAttribute("f0",&f0); + FDTD_Opts->QueryDoubleAttribute("f_max",&f_max); + exc->SetupCustomExcite(m_Excite_Elem->Attribute("Function"), f0, f_max); + break; + } + + if (FDTD_Opts->QueryIntAttribute("TimeStepMethod",&ihelp)==TIXML_SUCCESS) + this->SetTimeStepMethod(ihelp); + if (FDTD_Opts->QueryDoubleAttribute("TimeStep",&dhelp)==TIXML_SUCCESS) + this->SetTimeStep(dhelp); + if (FDTD_Opts->QueryDoubleAttribute("TimeStepFactor",&dhelp)==TIXML_SUCCESS) + this->SetTimeStepFactor(dhelp); +} + +void openEMS::SetGaussExcite(double f0, double fc) +{ + this->InitExcitation(); + m_Exc->SetupGaussianPulse(f0, fc); +} + +Excitation* openEMS::InitExcitation() +{ + delete m_Exc; + m_Exc = new Excitation(); + return m_Exc; +} + +void openEMS::SetCSX(ContinuousStructure* csx) +{ + delete m_CSX; + m_CSX = csx; +} + +int openEMS::SetupFDTD() +{ + timeval startTime; + gettimeofday(&startTime,NULL); + + if (m_CSX==NULL) + { + cerr << "openEMS::SetupFDTD: Error: CSXCAD is not set!" << endl; + return 3; + } + if (m_CSX==NULL) + { + cerr << "openEMS::SetupFDTD: Error: CSXCAD is not set!" << endl; + return 3; + } + std::string ec = m_CSX->Update(); + if (!ec.empty()) + cerr << ec << endl; + if (g_settings.GetVerboseLevel()>2) + m_CSX->ShowPropertyStatus(cerr); + + if (CylinderCoords) + if (m_CSX->GetCoordInputType()!=CYLINDRICAL) + { + cerr << "openEMS::SetupFDTD: Warning: Coordinate system found in the CSX file is not a cylindrical. Forcing to cylindrical coordinate system!" << endl; + m_CSX->SetCoordInputType(CYLINDRICAL); //tell CSX to use cylinder-coords + } + + if (m_debugCSX) + m_CSX->Write2XML("debugCSX.xml"); + + //*************** setup operator ************// + if (SetupOperator()==false) + return 2; + + // default material averaging is quarter cell averaging + FDTD_Op->SetQuarterCellMaterialAvg(); + + if (m_CellConstantMaterial) + { + FDTD_Op->SetCellConstantMaterial(); + if (g_settings.GetVerboseLevel()>0) + cout << "Enabling constant cell material assumption." << endl; + } + + if (m_Exc==NULL) + { + cerr << "openEMS::SetupFDTD: Error, excitation is not defined! Abort!" << endl; + return 3; + } + + FDTD_Op->SetExcitationSignal(m_Exc); + FDTD_Op->AddExtension(new Operator_Ext_Excitation(FDTD_Op)); + if (!CylinderCoords) + FDTD_Op->AddExtension(new Operator_Ext_TFSF(FDTD_Op)); + + if (FDTD_Op->SetGeometryCSX(m_CSX)==false) return(2); + + SetupBoundaryConditions(); + + FDTD_Op->SetTimeStepMethod(m_TS_method); + + if (m_TS>0) + FDTD_Op->SetTimestep(m_TS); + if (m_TS_fac<1) + FDTD_Op->SetTimestepFactor(m_TS_fac); + + // Is a steady state detection requested + Operator_Ext_SteadyState* Op_Ext_SSD = NULL; + if (m_Exc->GetSignalPeriod()>0) + { + cout << "Create a steady state detection using a period of " << m_Exc->GetSignalPeriod() << " s" << endl; + Op_Ext_SSD = new Operator_Ext_SteadyState(FDTD_Op, m_Exc->GetSignalPeriod()); + unsigned int pos[3]; + for (int p=0;p<3;++p) + pos[p] = FDTD_Op->GetNumberOfLines(p)/2; + Op_Ext_SSD->Add_E_Probe(pos, 0); + Op_Ext_SSD->Add_E_Probe(pos, 1); + Op_Ext_SSD->Add_E_Probe(pos, 2); + + for (int n=0;n<3;++n) + { + for (int p=0;p<3;++p) + pos[p] = FDTD_Op->GetNumberOfLines(p)/2; + + pos[n] *= 1/4; + Op_Ext_SSD->Add_E_Probe(pos, 0); + Op_Ext_SSD->Add_E_Probe(pos, 1); + Op_Ext_SSD->Add_E_Probe(pos, 2); + + pos[n] *= 3/4; + Op_Ext_SSD->Add_E_Probe(pos, 0); + Op_Ext_SSD->Add_E_Probe(pos, 1); + Op_Ext_SSD->Add_E_Probe(pos, 2); + } + FDTD_Op->AddExtension(Op_Ext_SSD); + } + + if ((m_CSX->GetQtyPropertyType(CSProperties::LORENTZMATERIAL)>0) || (m_CSX->GetQtyPropertyType(CSProperties::DEBYEMATERIAL)>0)) + FDTD_Op->AddExtension(new Operator_Ext_LorentzMaterial(FDTD_Op)); + if (m_CSX->GetQtyPropertyType(CSProperties::CONDUCTINGSHEET)>0) + FDTD_Op->AddExtension(new Operator_Ext_ConductingSheet(FDTD_Op, m_Exc->GetMaxFreq())); + + //check all properties to request material storage during operator creation... + SetupMaterialStorages(); + + /******************* create the EC-FDTD operator *****************************/ + Operator::DebugFlags debugFlags = Operator::None; + if (DebugMat) + debugFlags |= Operator::debugMaterial; + if (DebugOp) + debugFlags |= Operator::debugOperator; + if (m_debugPEC) + debugFlags |= Operator::debugPEC; + + FDTD_Op->CalcECOperator( debugFlags ); + /*******************************************************************************/ + + //reset flags for material storage, if no dump-box resets it to true, it will be cleaned up... + FDTD_Op->SetMaterialStoreFlags(0,false); + FDTD_Op->SetMaterialStoreFlags(1,false); + FDTD_Op->SetMaterialStoreFlags(2,false); + FDTD_Op->SetMaterialStoreFlags(3,false); + + unsigned int maxTime_TS = (unsigned int)(m_maxTime/FDTD_Op->GetTimestep()); + if ((m_maxTime>0) && (maxTime_TSbuildExcitationSignal(NrTS)) + exit(2); + m_Exc->DumpVoltageExcite("et"); + m_Exc->DumpCurrentExcite("ht"); + + timeval OpDoneTime; + gettimeofday(&OpDoneTime,NULL); + + if (g_settings.GetVerboseLevel()>0) + { + FDTD_Op->ShowStat(); + FDTD_Op->ShowExtStat(); + cout << "Creation time for operator: " << CalcDiffTime(OpDoneTime,startTime) << " s" << endl; + } + cout << "FDTD simulation size: " << FDTD_Op->GetNumberOfLines(0) << "x" << FDTD_Op->GetNumberOfLines(1) << "x" << FDTD_Op->GetNumberOfLines(2) << " --> " << FDTD_Op->GetNumberCells() << " FDTD cells " << endl; + cout << "FDTD timestep is: " <GetTimestep() << " s; Nyquist rate: " << m_Exc->GetNyquistNum() << " timesteps @" << CalcNyquistFrequency(m_Exc->GetNyquistNum(),FDTD_Op->GetTimestep()) << " Hz" << endl; + if (m_Exc->GetNyquistNum()>1000) + cerr << "openEMS::SetupFDTD: Warning, the timestep seems to be very small --> long simulation. Check your mesh!?" << endl; + + if (m_Exc->GetSignalPeriod()==0) + { + cout << "Excitation signal length is: " << m_Exc->GetLength() << " timesteps (" << m_Exc->GetLength()*FDTD_Op->GetTimestep() << "s)" << endl; + cout << "Max. number of timesteps: " << NrTS << " ( --> " << (double)NrTS/(double)(m_Exc->GetLength()) << " * Excitation signal length)" << endl; + if ( ((double)NrTS/(double)m_Exc->GetLength() < 3) && (m_Exc->GetExciteType()==0)) + cerr << "openEMS::SetupFDTD: Warning, max. number of timesteps is smaller than three times the excitation. " << endl << \ + "\tYou may want to choose a higher number of max. timesteps... " << endl; + } + else + { + int p = int(m_Exc->GetSignalPeriod()/FDTD_Op->GetTimestep()); + cout << "Excitation signal period is: " << p << " timesteps (" << m_Exc->GetSignalPeriod() << "s)" << endl; + cout << "Max. number of timesteps: " << NrTS << " ( --> " << (double)NrTS/(double)(m_Exc->GetLength()) << " * Excitation signal period)" << endl; + if (NrTS/p < 3) + cerr << "openEMS::SetupFDTD: Warning, max. number of timesteps is smaller than three times the excitation signal period. " << endl << \ + "\tYou may want to choose a higher number of max. timesteps... " << endl; + } + + if (m_no_simulation) + { + // simulation was disabled (to generate debug output only) + return 1; + } + + //create FDTD engine + FDTD_Eng = FDTD_Op->CreateEngine(); + + if (Op_Ext_SSD) + { + Eng_Ext_SSD = dynamic_cast(Op_Ext_SSD->GetEngineExtention()); + Eng_Ext_SSD->SetEngineInterface(this->NewEngineInterface()); + } + + //setup all processing classes + if (SetupProcessing()==false) + return 2; + + // Cleanup all unused material storages... + FDTD_Op->CleanupMaterialStorage(); + + //check and warn for unused properties and primitives + m_CSX->WarnUnusedPrimitves(cerr); + + // dump all boxes (voltage, current, fields, ...) + if (m_debugBox) + { + PA->DumpBoxes2File("box_dump_"); + } + + return 0; +} + +string FormatTime(int sec) +{ + stringstream ss; + if (sec<60) + { + ss << setw(9) << sec << "s"; + return ss.str(); + } + if (sec<3600) + { + ss << setw(6) << sec/60 << "m" << setw(2) << setfill('0') << sec%60 << "s"; + return ss.str(); + } + ss << setw(3) << sec/3600 << "h" << setw(2) << setfill('0') << (sec%3600)/60 << "m" << setw(2) << setfill('0') << sec%60 << "s"; + return ss.str(); +} + +bool openEMS::CheckAbortCond() +{ + if (m_Abort) //abort was set externally + return true; + + //check whether the file "ABORT" exist in current working directory + ifstream ifile("ABORT"); + if (ifile) + { + ifile.close(); + cerr << "openEMS::CheckAbortCond(): Found file \"ABORT\", aborting simulation..." << endl; + return true; + } + + return false; +} + +void openEMS::RunFDTD() +{ + cout << "Running FDTD engine... this may take a while... grab a cup of coffee?!?" << endl; + + //special handling of a field processing, needed to realize the end criteria... + ProcessFields* ProcField = new ProcessFields(NewEngineInterface()); + PA->AddProcessing(ProcField); + double maxE=0,currE=0; + + //init processings + PA->InitAll(); + + //add all timesteps to end-crit field processing with max excite amplitude + unsigned int maxExcite = FDTD_Op->GetExcitationSignal()->GetMaxExcitationTimestep(); +// for (unsigned int n=0; nExc->Volt_Count; ++n) +// ProcField->AddStep(FDTD_Op->Exc->Volt_delay[n]+maxExcite); + ProcField->AddStep(maxExcite); + + double change=1; + int prevTS=0,currTS=0; + double numCells = FDTD_Op->GetNumberCells(); + double speed = 0; + double t_diff; + double t_run; + + timeval currTime; + gettimeofday(&currTime,NULL); + timeval startTime = currTime; + timeval prevTime= currTime; + + if (m_DumpStats) + InitRunStatistics(__OPENEMS_RUN_STAT_FILE__); + //*************** simulate ************// + + PA->PreProcess(); + int step=PA->Process(); + if ((step<0) || (step>(int)NrTS)) step=NrTS; + while ((FDTD_Eng->GetNumberOfTimesteps()endCrit) && !CheckAbortCond()) + { + FDTD_Eng->IterateTS(step); + step=PA->Process(); + + if ((Eng_Ext_SSD==NULL) && ProcField->CheckTimestep()) + { + currE = ProcField->CalcTotalEnergyEstimate(); + if (currE>maxE) + maxE=currE; + } + +// cout << " do " << step << " steps; current: " << eng.GetNumberOfTimesteps() << endl; + currTS = FDTD_Eng->GetNumberOfTimesteps(); + if ((step<0) || (step>(int)(NrTS - currTS))) step=NrTS - currTS; + + gettimeofday(&currTime,NULL); + + t_diff = CalcDiffTime(currTime,prevTime); + + if (t_diff>4) + { + t_run = CalcDiffTime(currTime,startTime); + speed = numCells*(currTS-prevTS)/t_diff; + cout << "[@" << FormatTime(t_run) << "] Timestep: " << setw(12) << currTS ; + cout << " || Speed: " << setw(6) << setprecision(1) << std::fixed << speed*1e-6 << " MC/s (" << setw(4) << setprecision(3) << std::scientific << t_diff/(currTS-prevTS) << " s/TS)" ; + if (Eng_Ext_SSD==NULL) + { + currE = ProcField->CalcTotalEnergyEstimate(); + if (currE>maxE) + maxE=currE; + if (maxE) + change = currE/maxE; + cout << " || Energy: ~" << setw(6) << setprecision(2) << std::scientific << currE << " (-" << setw(5) << setprecision(2) << std::fixed << fabs(10.0*log10(change)) << "dB)" << endl; + } + else + { + change = Eng_Ext_SSD->GetLastDiff(); + cout << " || SteadyState: " << setw(6) << setprecision(2) << std::fixed << 10.0*log10(change) << " dB" << endl; + } + prevTime=currTime; + prevTS=currTS; + + PA->FlushNext(); + + if (m_DumpStats) + DumpRunStatistics(__OPENEMS_RUN_STAT_FILE__, t_run, currTS, speed, currE); + } + } + if ((change>endCrit) && (FDTD_Op->GetExcitationSignal()->GetExciteType()==0)) + cerr << "RunFDTD: Warning: Max. number of timesteps was reached before the end-criteria of -" << fabs(10.0*log10(endCrit)) << "dB was reached... " << endl << \ + "\tYou may want to choose a higher number of max. timesteps... " << endl; + + gettimeofday(&currTime,NULL); + t_diff = CalcDiffTime(currTime,startTime); + + cout << "Time for " << FDTD_Eng->GetNumberOfTimesteps() << " iterations with " << FDTD_Op->GetNumberCells() << " cells : " << t_diff << " sec" << endl; + cout << "Speed: " << numCells*(double)FDTD_Eng->GetNumberOfTimesteps()/t_diff*1e-6 << " MCells/s " << endl; + + if (m_DumpStats) + DumpStatistics(__OPENEMS_STAT_FILE__, t_diff); + + //*************** postproc ************// + PA->PostProcess(); +} + +bool openEMS::DumpStatistics(const string& filename, double time) +{ + ofstream stat_file; + stat_file.open(filename.c_str()); + + if (!stat_file.is_open()) + { + cerr << "openEMS::DumpStatistics: Error, opening file failed..." << endl; + return false; + } + stat_file << std::setprecision( 16 ); + stat_file << FDTD_Op->GetNumberCells() << "\t% number of cells" << endl; + stat_file << FDTD_Op->GetTimestep() << "\t% timestep (s)" << endl; + stat_file << FDTD_Eng->GetNumberOfTimesteps() << "\t% number of iterations" << endl; + stat_file << FDTD_Eng->GetNumberOfTimesteps()*FDTD_Op->GetTimestep() << "\t% total numercial time (s)" << endl; + stat_file << time << "\t% simulation time (s)" << endl; + stat_file << (double)FDTD_Op->GetNumberCells()*(double)FDTD_Eng->GetNumberOfTimesteps()/time << "\t% speed (cells/s)" << endl; + + stat_file.close(); + return true; +} + +bool openEMS::InitRunStatistics(const string& filename) +{ + ofstream stat_file; + stat_file.open(filename.c_str(), ios_base::out); + + if (!stat_file.is_open()) + { + cerr << "openEMS::InitRunStatistics: Error, opening file failed..." << endl; + return false; + } + stat_file << "%time\ttimestep\tspeed\tenergy" << endl; + stat_file.close(); + return true; +} + +bool openEMS::DumpRunStatistics(const string& filename, double time, unsigned int ts, double speed, double energy) +{ + ofstream stat_file; + stat_file.open(filename.c_str(), ios_base::app); + + if (!stat_file.is_open()) + { + cerr << "openEMS::DumpRunStatistics: Error, opening file failed..." << endl; + return false; + } + stat_file << time << "\t" << ts << "\t" << speed << "\t" << energy << endl; + stat_file.close(); + return true; +} -- cgit v1.2.3