diff --git a/FDTD/excitation.cpp b/FDTD/excitation.cpp index cc16078..8b0ceb6 100644 --- a/FDTD/excitation.cpp +++ b/FDTD/excitation.cpp @@ -28,18 +28,16 @@ Excitation::Excitation() Signal_volt = 0; Signal_curr = 0; - m_Excit_Type = -1; + this->Reset(0); - dT = 0; - m_nyquistTS = 0; - m_f_max = 0; - m_foi = 0; + m_Excite_Elem = NULL; + m_Excit_Type = -1; + m_SignalPeriod = 0; } Excitation::~Excitation() { - delete[] Signal_volt; - delete[] Signal_curr; + this->Reset(0); } void Excitation::Reset( double timestep ) @@ -49,15 +47,13 @@ void Excitation::Reset( double timestep ) delete[] Signal_curr; Signal_curr = 0; - m_Excit_Type = -1; - dT = timestep; m_nyquistTS = 0; m_f_max = 0; m_foi = 0; } -bool Excitation::setupExcitation( TiXmlElement* Excite, unsigned int maxTS ) +bool Excitation::setupExcitation( TiXmlElement* Excite) { if (!Excite) { @@ -65,6 +61,23 @@ bool Excitation::setupExcitation( TiXmlElement* Excite, unsigned int maxTS ) return false; } + m_Excite_Elem = Excite; + + double f0=0; + m_Excite_Elem->QueryIntAttribute("Type",&m_Excit_Type); + m_SignalPeriod = 0; + switch (m_Excit_Type) + { + case 1: // sinusoidal excite + m_Excite_Elem->QueryDoubleAttribute("f0",&f0); + m_SignalPeriod = 1/f0; + break; + } + return true; +} + +bool Excitation::buildExcitationSignal(unsigned int maxTS) +{ if (dT<=0) { cerr << "Excitation::setupExcitation: Error, invalid timestep... " << endl; @@ -73,17 +86,15 @@ bool Excitation::setupExcitation( TiXmlElement* Excite, unsigned int maxTS ) double f0=0; double fc=0; - Excite->QueryIntAttribute("Type",&m_Excit_Type); - switch (m_Excit_Type) { case 0: - Excite->QueryDoubleAttribute("f0",&f0); - Excite->QueryDoubleAttribute("fc",&fc); + m_Excite_Elem->QueryDoubleAttribute("f0",&f0); + m_Excite_Elem->QueryDoubleAttribute("fc",&fc); CalcGaussianPulsExcitation(f0,fc,maxTS); break; case 1: - Excite->QueryDoubleAttribute("f0",&f0); + m_Excite_Elem->QueryDoubleAttribute("f0",&f0); CalcSinusExcitation(f0,maxTS); break; case 2: @@ -93,18 +104,18 @@ bool Excitation::setupExcitation( TiXmlElement* Excite, unsigned int maxTS ) CalcStepExcitation(); break; case 10: - Excite->QueryDoubleAttribute("f0",&f0); - CalcCustomExcitation(f0,maxTS,Excite->Attribute("Function")); + m_Excite_Elem->QueryDoubleAttribute("f0",&f0); + CalcCustomExcitation(f0,maxTS,m_Excite_Elem->Attribute("Function")); break; default: - cerr << "Excitation::setupExcitation: Unknown excitation type: \"" << m_Excit_Type<< "\" !!" << endl; + cerr << "Excitation::buildExcitationSignal: Unknown excitation type: \"" << m_Excit_Type<< "\" !!" << endl; m_Excit_Type = -1; return false; } if (GetNyquistNum() == 0) { - cerr << "Excitation::setupExcitation: Unknown error... excitation setup failed!!" << endl; + cerr << "Excitation::buildExcitationSignal: Unknown error... excitation setup failed!!" << endl; return false; } diff --git a/FDTD/excitation.h b/FDTD/excitation.h index d27c2ef..8ab99d7 100644 --- a/FDTD/excitation.h +++ b/FDTD/excitation.h @@ -32,7 +32,8 @@ public: virtual void Reset( double timestep ); - bool setupExcitation( TiXmlElement* Excite, unsigned int maxTS ); + bool setupExcitation(TiXmlElement* Excite); + bool buildExcitationSignal(unsigned int maxTS); //! Get the excitation timestep with the (first) max amplitude virtual unsigned int GetMaxExcitationTimestep() const; @@ -61,14 +62,20 @@ public: //! Get the frequency of interest double GetFrequencyOfInterest() const {return m_foi;} + //! Get the signal period, 0 if not a periodical signal + double GetSignalPeriod() const {return m_SignalPeriod;} + FDTD_FLOAT* GetVoltageSignal() const {return Signal_volt;} FDTD_FLOAT* GetCurrentSignal() const {return Signal_curr;} protected: double dT; unsigned int m_nyquistTS; + double m_SignalPeriod; int m_Excit_Type; + TiXmlElement* m_Excite_Elem; + //Excitation time-signal unsigned int Length; FDTD_FLOAT* Signal_volt; diff --git a/FDTD/extensions/CMakeLists.txt b/FDTD/extensions/CMakeLists.txt index 575c6a8..e725f1e 100644 --- a/FDTD/extensions/CMakeLists.txt +++ b/FDTD/extensions/CMakeLists.txt @@ -21,6 +21,8 @@ set(SOURCES engine_ext_excitation.cpp operator_ext_tfsf.cpp engine_ext_tfsf.cpp + operator_ext_steadystate.cpp + engine_ext_steadystate.cpp ) # FDTD/extensions lib diff --git a/FDTD/extensions/engine_ext_steadystate.cpp b/FDTD/extensions/engine_ext_steadystate.cpp new file mode 100644 index 0000000..4b9575c --- /dev/null +++ b/FDTD/extensions/engine_ext_steadystate.cpp @@ -0,0 +1,104 @@ +/* +* Copyright (C) 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 "engine_ext_steadystate.h" +#include "operator_ext_steadystate.h" +#include "FDTD/engine_sse.h" +#include "FDTD/engine_interface_fdtd.h" + +Engine_Ext_SteadyState::Engine_Ext_SteadyState(Operator_Ext_SteadyState* op_ext): Engine_Extension(op_ext) +{ + m_Op_SS = op_ext; + m_Priority = ENG_EXT_PRIO_STEADYSTATE; + + for (size_t n=0;nm_E_probe_dir.size();++n) + { + double* rec = new double[m_Op_SS->m_TS_period*2]; + m_E_records.push_back(rec); + } + m_last_max_diff = 1; + last_total_energy = 0; + m_Eng_Interface = NULL; +} + +Engine_Ext_SteadyState::~Engine_Ext_SteadyState() +{ + for (size_t n=0;nm_TS_period; + unsigned int TS = m_Eng->GetNumberOfTimesteps(); + unsigned int rel_pos = m_Eng->GetNumberOfTimesteps()%(2*p); + for (size_t n=0;nGetVolt(m_Op_SS->m_E_probe_dir.at(n), m_Op_SS->m_E_probe_pos[0].at(n), m_Op_SS->m_E_probe_pos[1].at(n), m_Op_SS->m_E_probe_pos[2].at(n)); + if ((TS%(m_Op_SS->m_TS_period)==0) && (TS>=2*p)) + { + bool no_valid = true; + m_last_max_diff = 0; + double curr_total_energy = m_Eng_Interface->CalcFastEnergy(); + if (last_total_energy>0) + { + m_last_max_diff = abs(curr_total_energy-last_total_energy)/last_total_energy; + no_valid = false; + } +// cerr << curr_total_energy << "/" << last_total_energy << "=" << abs(curr_total_energy-last_total_energy)/last_total_energy << endl; + last_total_energy = curr_total_energy; + unsigned int old_pos = 0; + unsigned int new_pos = p; + if (rel_pos<=p) + { + new_pos = 0; + old_pos = p; + } + //cerr << TS << "/" << rel_pos << ": one period complete, new_pos" << new_pos << " old pos: " << old_pos << endl; + for (size_t n=0;n0) + { + m_last_max_diff = max(m_last_max_diff, diff_pow/curr_pow); + no_valid = false; + } +// cerr << "curr_pow: " << curr_pow << " diff_pow: " << diff_pow << " diff: " << diff_pow/curr_pow << endl; +// cerr << m_last_max_diff << endl; + } + if ((no_valid) || (m_last_max_diff>1)) + m_last_max_diff = 1; +// cerr << m_last_max_diff << endl; + } +} + +void Engine_Ext_SteadyState::Apply2Current() +{ + +} diff --git a/FDTD/extensions/engine_ext_steadystate.h b/FDTD/extensions/engine_ext_steadystate.h new file mode 100644 index 0000000..b66596e --- /dev/null +++ b/FDTD/extensions/engine_ext_steadystate.h @@ -0,0 +1,51 @@ +/* +* Copyright (C) 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 . +*/ + +#ifndef ENGINE_EXT_STEADYSTATE_H +#define ENGINE_EXT_STEADYSTATE_H + +#include "engine_extension.h" +#include "FDTD/engine.h" +#include "FDTD/operator.h" + +class Operator_Ext_SteadyState; +class Engine_Interface_FDTD; + +class Engine_Ext_SteadyState : public Engine_Extension +{ +public: + Engine_Ext_SteadyState(Operator_Ext_SteadyState* op_ext); + virtual ~Engine_Ext_SteadyState(); + + virtual void Apply2Voltages(); + virtual void Apply2Current(); + + void SetEngineInterface(Engine_Interface_FDTD* eng_if) {m_Eng_Interface=eng_if;} + double GetLastDiff() {return m_last_max_diff;} + +protected: + Operator_Ext_SteadyState* m_Op_SS; + double m_last_max_diff; + vector m_E_records; + vector m_H_records; + + double last_total_energy; + Engine_Interface_FDTD* m_Eng_Interface; +}; + + +#endif // ENGINE_EXT_STEADYSTATE_H diff --git a/FDTD/extensions/engine_extension.h b/FDTD/extensions/engine_extension.h index 01d3194..12ef712 100644 --- a/FDTD/extensions/engine_extension.h +++ b/FDTD/extensions/engine_extension.h @@ -21,6 +21,7 @@ #define ENG_EXT_PRIO_DEFAULT 0 //default engine extension priority // priority definitions for some important extensions +#define ENG_EXT_PRIO_STEADYSTATE +2e6 //steady state extension priority #define ENG_EXT_PRIO_UPML +1e6 //unaxial pml extension priority #define ENG_EXT_PRIO_CYLINDER +1e5 //cylindrial extension priority #define ENG_EXT_PRIO_TFSF +5e4 //total-field/scattered-field extension priority diff --git a/FDTD/extensions/operator_ext_steadystate.cpp b/FDTD/extensions/operator_ext_steadystate.cpp new file mode 100644 index 0000000..b8395a8 --- /dev/null +++ b/FDTD/extensions/operator_ext_steadystate.cpp @@ -0,0 +1,104 @@ +/* +* Copyright (C) 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 "operator_ext_steadystate.h" +#include "engine_ext_steadystate.h" + +Operator_Ext_SteadyState::Operator_Ext_SteadyState(Operator* op, double period): Operator_Extension(op) +{ + this->Reset(); + m_T_period = period; +} + +Operator_Ext_SteadyState::Operator_Ext_SteadyState(Operator* op, Operator_Ext_SteadyState* op_ext): Operator_Extension(op, op_ext) +{ + this->Reset(); + m_T_period = op_ext->m_T_period; +} + +Operator_Ext_SteadyState::~Operator_Ext_SteadyState() +{ +} + +Operator_Extension* Operator_Ext_SteadyState::Clone(Operator* op) +{ + //disable cloning, only the main operator needs to have a steady state detection + UNUSED(op); + return NULL; +} + +bool Operator_Ext_SteadyState::BuildExtension() +{ + double dT = m_Op->GetTimestep(); + m_TS_period = round(m_T_period/dT); + m_T_period = m_TS_period*dT; + return true; +} + +void Operator_Ext_SteadyState::Reset() +{ + for (int n=0;n<3;++n) + { + m_E_probe_pos[n].clear(); + m_H_probe_pos[n].clear(); + } + m_E_probe_dir.clear(); + m_H_probe_dir.clear(); + m_T_period = 0; + m_TS_period = 0; +} + +bool Operator_Ext_SteadyState::Add_E_Probe(unsigned int pos[3], int dir) +{ + if ((dir<0) || (dir>2)) + return false; + for (int n=0;n<3;++n) + if (pos[n]>=m_Op->GetNumberOfLines(n)) + return false; + for (int n=0;n<3;++n) + m_E_probe_pos[n].push_back(pos[n]); + m_E_probe_dir.push_back(dir); + return true; +} + +bool Operator_Ext_SteadyState::Add_H_Probe(unsigned int pos[3], int dir) +{ + if ((dir<0) || (dir>2)) + return false; + for (int n=0;n<3;++n) + if (pos[n]>=m_Op->GetNumberOfLines(n)) + return false; + for (int n=0;n<3;++n) + m_H_probe_pos[n].push_back(pos[n]); + m_H_probe_dir.push_back(dir); + return true; +} + +Engine_Extension* Operator_Ext_SteadyState::CreateEngineExtention() +{ + m_Eng_Ext = new Engine_Ext_SteadyState(this); + return m_Eng_Ext; +} + +void Operator_Ext_SteadyState::ShowStat(ostream &ostr) const +{ + Operator_Extension::ShowStat(ostr); + cout << "Period time (s): " << m_T_period << "\t Period TS: " << m_TS_period << endl; + cout << "Number of E probes\t: " << m_E_probe_dir.size() << endl; + cout << "Number of H probes\t: " << m_H_probe_dir.size() << endl; +} + diff --git a/FDTD/extensions/operator_ext_steadystate.h b/FDTD/extensions/operator_ext_steadystate.h new file mode 100644 index 0000000..b735ca3 --- /dev/null +++ b/FDTD/extensions/operator_ext_steadystate.h @@ -0,0 +1,61 @@ +/* +* Copyright (C) 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 . +*/ + +#ifndef OPERATOR_EXT_STEADYSTATE_H +#define OPERATOR_EXT_STEADYSTATE_H + +#include "operator_extension.h" +#include "FDTD/operator.h" + +class Engine_Ext_SteadyState; + +class Operator_Ext_SteadyState : public Operator_Extension +{ + friend class Engine_Ext_SteadyState; +public: + Operator_Ext_SteadyState(Operator* op, double period); + virtual ~Operator_Ext_SteadyState(); + + virtual Operator_Extension* Clone(Operator* op); + + virtual bool BuildExtension(); + virtual Engine_Extension* CreateEngineExtention(); + + virtual bool IsCylinderCoordsSave(bool closedAlpha, bool R0_included) const {UNUSED(closedAlpha); UNUSED(R0_included); return true;} + virtual bool IsCylindricalMultiGridSave(bool child) const {UNUSED(child); return true;} + virtual bool IsMPISave() const {return true;} + + virtual string GetExtensionName() const {return string("Steady-State Detection Extension");} + + virtual void ShowStat(ostream &ostr) const; + + virtual void Reset(); + + bool Add_E_Probe(unsigned int pos[3], int dir); + bool Add_H_Probe(unsigned int pos[3], int dir); + +protected: + Operator_Ext_SteadyState(Operator* op, Operator_Ext_SteadyState* op_ext); + double m_T_period; + unsigned int m_TS_period; + vector m_E_probe_pos[3]; + vector m_E_probe_dir; + vector m_H_probe_pos[3]; + vector m_H_probe_dir; +}; + +#endif // OPERATOR_EXT_STEADYSTATE_H diff --git a/FDTD/extensions/operator_extension.cpp b/FDTD/extensions/operator_extension.cpp index 000f938..27b2b93 100644 --- a/FDTD/extensions/operator_extension.cpp +++ b/FDTD/extensions/operator_extension.cpp @@ -28,6 +28,7 @@ Operator_Extension::Operator_Extension(Operator* op) m_Op_Cyl = dynamic_cast(op); if (m_Op_Cyl) m_CC_R0_included=m_Op_Cyl->GetR0Included(); + m_Eng_Ext = NULL; } Operator_Extension::~Operator_Extension() @@ -42,6 +43,7 @@ Operator_Extension::Operator_Extension(Operator* op, Operator_Extension* op_ext) m_Active = op_ext->m_Active; if (m_Op_Cyl) m_CC_R0_included=m_Op_Cyl->GetR0Included(); + m_Eng_Ext = NULL; } void Operator_Extension::ShowStat(ostream &ostr) const diff --git a/FDTD/extensions/operator_extension.h b/FDTD/extensions/operator_extension.h index da20658..eea7193 100644 --- a/FDTD/extensions/operator_extension.h +++ b/FDTD/extensions/operator_extension.h @@ -49,6 +49,7 @@ public: virtual bool BuildExtension() {return true;} virtual Engine_Extension* CreateEngineExtention() {return 0;} + virtual Engine_Extension* GetEngineExtention() {return m_Eng_Ext;} //! The cylindrical operator will check whether the extension is save to use. Default is false. Derive this method to override. virtual bool IsCylinderCoordsSave(bool closedAlpha, bool R0_included) const {UNUSED(closedAlpha); UNUSED(R0_included); return false;} @@ -78,6 +79,7 @@ protected: //FDTD Operator Operator* m_Op; + Engine_Extension* m_Eng_Ext; //Cylindrical FDTD Operator (not NULL if a cylindrical FDTD is used) Operator_Cylinder* m_Op_Cyl; diff --git a/FDTD/operator.cpp b/FDTD/operator.cpp index 42f9f74..ad932ba 100644 --- a/FDTD/operator.cpp +++ b/FDTD/operator.cpp @@ -999,6 +999,14 @@ int Operator::CalcECOperator( DebugFlags debugFlags ) dT*=m_TimeStepFactor; + if (m_Exc->GetSignalPeriod()>0) + { + unsigned int TS = floor(m_Exc->GetSignalPeriod()/dT); + double new_dT = m_Exc->GetSignalPeriod()/TS; + cerr << "Operartor::CalcECOperator: Increasing timestep by " << round((new_dT-dT)/dT*1000)/10.0 << "% to " << new_dT << " (" << dT << ") to match periodic signal" << endl; + dT = new_dT; + } + m_Exc->Reset(dT); InitOperator(); diff --git a/openEMS.pro b/openEMS.pro index d8e6e33..5981ee9 100644 --- a/openEMS.pro +++ b/openEMS.pro @@ -237,7 +237,9 @@ SOURCES += FDTD/extensions/engine_extension.cpp \ FDTD/extensions/operator_ext_excitation.cpp \ FDTD/extensions/engine_ext_excitation.cpp \ FDTD/extensions/operator_ext_tfsf.cpp \ - FDTD/extensions/engine_ext_tfsf.cpp + FDTD/extensions/engine_ext_tfsf.cpp \ + FDTD/extensions/operator_ext_steadystate.cpp \ + FDTD/extensions/engine_ext_steadystate.cpp # Common source files SOURCES += Common/operator_base.cpp \ @@ -303,7 +305,9 @@ HEADERS += FDTD/extensions/operator_extension.h \ FDTD/extensions/operator_ext_excitation.h \ FDTD/extensions/engine_ext_excitation.h \ FDTD/extensions/operator_ext_tfsf.h \ - FDTD/extensions/engine_ext_tfsf.h + FDTD/extensions/engine_ext_tfsf.h \ + FDTD/extensions/operator_ext_steadystate.h \ + FDTD/extensions/engine_ext_steadystate.h # Common header files HEADERS += Common/operator_base.h \ diff --git a/openems.cpp b/openems.cpp index b23f8a0..39a6ee2 100644 --- a/openems.cpp +++ b/openems.cpp @@ -31,6 +31,8 @@ #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" @@ -61,6 +63,7 @@ openEMS::openEMS() { FDTD_Op=NULL; FDTD_Eng=NULL; + Eng_Ext_SSD=NULL; m_CSX=NULL; PA=NULL; CylinderCoords = false; @@ -313,7 +316,7 @@ Engine_Interface_FDTD* openEMS::NewEngineInterface(int multigridlevel) if (mgl==multigridlevel) { if (g_settings.GetVerboseLevel()>0) - cerr << __func__ << ": Operator with requested multi-grid level found." << endl; + 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(); @@ -321,7 +324,7 @@ Engine_Interface_FDTD* openEMS::NewEngineInterface(int multigridlevel) if (op_cyl_mg==NULL) //inner most operator reached { if (g_settings.GetVerboseLevel()>0) - cerr << __func__ << ": Operator with highest multi-grid level chosen." << endl; + cout << __func__ << ": Operator with highest multi-grid level chosen." << endl; return new Engine_Interface_Cylindrical_FDTD(op_cyl_inner); } // try next level @@ -641,9 +644,6 @@ int openEMS::SetupFDTD(const char* file) if (m_OverSampling<2) m_OverSampling=2; - double maxTime=0; - FDTD_Opts->QueryDoubleAttribute("MaxTime",&maxTime); - TiXmlElement* BC = FDTD_Opts->FirstChildElement("BoundaryCond"); if (BC==NULL) { @@ -688,10 +688,12 @@ int openEMS::SetupFDTD(const char* file) { FDTD_Op->SetCellConstantMaterial(); if (g_settings.GetVerboseLevel()>0) - cerr << "Enabling constant cell material assumption." << endl; + cout << "Enabling constant cell material assumption." << endl; } m_Exc = new Excitation(); + if (!m_Exc->setupExcitation(FDTD_Opts->FirstChildElement("Excitation"))) + exit(2); FDTD_Op->SetExcitationSignal(m_Exc); FDTD_Op->AddExtension(new Operator_Ext_Excitation(FDTD_Op)); if (!CylinderCoords) @@ -713,6 +715,36 @@ int openEMS::SetupFDTD(const char* file) if (FDTD_Opts->QueryDoubleAttribute("TimeStepFactor",×tepfactor)==TIXML_SUCCESS) FDTD_Op->SetTimestepFactor(timestepfactor); + // 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] = 10; + 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] = FDTD_Op->GetNumberOfLines(n)-10; + 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) @@ -746,11 +778,13 @@ int openEMS::SetupFDTD(const char* file) FDTD_Op->SetMaterialStoreFlags(2,false); FDTD_Op->SetMaterialStoreFlags(3,false); + double maxTime=0; + FDTD_Opts->QueryDoubleAttribute("MaxTime",&maxTime); unsigned int maxTime_TS = (unsigned int)(maxTime/FDTD_Op->GetTimestep()); if ((maxTime_TS>0) && (maxTime_TSsetupExcitation(FDTD_Opts->FirstChildElement("Excitation"),NrTS)) + if (!m_Exc->buildExcitationSignal(NrTS)) exit(2); m_Exc->DumpVoltageExcite("et"); m_Exc->DumpCurrentExcite("ht"); @@ -784,6 +818,12 @@ int openEMS::SetupFDTD(const char* file) //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; @@ -879,7 +919,7 @@ void openEMS::RunFDTD() FDTD_Eng->IterateTS(step); step=PA->Process(); - if (ProcField->CheckTimestep()) + if ((Eng_Ext_SSD==NULL) && ProcField->CheckTimestep()) { currE = ProcField->CalcTotalEnergyEstimate(); if (currE>maxE) @@ -896,16 +936,24 @@ void openEMS::RunFDTD() if (t_diff>4) { - currE = ProcField->CalcTotalEnergyEstimate(); - if (currE>maxE) - maxE=currE; 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 (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; + 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; diff --git a/openems.h b/openems.h index 4072850..f1433d0 100644 --- a/openems.h +++ b/openems.h @@ -35,6 +35,7 @@ class TiXmlElement; class ContinuousStructure; class Engine_Interface_FDTD; class Excitation; +class Engine_Ext_SteadyState; double CalcDiffTime(timeval t1, timeval t2); string FormatTime(int sec); @@ -91,6 +92,7 @@ protected: bool m_CellConstantMaterial; Operator* FDTD_Op; Engine* FDTD_Eng; + Engine_Ext_SteadyState* Eng_Ext_SSD; ProcessingArray* PA; Excitation* m_Exc;