1 /*
2 * Copyright (C) 2015 Thorsten Liebig (Thorsten.Liebig@gmx.de)
3 *
4 * This program is free software: you can redistribute it and/or modify
5 * it under the terms of the GNU General Public License as published by
6 * the Free Software Foundation, either version 3 of the License, or
7 * (at your option) any later version.
8 *
9 * This program is distributed in the hope that it will be useful,
10 * but WITHOUT ANY WARRANTY; without even the implied warranty of
11 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12 * GNU General Public License for more details.
13 *
14 * You should have received a copy of the GNU General Public License
15 * along with this program. If not, see <http://www.gnu.org/licenses/>.
16 */
17
18 #include "engine_ext_steadystate.h"
19 #include "operator_ext_steadystate.h"
20 #include "FDTD/engine_sse.h"
21 #include "FDTD/engine_interface_fdtd.h"
22
Engine_Ext_SteadyState(Operator_Ext_SteadyState * op_ext)23 Engine_Ext_SteadyState::Engine_Ext_SteadyState(Operator_Ext_SteadyState* op_ext): Engine_Extension(op_ext)
24 {
25 m_Op_SS = op_ext;
26 m_Priority = ENG_EXT_PRIO_STEADYSTATE;
27
28 for (size_t n=0;n<m_Op_SS->m_E_probe_dir.size();++n)
29 {
30 double* rec = new double[m_Op_SS->m_TS_period*2];
31 m_E_records.push_back(rec);
32 }
33 m_last_max_diff = 1;
34 last_total_energy = 0;
35 m_Eng_Interface = NULL;
36 }
37
~Engine_Ext_SteadyState()38 Engine_Ext_SteadyState::~Engine_Ext_SteadyState()
39 {
40 for (size_t n=0;n<m_E_records.size();++n)
41 {
42 delete[] m_E_records.at(n);
43 m_E_records.at(n) = NULL;
44 }
45 m_E_records.clear();
46 delete m_Eng_Interface;
47 m_Eng_Interface = NULL;
48 }
49
Apply2Voltages()50 void Engine_Ext_SteadyState::Apply2Voltages()
51 {
52 unsigned int p = m_Op_SS->m_TS_period;
53 unsigned int TS = m_Eng->GetNumberOfTimesteps();
54 unsigned int rel_pos = m_Eng->GetNumberOfTimesteps()%(2*p);
55 for (size_t n=0;n<m_E_records.size();++n)
56 m_E_records.at(n)[rel_pos] = m_Eng->GetVolt(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));
57 if ((TS%(m_Op_SS->m_TS_period)==0) && (TS>=2*p))
58 {
59 bool no_valid = true;
60 m_last_max_diff = 0;
61 double curr_total_energy = m_Eng_Interface->CalcFastEnergy();
62 if (last_total_energy>0)
63 {
64 m_last_max_diff = abs(curr_total_energy-last_total_energy)/last_total_energy;
65 no_valid = false;
66 }
67 //cerr << curr_total_energy << "/" << last_total_energy << "=" << abs(curr_total_energy-last_total_energy)/last_total_energy << endl;
68 last_total_energy = curr_total_energy;
69 unsigned int old_pos = 0;
70 unsigned int new_pos = p;
71 if (rel_pos<=p)
72 {
73 new_pos = 0;
74 old_pos = p;
75 }
76 //cerr << TS << "/" << rel_pos << ": one period complete, new_pos" << new_pos << " old pos: " << old_pos << endl;
77 double *curr_pow = new double[m_E_records.size()];
78 double *diff_pow = new double[m_E_records.size()];
79 double max_pow = 0;
80 for (size_t n=0;n<m_E_records.size();++n)
81 {
82 double *buf = m_E_records.at(n);
83 curr_pow[n] = 0;
84 diff_pow[n] = 0;
85 for (unsigned int nt=0;nt<p;++nt)
86 {
87 curr_pow[n] += buf[nt+new_pos]*buf[nt+new_pos];
88 diff_pow[n] += (buf[nt+old_pos]-buf[nt+new_pos])*(buf[nt+old_pos]-buf[nt+new_pos]);
89 }
90 max_pow = max(max_pow, curr_pow[n]);
91 }
92 for (size_t n=0;n<m_E_records.size();++n)
93 {
94 //cerr << "curr_pow: " << curr_pow[n] << " diff_pow: " << diff_pow[n] << " diff: " << diff_pow[n]/curr_pow[n] << endl;
95 if (curr_pow[n]>max_pow*1e-2)
96 {
97 m_last_max_diff = max(m_last_max_diff, diff_pow[n]/curr_pow[n]);
98 //cerr << m_last_max_diff << endl;
99 no_valid = false;
100 }
101 }
102 if ((no_valid) || (m_last_max_diff>1))
103 m_last_max_diff = 1;
104 delete[] curr_pow; curr_pow = NULL;
105 //cerr << m_last_max_diff << endl;
106 }
107 }
108
Apply2Current()109 void Engine_Ext_SteadyState::Apply2Current()
110 {
111
112 }
113