From 6f06497dabb2839492cadafb37f04396d0b591b7 Mon Sep 17 00:00:00 2001 From: Thorsten Liebig Date: Mon, 28 Jun 2010 18:05:03 +0200 Subject: [PATCH] Added frequency domain probe support Signed-off-by: Thorsten Liebig --- FDTD/excitation.cpp | 14 ++----- FDTD/excitation.h | 1 - FDTD/processcurrent.cpp | 47 +++++++++++++++++++++-- FDTD/processcurrent.h | 5 +++ FDTD/processing.cpp | 85 ++++++++++++++++++++++++++++++++++++++--- FDTD/processing.h | 25 ++++++++++++ FDTD/processvoltage.cpp | 44 +++++++++++++++++++-- FDTD/processvoltage.h | 6 +++ openEMS.pro | 6 ++- openems.cpp | 4 ++ tools/useful.cpp | 31 +++++++++++++++ tools/useful.h | 24 ++++++++++++ 12 files changed, 266 insertions(+), 26 deletions(-) create mode 100644 tools/useful.cpp create mode 100644 tools/useful.h diff --git a/FDTD/excitation.cpp b/FDTD/excitation.cpp index 290e254..5fc15e7 100644 --- a/FDTD/excitation.cpp +++ b/FDTD/excitation.cpp @@ -16,6 +16,7 @@ */ #include "tools/array_ops.h" +#include "tools/useful.h" #include "fparser.hh" #include "tinyxml.h" #include "excitation.h" @@ -100,13 +101,6 @@ bool Excitation::setupExcitation( TiXmlElement* Excite, unsigned int maxTS ) return true; } -unsigned int Excitation::CalcNyquistNum(double fmax) const -{ - if (dT==0) return 1; - double T0 = 1/fmax; - return floor(T0/2/dT); -} - unsigned int Excitation::GetMaxExcitationTimestep() const { FDTD_FLOAT maxAmp=0; @@ -147,7 +141,7 @@ void Excitation::CalcGaussianPulsExcitation(double f0, double fc, int nTS) Signal_curr[n] = cos(2.0*PI*f0*(t-9.0/(2.0*PI*fc)))*exp(-1*pow(2.0*PI*fc*t/3.0-3,2)); } - SetNyquistNum( CalcNyquistNum(f0+fc) ); + SetNyquistNum( CalcNyquistNum(f0+fc,dT) ); } void Excitation::CalcDiracPulsExcitation() @@ -216,7 +210,7 @@ void Excitation::CalcCustomExcitation(double f0, int nTS, string signal) Signal_curr[n] = fParse.Eval(vars); } - SetNyquistNum( CalcNyquistNum(f0) ); + SetNyquistNum( CalcNyquistNum(f0,dT) ); } void Excitation::CalcSinusExcitation(double f0, int nTS) @@ -240,7 +234,7 @@ void Excitation::CalcSinusExcitation(double f0, int nTS) Signal_curr[n] = sin(2.0*PI*f0*t); } - SetNyquistNum( CalcNyquistNum(f0) ); + SetNyquistNum( CalcNyquistNum(f0,dT) ); } void Excitation::setupVoltageExcitation( vector const volt_vIndex[3], vector const& volt_vExcit, diff --git a/FDTD/excitation.h b/FDTD/excitation.h index f35e5ba..bc31df2 100644 --- a/FDTD/excitation.h +++ b/FDTD/excitation.h @@ -34,7 +34,6 @@ public: //! Get the excitation timestep with the (first) max amplitude virtual unsigned int GetMaxExcitationTimestep() const; - unsigned int CalcNyquistNum(double fmax) const; void setupVoltageExcitation( vector const volt_vIndex[3], vector const& volt_vExcit, vector const& volt_vDelay, vector const& volt_vDir ); diff --git a/FDTD/processcurrent.cpp b/FDTD/processcurrent.cpp index f0ea782..e0efb07 100644 --- a/FDTD/processcurrent.cpp +++ b/FDTD/processcurrent.cpp @@ -18,6 +18,7 @@ #include "tools/global.h" #include "processcurrent.h" #include +#include ProcessCurrent::ProcessCurrent(Operator* op, Engine* eng) : Processing(op, eng) { @@ -25,6 +26,7 @@ ProcessCurrent::ProcessCurrent(Operator* op, Engine* eng) : Processing(op, eng) ProcessCurrent::~ProcessCurrent() { + ProcessCurrent::FlushData(); } void ProcessCurrent::DefineStartStopCoord(double* dstart, double* dstop) @@ -44,6 +46,14 @@ void ProcessCurrent::DefineStartStopCoord(double* dstart, double* dstop) } } +void ProcessCurrent::Init() +{ + FD_currents.clear(); + for (size_t n=0;nnumTS << " i: " << current << endl; + // cerr << "ts: " << Eng->numTS << " i: " << current << endl; current*=m_weight; - v_current.push_back(current); - //current is sampled half a timestep later then the voltages - file << setprecision(m_precision) << (0.5 + (double)Eng->GetNumberOfTimesteps())*Op->GetTimestep() << "\t" << current << endl; + + if (ProcessInterval) + { + if (Eng->GetNumberOfTimesteps()%ProcessInterval==0) + { + v_current.push_back(current); + //current is sampled half a timestep later then the voltages + file << setprecision(m_precision) << (0.5 + (double)Eng->GetNumberOfTimesteps())*Op->GetTimestep() << "\t" << current << endl; + } + } + + if (m_FD_Interval) + { + if (Eng->GetNumberOfTimesteps()%m_FD_Interval==0) + { + double T = ((double)Eng->GetNumberOfTimesteps() + 0.5) * Op->GetTimestep(); + for (size_t n=0;n v_current; + + vector<_Complex double> FD_currents; + void WriteFDCurrents(); }; #endif // PROCESSCURRENT_H diff --git a/FDTD/processing.cpp b/FDTD/processing.cpp index 931fcfa..1856670 100644 --- a/FDTD/processing.cpp +++ b/FDTD/processing.cpp @@ -16,7 +16,10 @@ */ #include "tools/global.h" +#include "tools/useful.h" #include "processing.h" +#include +#include Processing::Processing(Operator* op, Engine* eng) { @@ -26,7 +29,10 @@ Processing::Processing(Operator* op, Engine* eng) m_PS_pos = 0; SetPrecision(12); ProcessInterval=0; + m_FD_SampleCount=0; + m_FD_Interval=0; m_weight=1; + m_Flush = false; } Processing::~Processing() @@ -53,21 +59,36 @@ bool Processing::CheckTimestep() { if (Eng->GetNumberOfTimesteps()%ProcessInterval==0) return true; } + + if (m_FD_Interval) + { + if (Eng->GetNumberOfTimesteps()%m_FD_Interval==0) return true; + } return false; } int Processing::GetNextInterval() const { if (Enabled==false) return -1; - unsigned int next=-1; + int next=INT_MAX; if (m_ProcessSteps.size()>m_PS_pos) { - next = m_ProcessSteps.at(m_PS_pos)-Eng->GetNumberOfTimesteps(); + next = (int)m_ProcessSteps.at(m_PS_pos)-(int)Eng->GetNumberOfTimesteps(); + } + if (ProcessInterval!=0) + { + int next_Interval = (int)ProcessInterval - (int)Eng->GetNumberOfTimesteps()%ProcessInterval; + if (next_IntervalGetNumberOfTimesteps()%m_FD_Interval; + if (next_IntervalGetNumberOfTimesteps()%ProcessInterval; - if (next_Interval steps) } } +void Processing::AddFrequency(double freq) +{ + unsigned int nyquistTS = CalcNyquistNum(freq,Op->GetTimestep()); + + if (nyquistTS == 0) + { + cerr << "Processing::AddFrequency: Requested frequency " << freq << " is too high for the current timestep used... skipping..." << endl; + return; + } + else if (nyquistTSExc->GetNyquistNum()) + { + cerr << "Processing::AddFrequency: Warning: Requested frequency " << freq << " is higher than maximum excited frequency..." << endl; + } + + if (m_FD_Interval==0) + m_FD_Interval = Op->Exc->GetNyquistNum(); + if (m_FD_Interval>nyquistTS) + m_FD_Interval = nyquistTS; + + m_FD_Samples.push_back(freq); +} + +void Processing::AddFrequency(vector freqs) +{ + for (size_t n=0;nSnapToMesh(dstart,start)==false) @@ -222,11 +273,33 @@ void Processing::DumpBox2File( string vtkfilenameprefix, bool dualMesh ) const file.close(); } +void Processing::Dump_FD_Data(vector<_Complex double> value, double factor, string filename) +{ + ofstream file; + file.open( filename.c_str() ); + if (!file.is_open()) + cerr << "Can't open file: " << filename << endl; + + for (size_t n=0;nFlushNext(); + } +} + void ProcessingArray::Reset() { for (size_t i=0;i #include +#include +#include #include "operator.h" #include "engine.h" @@ -29,6 +31,7 @@ public: Processing(Operator* op, Engine* eng); virtual ~Processing(); + virtual void Init() {}; virtual void Reset(); virtual void DefineStartStopCoord(double* dstart, double* dstop); @@ -38,6 +41,9 @@ public: void AddStep(unsigned int step); void AddSteps(vector steps); + void AddFrequency(double freq); + void AddFrequency(vector freqs); + bool CheckTimestep(); virtual int Process() {return GetNextInterval();} @@ -49,17 +55,24 @@ public: virtual void SetWeight(double weight) {m_weight=weight;} virtual double GetWeight() {return m_weight;} + //! Invoke this flag to flush all stored data to disk + virtual void FlushNext() {m_Flush = true;} + virtual void FlushData() {}; + //! Set the dump precision void SetPrecision(unsigned int val) {m_precision = val;} virtual void OpenFile( string outfile ); virtual void DumpBox2File( string vtkfilenameprefix, bool dualMesh = false ) const; //!< dump geometry to file + protected: Operator* Op; Engine* Eng; unsigned int m_precision; + bool m_Flush; + double m_weight; bool Enabled; @@ -70,6 +83,15 @@ protected: size_t m_PS_pos; //! current position in list of processing steps vector m_ProcessSteps; //! list of processing steps + //! Vector of frequency samples + vector m_FD_Samples; + //! Number of samples already processed + unsigned int m_FD_SampleCount; + //! Sampling interval needed for the FD_Samples + unsigned int m_FD_Interval; + + void Dump_FD_Data(vector<_Complex double> value, double factor, string filename); + //! define/store snapped start/stop coords as mesh index unsigned int start[3]; unsigned int stop[3]; @@ -94,6 +116,9 @@ public: void AddProcessing(Processing* proc); + //! Invoke this flag to flush all stored data to disk for all processings on next Process() + void FlushNext(); + void Reset(); //! Deletes all given processing's, can be helpful, but use carefull!!! diff --git a/FDTD/processvoltage.cpp b/FDTD/processvoltage.cpp index ccec3ae..a847f85 100644 --- a/FDTD/processvoltage.cpp +++ b/FDTD/processvoltage.cpp @@ -16,6 +16,7 @@ */ #include "processvoltage.h" +#include #include ProcessVoltage::ProcessVoltage(Operator* op, Engine* eng) : Processing(op, eng) @@ -24,16 +25,53 @@ ProcessVoltage::ProcessVoltage(Operator* op, Engine* eng) : Processing(op, eng) ProcessVoltage::~ProcessVoltage() { + ProcessVoltage::FlushData(); +} + +void ProcessVoltage::Init() +{ + FD_voltages.clear(); + for (size_t n=0;nGetNumberOfTimesteps()*Op->GetTimestep() << "\t" << voltage << endl; + + if (ProcessInterval) + { + if (Eng->GetNumberOfTimesteps()%ProcessInterval==0) + { + voltages.push_back(voltage); + file << setprecision(m_precision) << (double)Eng->GetNumberOfTimesteps()*Op->GetTimestep() << "\t" << voltage << endl; + } + } + + if (m_FD_Interval) + { + if (Eng->GetNumberOfTimesteps()%m_FD_Interval==0) + { + double T = (double)Eng->GetNumberOfTimesteps() * Op->GetTimestep(); + for (size_t n=0;n voltages; + + vector<_Complex double> FD_voltages; + void WriteFDVoltages(); }; #endif // PROCESSVOLTAGE_H diff --git a/openEMS.pro b/openEMS.pro index 933c5ea..3df9bad 100644 --- a/openEMS.pro +++ b/openEMS.pro @@ -56,7 +56,8 @@ SOURCES += main.cpp \ FDTD/operator_sse_compressed.cpp \ FDTD/engine_sse_compressed.cpp \ FDTD/operator_multithread.cpp \ - tools/global.cpp + tools/global.cpp \ + tools/useful.cpp HEADERS += tools/ErrorMsg.h \ tools/AdrOp.h \ tools/constants.h \ @@ -83,7 +84,8 @@ HEADERS += tools/ErrorMsg.h \ FDTD/operator_sse_compressed.h \ FDTD/engine_sse_compressed.h \ FDTD/operator_multithread.h \ - tools/global.h + tools/global.h \ + tools/useful.h QMAKE_CXXFLAGS_RELEASE = -O3 \ -g \ -march=native diff --git a/openems.cpp b/openems.cpp index c210794..1f0293c 100644 --- a/openems.cpp +++ b/openems.cpp @@ -330,8 +330,10 @@ int openEMS::SetupFDTD(const char* file) continue; } proc->SetProcessInterval(Nyquist/m_OverSampling); + proc->AddFrequency(pb->GetFDSamples()); proc->DefineStartStopCoord(start,stop); proc->SetWeight(pb->GetWeighting()); + proc->Init(); PA->AddProcessing(proc); prim->SetPrimitiveUsed(true); } @@ -467,6 +469,8 @@ void openEMS::RunFDTD() cout << " --- Energy: ~" << setw(6) << setprecision(2) << std::scientific << currE << " (decrement: " << setw(6) << setprecision(2) << std::fixed << fabs(10.0*log10(change)) << "dB)" << endl; prevTime=currTime; prevTS=currTS; + + PA->FlushNext(); } } diff --git a/tools/useful.cpp b/tools/useful.cpp new file mode 100644 index 0000000..8baac9b --- /dev/null +++ b/tools/useful.cpp @@ -0,0 +1,31 @@ +/* +* Copyright (C) 2010 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 "useful.h" +#include +#include +#include +#include + +unsigned int CalcNyquistNum(double fmax, double dT) +{ + if (fmax==0) return UINT_MAX; + if (dT==0) return 1; + double T0 = 1/fmax; + return floor(T0/2/dT); +} + diff --git a/tools/useful.h b/tools/useful.h new file mode 100644 index 0000000..a5a36a3 --- /dev/null +++ b/tools/useful.h @@ -0,0 +1,24 @@ +/* +* Copyright (C) 2010 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 USEFUL_H +#define USEFUL_H + +//! Calc the nyquist number of timesteps for a given frequency and timestep +unsigned int CalcNyquistNum(double fmax, double dT); + +#endif // USEFUL_H