diff --git a/FDTD/extensions/CMakeLists.txt b/FDTD/extensions/CMakeLists.txt index 3566944..4346092 100644 --- a/FDTD/extensions/CMakeLists.txt +++ b/FDTD/extensions/CMakeLists.txt @@ -24,6 +24,8 @@ set(SOURCES ${CMAKE_CURRENT_SOURCE_DIR}/engine_ext_tfsf.cpp ${CMAKE_CURRENT_SOURCE_DIR}/operator_ext_steadystate.cpp ${CMAKE_CURRENT_SOURCE_DIR}/engine_ext_steadystate.cpp + ${CMAKE_CURRENT_SOURCE_DIR}/operator_ext_lumpedRLC.cpp + ${CMAKE_CURRENT_SOURCE_DIR}/engine_ext_lumpedRLC.cpp PARENT_SCOPE ) diff --git a/FDTD/extensions/engine_ext_lumpedRLC.cpp b/FDTD/extensions/engine_ext_lumpedRLC.cpp new file mode 100644 index 0000000..6eac857 --- /dev/null +++ b/FDTD/extensions/engine_ext_lumpedRLC.cpp @@ -0,0 +1,195 @@ +/* +* Additional +* Copyright (C) 2023 Gadi Lahav (gadi@rfwithcare.com) +* +* 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_lumpedRLC.h" +#include "operator_ext_lumpedRLC.h" + +#include "FDTD/engine_sse.h" + +Engine_Ext_LumpedRLC::Engine_Ext_LumpedRLC(Operator_Ext_LumpedRLC* op_ext_RLC) : Engine_Extension(op_ext_RLC) +{ + // Local pointer of the operator. + m_Op_Ext_RLC = op_ext_RLC; + + v_Vdn = new FDTD_FLOAT*[3]; + v_Jn = new FDTD_FLOAT*[3]; + + // No additional allocations are required if there are no actual lumped elements. + if (!(m_Op_Ext_RLC->RLC_count)) + return; + + // Initialize ADE containers for currents and voltages + v_Il = new FDTD_FLOAT[m_Op_Ext_RLC->RLC_count]; + + for (uint posIdx = 0 ; posIdx < m_Op_Ext_RLC->RLC_count ; ++posIdx) + v_Il[posIdx] = 0.0; + + for (uint k = 0 ; k < 3 ; k++) + { + v_Vdn[k] = new FDTD_FLOAT[m_Op_Ext_RLC->RLC_count]; + v_Jn[k] = new FDTD_FLOAT[m_Op_Ext_RLC->RLC_count]; + + for (uint posIdx = 0 ; posIdx < m_Op_Ext_RLC->RLC_count ; ++posIdx) + { + v_Jn[k][posIdx] = 0.0; + v_Vdn[k][posIdx] = 0.0;; + } + } + +} + +Engine_Ext_LumpedRLC::~Engine_Ext_LumpedRLC() +{ + // Only delete if values were allocated in the first place + if (m_Op_Ext_RLC->RLC_count) + { + delete[] v_Il; + + for (uint k = 0 ; k < 3 ; k++) + { + delete[] v_Vdn[k]; + delete[] v_Jn[k]; + } + } + + delete[] v_Vdn; + delete[] v_Jn; + + v_Il = NULL; + + v_Vdn = NULL; + v_Jn = NULL; + + m_Op_Ext_RLC = NULL; + + +} + +void Engine_Ext_LumpedRLC::DoPreVoltageUpdates() +{ + uint **pos = m_Op_Ext_RLC->v_RLC_pos; + int *dir = m_Op_Ext_RLC->v_RLC_dir; + + // Iterate Vd containers + FDTD_FLOAT *v_temp; + v_temp = v_Vdn[2]; + v_Vdn[2] = v_Vdn[1]; + v_Vdn[1] = v_Vdn[0]; + v_Vdn[0] = v_temp; + + // In pre-process, only update the parallel inductor current: + for (uint pIdx = 0 ; pIdx < m_Op_Ext_RLC->RLC_count ; pIdx++) + v_Il[pIdx] += (m_Op_Ext_RLC->v_RLC_i2v[pIdx])*(m_Op_Ext_RLC->v_RLC_ilv[pIdx])*v_Vdn[1][pIdx]; + + return; +} + +void Engine_Ext_LumpedRLC::Apply2Voltages() +{ + uint **pos = m_Op_Ext_RLC->v_RLC_pos; + int *dir = m_Op_Ext_RLC->v_RLC_dir; + + // Iterate J containers + FDTD_FLOAT *v_temp; + v_temp = v_Jn[2]; + v_Jn[2] = v_Jn[1]; + v_Jn[1] = v_Jn[0]; + v_Jn[0] = v_temp; + + + // Read engine calculated node voltage + switch (m_Eng->GetType()) + { + case Engine::BASIC: + { + for (uint pIdx = 0 ; pIdx < m_Op_Ext_RLC->RLC_count ; pIdx++) + v_Vdn[0][pIdx] = m_Eng->Engine::GetVolt(dir[pIdx],pos[0][pIdx],pos[1][pIdx],pos[2][pIdx]); + + break; + } + case Engine::SSE: + { + Engine_sse* eng_sse = (Engine_sse*)m_Eng; + for (uint pIdx = 0 ; pIdx < m_Op_Ext_RLC->RLC_count ; pIdx++) + v_Vdn[0][pIdx] = eng_sse->Engine_sse::GetVolt(dir[pIdx],pos[0][pIdx],pos[1][pIdx],pos[2][pIdx]); + + break; + } + default: + { + for (uint pIdx = 0 ; pIdx < m_Op_Ext_RLC->RLC_count ; pIdx++) + v_Vdn[0][pIdx] = m_Eng->GetVolt(dir[pIdx],pos[0][pIdx],pos[1][pIdx],pos[2][pIdx]);; + + break; + } + } + + // Post process: Calculate node voltage with respect to the lumped RLC auxilliary quantity, J + for (uint pIdx = 0 ; pIdx < m_Op_Ext_RLC->RLC_count ; pIdx++) + { + // Calculate updated node voltage, with series and parallel additions + v_Vdn[0][pIdx] = (m_Op_Ext_RLC->v_RLC_vvd[pIdx])*( + v_Vdn[0][pIdx] - v_Il[pIdx] // Addition for Parallel inductor + + + (m_Op_Ext_RLC->v_RLC_vv2[pIdx])*v_Vdn[2][pIdx] // Vd[n-2] addition + + + (m_Op_Ext_RLC->v_RLC_vj1[pIdx])*v_Jn[1][pIdx] // J[n-1] addition + + + (m_Op_Ext_RLC->v_RLC_vj2[pIdx])*v_Jn[2][pIdx]); // J[n-2] addition + + // Update J[0] + v_Jn[0][pIdx] = (m_Op_Ext_RLC->v_RLC_ib0[pIdx])*(v_Vdn[0][pIdx] - v_Vdn[2][pIdx]) + - + ((m_Op_Ext_RLC->v_RLC_b1[pIdx])*(m_Op_Ext_RLC->v_RLC_ib0[pIdx]))*v_Jn[1][pIdx] + - + ((m_Op_Ext_RLC->v_RLC_b2[pIdx])*(m_Op_Ext_RLC->v_RLC_ib0[pIdx]))*v_Jn[2][pIdx]; + } + + + // Update node voltage + switch (m_Eng->GetType()) + { + case Engine::BASIC: + { + for (uint pIdx = 0 ; pIdx < m_Op_Ext_RLC->RLC_count ; pIdx++) + m_Eng->Engine::SetVolt(dir[pIdx],pos[0][pIdx],pos[1][pIdx],pos[2][pIdx],v_Vdn[0][pIdx]); + + break; + } + case Engine::SSE: + { + Engine_sse* eng_sse = (Engine_sse*)m_Eng; + for (uint pIdx = 0 ; pIdx < m_Op_Ext_RLC->RLC_count ; pIdx++) + eng_sse->Engine_sse::SetVolt(dir[pIdx],pos[0][pIdx],pos[1][pIdx],pos[2][pIdx],v_Vdn[0][pIdx]); + + break; + } + default: + { + for (uint pIdx = 0 ; pIdx < m_Op_Ext_RLC->RLC_count ; pIdx++) + m_Eng->SetVolt(dir[pIdx],pos[0][pIdx],pos[1][pIdx],pos[2][pIdx],v_Vdn[0][pIdx]); + + break; + } + } + + + return; +} + + diff --git a/FDTD/extensions/engine_ext_lumpedRLC.h b/FDTD/extensions/engine_ext_lumpedRLC.h new file mode 100644 index 0000000..af44af4 --- /dev/null +++ b/FDTD/extensions/engine_ext_lumpedRLC.h @@ -0,0 +1,54 @@ +/* +* Copyright (C) 2023 Gadi Lahav (gadi@rfwithcare.com) +* +* 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_LUMPEDRLC_H +#define ENGINE_EXT_LUMPEDRLC_H + +#include "engine_extension.h" +#include "FDTD/engine.h" +#include "FDTD/operator.h" + +class Operator_Ext_LumpedRLC; + +class Engine_Ext_LumpedRLC : public Engine_Extension +{ + friend class Operator_Ext_LumpedRLC; + friend class Operator; + friend class ContinuousStructure; + +public: + + Engine_Ext_LumpedRLC(Operator_Ext_LumpedRLC *op_ext_RLC); + virtual ~Engine_Ext_LumpedRLC(); + + virtual void DoPreVoltageUpdates(); + virtual void Apply2Voltages(); + +protected: + Operator_Ext_LumpedRLC* m_Op_Ext_RLC; + + // Auxilliary containers + + // Array setup: volt_C_ADE[mesh_pos] + FDTD_FLOAT *v_Il; // Container for current on inductor- Parallel RLC + + FDTD_FLOAT **v_Vdn; // Container for nodal vd at [n],[n-1],[n-2] + FDTD_FLOAT **v_Jn; // Container for nodal J at [n],[n-1],[n-2] + +}; + +#endif // ENGINE_EXT_LORENTZMATERIAL_H diff --git a/FDTD/extensions/operator_ext_lumpedRLC.cpp b/FDTD/extensions/operator_ext_lumpedRLC.cpp new file mode 100644 index 0000000..c4a915a --- /dev/null +++ b/FDTD/extensions/operator_ext_lumpedRLC.cpp @@ -0,0 +1,514 @@ +/* +* Copyright (C) 2023 Gadi Lahav (gadi@rfwithcare.com) +* +* 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.h" +#include "tools/array_ops.h" +#include "tools/constants.h" +//#include "cond_sheet_parameter.h" +#include "tools/AdrOp.h" + +#include "operator_ext_lumpedRLC.h" +#include "engine_ext_lumpedRLC.h" + +#include "CSPrimBox.h" +#include "CSProperties.h" +#include "CSPropLumpedElement.h" + +#define COPY_V2A(V,A) std::copy(V.begin(),V.end(),A) + +Operator_Ext_LumpedRLC::Operator_Ext_LumpedRLC(Operator* op) : Operator_Extension(op) +{ + // Parallel circuit coefficients + v_RLC_ilv = NULL; + v_RLC_i2v = NULL; + + // Series circuit coefficients + v_RLC_vv2 = NULL; // Coefficient for [n-2] time of Vd update in Vd equation + v_RLC_vj1 = NULL; // Coefficient for [n-1] time of J update in Vd equation + v_RLC_vj2 = NULL; // Coefficient for [n-2] time of J update in Vd equation + v_RLC_vvd = NULL; // Coefficient to multiply all Vd in the Vd update equation + v_RLC_ib0 = NULL; // Inverse of beta_0 + v_RLC_b1 = NULL; // beta_1 + v_RLC_b2 = NULL; // beta_2 + + // Additional containers + v_RLC_dir = NULL; + v_RLC_pos = NULL; + + RLC_count = 0; +} + +Operator_Ext_LumpedRLC::Operator_Ext_LumpedRLC(Operator* op, Operator_Ext_LumpedRLC* op_ext) : Operator_Extension(op,op_ext) +{ + // Parallel circuit coefficients + v_RLC_ilv = NULL; + v_RLC_i2v = NULL; + + // Series circuit coefficients + v_RLC_vv2 = NULL; // Coefficient for [n-2] time of Vd update in Vd equation + v_RLC_vj1 = NULL; // Coefficient for [n-1] time of J update in Vd equation + v_RLC_vj2 = NULL; // Coefficient for [n-2] time of J update in Vd equation + v_RLC_vvd = NULL; // Coefficient to multiply all Vd in the Vd update equation + v_RLC_ib0 = NULL; // Inverse of beta_0 + v_RLC_b1 = NULL; // beta_1 + v_RLC_b2 = NULL; // beta_2 + + // Additional containers + v_RLC_dir = NULL; + v_RLC_pos = NULL; + + RLC_count = 0; +} + +Operator_Ext_LumpedRLC::~Operator_Ext_LumpedRLC() +{ + if (this->RLC_count) + { + // Parallel circuit coefficients + delete[] v_RLC_ilv; + delete[] v_RLC_i2v; + + // Series circuit coefficients + delete[] v_RLC_vv2; + delete[] v_RLC_vj1; + delete[] v_RLC_vj2; + delete[] v_RLC_vvd; + delete[] v_RLC_ib0; + delete[] v_RLC_b1; + delete[] v_RLC_b2; + + // Additional containers + delete[] v_RLC_dir; + + for (uint dIdx = 0 ; dIdx < 3 ; dIdx++) + delete[] v_RLC_pos[dIdx]; + + delete[] v_RLC_pos; + + } +} + +Operator_Extension* Operator_Ext_LumpedRLC::Clone(Operator* op) +{ + if (dynamic_cast(this)==NULL) + return NULL; + return new Operator_Ext_LumpedRLC(op, this); +} + +bool Operator_Ext_LumpedRLC::BuildExtension() +{ + double dT = m_Op->GetTimestep(); + + double fMax = m_Op->GetExcitationSignal()->GetCenterFreq() + + + m_Op->GetExcitationSignal()->GetCutOffFreq(); + + uint pos[] = {0,0,0}; + + vector cs_props; + + int dir; + CSPropLumpedElement::LEtype lumpedType; + + vector v_pos[3]; + + vector v_dir; + + vector v_ilv; + vector v_i2v; + + vector v_vv2; + vector v_vj1; + vector v_vj2; + vector v_vvd; + vector v_ib0; + vector v_b1; + vector v_b2; + + // Lumped RLC parameters + double R, L, C; + + // clear all vectors to initialize them + for (uint dIdx = 0 ; dIdx < 3 ; dIdx++) + v_pos[dIdx].clear(); + + v_dir.clear(); + + v_ilv.clear(); + v_i2v.clear(); + + v_vv2.clear(); + v_vj1.clear(); + v_vj2.clear(); + v_vvd.clear(); + v_ib0.clear(); + v_b1.clear(); + v_b2.clear(); + + // Obtain from CSX (continuous structure) all the lumped RLC properties + // Properties are material properties, not the objects themselves + cs_props = m_Op->CSX->GetPropertyByType(CSProperties::LUMPED_ELEMENT); + + // Iterate through various properties. In theory, there should be a property set per- + // primitive, as each lumped element should have it's own unique properties. + for(size_t n = 0 ; n < cs_props.size() ; ++n) + { + // Cast current property to lumped RLC property continuous structure properties + CSPropLumpedElement* cs_RLC_props = dynamic_cast(cs_props.at(n)); + if (cs_RLC_props==NULL) + return false; //sanity check: this should never happen! + + // Store direction and type + dir = cs_RLC_props->GetDirection(); + lumpedType = cs_RLC_props->GetLEtype(); + + // if (lumpedType == LEtype::INVALID + if (lumpedType == CSPropLumpedElement::INVALID) + { + cerr << "Operator_Ext_LumpedRLC::BuildExtension(): Warning: RLCtype is invalid! considering as parallel. " + << " ID: " << cs_RLC_props->GetID() << " @ Property: " << cs_RLC_props->GetName() << endl; + lumpedType = CSPropLumpedElement::PARALLEL; + } + + // Extract R, L and C from property class + C = cs_RLC_props->GetCapacity(); + if (C < 0) + C = NAN; + R = cs_RLC_props->GetResistance(); + if (R < 0) + R = NAN; + L = cs_RLC_props->GetInductance(); + if (L < 0) + L = NAN; + + // Check that this is a lumped RLC + if (!(this->IsLElumpedRLC(cs_RLC_props))) + continue; + + if ((dir < 0) || (dir > 2)) + { + cerr << "Operator_Ext_LumpedRLC::Calc_LumpedElements(): Warning: Lumped Element direction is invalid! skipping. " + << " ID: " << cs_RLC_props->GetID() << " @ Property: " << cs_RLC_props->GetName() << endl; + continue; + } + + // Initialize other two direction containers + int dir_p1 = (dir + 1) % 3; + int dir_p2 = (dir + 2) % 3; + + // Now iterate through primitive(s). I still think there should be only one per- + // material definition, but maybe I'm wrong... + vector cs_RLC_prims = cs_RLC_props->GetAllPrimitives(); + + for (size_t boxIdx = 0 ; boxIdx < cs_RLC_prims.size() ; ++boxIdx) + { + CSPrimBox* cBox = dynamic_cast(cs_RLC_prims.at(boxIdx)); + + if (cBox) + { + + // Get box start and stop positions + unsigned int uiStart[3], + uiStop[3]; + + + // snap to the native coordinate system + int Snap_Dimension = + m_Op->SnapBox2Mesh( + cBox->GetStartCoord()->GetCoords(m_Op->m_MeshType), // Start Coord + cBox->GetStopCoord()->GetCoords(m_Op->m_MeshType), // Stop Coord + uiStart, // Start Index + uiStop, // Stop Index + false, // Dual (doublet) Grid? + true); // Full mesh? + + // Verify that snapped dimension is correct + if (Snap_Dimension<=0) + { + if (Snap_Dimension>=-1) + cerr << "Operator_Ext_LumpedRLC::BuildExtension(): Warning: Lumped RLC snapping failed! Dimension is: " << Snap_Dimension << " skipping. " + << " ID: " << cs_RLC_prims.at(boxIdx)->GetID() << " @ Property: " << cs_RLC_props->GetName() << endl; + // Snap_Dimension == -2 means outside the simulation domain --> no special warning, but box probably marked as unused! + continue; + } + + // Verify that in the direction of the current propagation, the size isn't zero. + if (uiStart[dir]==uiStop[dir]) + { + cerr << "Operator_Ext_LumpedRLC::BuildExtension(): Warning: Lumped RLC with zero (snapped) length is invalid! skipping. " + << " ID: " << cs_RLC_prims.at(boxIdx)->GetID() << " @ Property: " << cs_RLC_props->GetName() << endl; + continue; + } + + // Calculate number of cells per-direction + unsigned int Ncells_0 = uiStop[dir] - uiStart[dir], + Ncells_1 = uiStop[dir_p1] - uiStart[dir_p1] + 1, + Ncells_2 = uiStop[dir_p2] - uiStart[dir_p2] + 1; + + // All cells in directions 1 and 2 are considered parallel connection + unsigned int Npar = Ncells_1*Ncells_2; + + // Separate elements such that individual elements can be calculated. + double dL = L*Npar/Ncells_0, + dR = R*Npar/Ncells_0, + dG = (1/R)*Ncells_0/Npar, + dC = C*Ncells_0/Npar; + + double ib0 = 2.0*dT*dC/(4.0*dL*dC + 2.0*dT*dR*dC + dT*dT), + b1 = (dT*dT - 4.0*dL*dC)/(dT*dC), + b2 = (4.0*dL*dC - 2.0*dT*dR*dC + dT*dT)/(2.0*dT*dC); + + // Special case: If this is a parallel resonant circuit, and there is no + // parallel resistor, use zero conductivity. May be risky when low-loss + // simulations are involved + if (lumpedType == CSPropLumpedElement::PARALLEL) + if (R == 0.0) + dG = 0.0; + + int iPos = 0; + + double Zmin,Zcd_min; + + // In the various positions, update the capacitors and "inverse" resistors + for (pos[dir] = uiStart[dir] ; pos[dir] < uiStop[dir] ; ++pos[dir]) + { + for (pos[dir_p1] = uiStart[dir_p1] ; pos[dir_p1] <= uiStop[dir_p1] ; ++pos[dir_p1]) + { + for (pos[dir_p2] = uiStart[dir_p2] ; pos[dir_p2] <= uiStop[dir_p2] ; ++pos[dir_p2]) + { + iPos = m_Op->MainOp->SetPos(pos[0],pos[1],pos[2]); + + + // Separate to two different cases. Parallel and series + switch (lumpedType) + { + case CSPropLumpedElement::PARALLEL: + // Update capacitor either way. + if (dC > 0) + m_Op->EC_C[dir][iPos] = dC; + else + // This case takes the "natural" capacitor into account. + dC = m_Op->EC_C[dir][iPos]; + + v_i2v.push_back((dT/dC)/(1.0 + dT*dG/(2.0*dC))); + + // Update conductivity + if (R >= 0) + m_Op->EC_G[dir][iPos] = dG; + + // Update coefficients with respect to the parallel inductance + if (L > 0) + v_ilv.push_back(dT/dL); + else + v_ilv.push_back(0.0); + + // Take into account the case that the "natural" capacitor is too small + // with respect to the inductor or the resistor, and add a warning. + if (dC == 0) + { + double Cd = m_Op->EC_C[dir][iPos]; + Zmin = max(dR,2*PI*fMax*dL); + Zcd_min = 1.0/(2.0*PI*fMax*Cd); + + // Check if the "parasitic" capcitance is not small enough + if (Zcd_min < LUMPED_RLC_Z_FACT*Zmin) + { + Cd = 1.0/(2*PI*fMax*Zmin*LUMPED_RLC_Z_FACT); + m_Op->EC_C[dir][iPos] = Cd; + } + } + + v_vv2.push_back(0.0); + v_vj1.push_back(0.0); + v_vj2.push_back(0.0); + v_vvd.push_back(1.0); + v_ib0.push_back(0.0); + v_b1.push_back(0.0); + v_b2.push_back(0.0); + + // Update with discrete component values of + m_Op->Calc_ECOperatorPos(dir,pos); + + v_dir.push_back(dir); + + break; + + case CSPropLumpedElement::SERIES: + m_Op->EC_G[dir][iPos] = 0.0; + + // is a series inductor, modeled separately. + FDTD_FLOAT Cd = m_Op->EC_C[dir][iPos]; + + // Calculate minimum impedance, at maximum frequency + Zmin = sqrt(pow(dR,2) + pow(2*PI*fMax*dL - 1.0/(dC*2*PI*fMax),2)); + Zcd_min = 1.0/(2.0*PI*fMax*Cd); + + // Check if the "parasitic" capcitance is not small enough + if (Zcd_min < LUMPED_RLC_Z_FACT*Zmin) + { + Cd = 1.0/(2*PI*fMax*Zmin*LUMPED_RLC_Z_FACT); + m_Op->EC_C[dir][iPos] = Cd; + } + + // No contribution from parallel inductor + v_ilv.push_back(0.0); + v_i2v.push_back(0.0); + + // Contributions from series resistor and inductor + v_vv2.push_back(0.5*dT*ib0/Cd); + v_vj1.push_back(0.5*dT*(b1*ib0 - 1.0)/Cd); + v_vj2.push_back(0.5*dT*b2*ib0/Cd); + v_vvd.push_back(1.0/(1.0 + 0.5*dT*ib0/Cd)); + v_ib0.push_back(ib0); + v_b1.push_back(b1); + v_b2.push_back(b2); + + m_Op->Calc_ECOperatorPos(dir,pos); + + v_dir.push_back(dir); + + break; + } + + // Store position and direction + for (uint dIdx = 0 ; dIdx < 3 ; ++dIdx) + v_pos[dIdx].push_back(pos[dIdx]); + + } + } + } + + // Build metallic caps + if (cs_RLC_props->GetCaps()) + for (pos[dir_p1] = uiStart[dir_p1] ; pos[dir_p1] <= uiStop[dir_p1] ; ++pos[dir_p1]) + { + for (pos[dir_p2] = uiStart[dir_p2] ; pos[dir_p2] <= uiStop[dir_p2] ; ++pos[dir_p2]) + { + pos[dir]=uiStart[dir]; + if (pos[dir_p1]SetVV(dir_p1,pos[0],pos[1],pos[2], 0 ); + m_Op->SetVI(dir_p1,pos[0],pos[1],pos[2], 0 ); + ++(m_Op->m_Nr_PEC[dir_p1]); + } + + if (pos[dir_p2]SetVV(dir_p2,pos[0],pos[1],pos[2], 0 ); + m_Op->SetVI(dir_p2,pos[0],pos[1],pos[2], 0 ); + ++(m_Op->m_Nr_PEC[dir_p2]); + } + + pos[dir]=uiStop[dir]; + if (pos[dir_p1]SetVV(dir_p1,pos[0],pos[1],pos[2], 0 ); + m_Op->SetVI(dir_p1,pos[0],pos[1],pos[2], 0 ); + ++(m_Op->m_Nr_PEC[dir_p1]); + } + + if (pos[dir_p2]SetVV(dir_p2,pos[0],pos[1],pos[2], 0 ); + m_Op->SetVI(dir_p2,pos[0],pos[1],pos[2], 0 ); + ++(m_Op->m_Nr_PEC[dir_p2]); + } + } + } + + + // Mark as used + cBox->SetPrimitiveUsed(true); + } + } + + } + + // Start data storage + RLC_count = v_dir.size(); + + // values + if (RLC_count) + { + // Allocate space to all variables + v_RLC_dir = new int[RLC_count]; + + // Parallel circuit coefficients + v_RLC_ilv = new FDTD_FLOAT[RLC_count]; + v_RLC_i2v = new FDTD_FLOAT[RLC_count]; + + // Series circuit coefficients + v_RLC_vv2 = new FDTD_FLOAT[RLC_count]; + v_RLC_vj1 = new FDTD_FLOAT[RLC_count]; + v_RLC_vj2 = new FDTD_FLOAT[RLC_count]; + v_RLC_vvd = new FDTD_FLOAT[RLC_count]; + v_RLC_ib0 = new FDTD_FLOAT[RLC_count]; + v_RLC_b1 = new FDTD_FLOAT[RLC_count]; + v_RLC_b2 = new FDTD_FLOAT[RLC_count]; + + v_RLC_pos = new uint*[3]; + for (uint dIdx = 0 ; dIdx < 3 ; ++dIdx) + v_RLC_pos[dIdx] = new uint[RLC_count]; + + // Copy all vectors to arrays + COPY_V2A(v_dir, v_RLC_dir); + + COPY_V2A(v_ilv, v_RLC_ilv); + COPY_V2A(v_i2v, v_RLC_i2v); + + COPY_V2A(v_vv2,v_RLC_vv2); + COPY_V2A(v_vj1,v_RLC_vj1); + COPY_V2A(v_vj2,v_RLC_vj2); + COPY_V2A(v_vvd,v_RLC_vvd); + COPY_V2A(v_ib0,v_RLC_ib0); + COPY_V2A(v_b1,v_RLC_b1); + COPY_V2A(v_b2,v_RLC_b2); + + for (uint dIdx = 0 ; dIdx < 3 ; ++dIdx) + COPY_V2A(v_pos[dIdx],v_RLC_pos[dIdx]); + } + + return true; +} + +Engine_Extension* Operator_Ext_LumpedRLC::CreateEngineExtention() +{ + Engine_Ext_LumpedRLC* eng_ext_RLC = new Engine_Ext_LumpedRLC(this); + return eng_ext_RLC; +} + +void Operator_Ext_LumpedRLC::ShowStat(ostream &ostr) const +{ + Operator_Extension::ShowStat(ostr); + string On_Off[2] = {"Off", "On"}; + + ostr << "Active cells\t\t: " << RLC_count << endl; +} + +bool Operator_Ext_LumpedRLC::IsLElumpedRLC(const CSPropLumpedElement* const p_prop) +{ + CSPropLumpedElement::LEtype lumpedType = p_prop->GetLEtype(); + + double L = p_prop->GetInductance(); + + bool isParallelRLC = (lumpedType == CSPropLumpedElement::PARALLEL) && (L > 0.0); + bool isSeriesRLC = lumpedType == CSPropLumpedElement::SERIES; + + // This needs to be something that isn't a parallel RC circuit to add data to this extension. + return isParallelRLC || isSeriesRLC; +} + diff --git a/FDTD/extensions/operator_ext_lumpedRLC.h b/FDTD/extensions/operator_ext_lumpedRLC.h new file mode 100644 index 0000000..cb7d8ff --- /dev/null +++ b/FDTD/extensions/operator_ext_lumpedRLC.h @@ -0,0 +1,88 @@ +/* +* Copyright (C) 2023 Gadi Lahav (gadi@rfwithcare.com) +* +* 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_LUMPEDRLC_H_ +#define OPERATOR_EXT_LUMPEDRLC_H_ + +#include "vector" + +#include "FDTD/operator.h" +#include "operator_extension.h" +#include "operator_ext_cylinder.h" + +#include "engine_ext_lumpedRLC.h" + +#define LUMPED_RLC_Z_FACT 20.0 + +class Operator_Ext_LumpedRLC : public Operator_Extension +{ + friend class Engine_Ext_LumpedRLC; +public: + Operator_Ext_LumpedRLC(Operator* op); + + virtual ~Operator_Ext_LumpedRLC(); + + 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("Series\\Parallel Lumped RLC load");} + + virtual void ShowStat(ostream &ostr) const; + + virtual bool IsLElumpedRLC(const CSPropLumpedElement* const p_prop); + +protected: + //! Copy constructor + Operator_Ext_LumpedRLC(Operator* op, Operator_Ext_LumpedRLC* op_ext); + + // ADE update coefficients, array setup: coeff[mesh_pos_index] + + // Parallel circuit coefficients + FDTD_FLOAT *v_RLC_ilv; + FDTD_FLOAT *v_RLC_i2v; + + // Series circuit coefficients + FDTD_FLOAT *v_RLC_vv2; // Coefficient for [n-2] time of Vd update in Vd equation + FDTD_FLOAT *v_RLC_vj1; // Coefficient for [n-1] time of J update in Vd equation + FDTD_FLOAT *v_RLC_vj2; // Coefficient for [n-2] time of J update in Vd equation + FDTD_FLOAT *v_RLC_vvd; // Coefficient to multiply all Vd in the Vd update equation + FDTD_FLOAT *v_RLC_ib0; // Inverse of beta_0 + FDTD_FLOAT *v_RLC_b1; // beta_1 + FDTD_FLOAT *v_RLC_b2; // beta_2 + + // Additional containers + int *v_RLC_dir; + uint **v_RLC_pos; + + // Vector length indicator + uint RLC_count; + + + + +}; + + + +#endif /* OPERATOR_EXT_LUMPEDRLC_H_ */ diff --git a/FDTD/operator.cpp b/FDTD/operator.cpp index d435c45..d3a9e49 100644 --- a/FDTD/operator.cpp +++ b/FDTD/operator.cpp @@ -379,7 +379,6 @@ int Operator::SnapLine2Mesh(const double* start, const double* stop, unsigned in return ret; } - Grid_Path Operator::FindPath(double start[], double stop[]) { Grid_Path path; @@ -790,7 +789,7 @@ void Operator::DumpMaterial2File(string filename) delete vtk_Writer; } - bool Operator::SetupCSXGrid(CSRectGrid* grid) +bool Operator::SetupCSXGrid(CSRectGrid* grid) { for (int n=0; n<3; ++n) { @@ -1536,11 +1535,15 @@ bool Operator::Calc_EffMatPos(int ny, const unsigned int* pos, double* EffMat, v bool Operator::Calc_LumpedElements() { vector props = CSX->GetPropertyByType(CSProperties::LUMPED_ELEMENT); + for (size_t i=0;i(props.at(i)); + if (PLE==NULL) return false; //sanity check: this should never happen! + vector prims = PLE->GetAllPrimitives(); for (size_t bn=0;bnIsLEparRC(PLE))) + continue; + + if ((std::isnan(R)) && (std::isnan(C))) { cerr << "Operator::Calc_LumpedElements(): Warning: Lumped Element R or C not specified! skipping. " @@ -1703,6 +1711,18 @@ bool Operator::Calc_LumpedElements() return true; } +bool Operator::IsLEparRC(const CSPropLumpedElement* const p_prop) +{ + CSPropLumpedElement::LEtype lumpedType = p_prop->GetLEtype(); + + double L = p_prop->GetInductance(); + + bool IsParallelRC = (lumpedType == CSPropLumpedElement::PARALLEL) && !(L > 0.0); + + // This needs to be something that isn't a parallel RC circuit to add data to this extension. + return IsParallelRC; +} + void Operator::Init_EC() { for (int n=0; n<3; ++n) diff --git a/FDTD/operator.h b/FDTD/operator.h index 233ecf6..7edccd0 100644 --- a/FDTD/operator.h +++ b/FDTD/operator.h @@ -33,12 +33,16 @@ class Operator : public Operator_Base { friend class Engine; friend class Engine_Interface_FDTD; - friend class Operator_Ext_LorentzMaterial; //we need to find a way around this... friend class Operator_Extension only would be nice - friend class Operator_Ext_ConductingSheet; //we need to find a way around this... friend class Operator_Extension only would be nice + friend class Operator_Ext_LorentzMaterial; // We need to find a way around this... friend class Operator_Extension only would be nice + friend class Operator_Ext_ConductingSheet; // We need to find a way around this... friend class Operator_Extension only would be nice friend class Operator_Ext_PML_SF_Plane; friend class Operator_Ext_Excitation; friend class Operator_Ext_UPML; friend class Operator_Ext_Cylinder; + friend class Operator_Ext_LumpedRLC; // Gadi: I now know why the two previous remarks are here. + + // So apparaently I have to use functionality from operator + // in my "lumpedRLC" class. This is ugly... public: enum DebugFlags {None=0,debugMaterial=1,debugOperator=2,debugPEC=4}; @@ -244,6 +248,9 @@ protected: //! Calculate and setup lumped elements virtual bool Calc_LumpedElements(); + //! Condition to determine if this is a lumped RC, to invoke Calc_LumpedElements + virtual bool IsLEparRC(const CSPropLumpedElement* const p_prop); + //! Store the size of the applied boundary conditions int m_BC_Size[6]; diff --git a/openems.cpp b/openems.cpp index affb25c..120d250 100644 --- a/openems.cpp +++ b/openems.cpp @@ -30,6 +30,7 @@ #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_lumpedRLC.h" #include "FDTD/extensions/operator_ext_conductingsheet.h" #include "FDTD/extensions/operator_ext_steadystate.h" #include "FDTD/extensions/engine_ext_steadystate.h" @@ -292,7 +293,7 @@ void openEMS::WelcomeScreen() #endif cout << " ---------------------------------------------------------------------- " << endl; - cout << " | openEMS " << bits << " -- version " GIT_VERSION << endl; + cout << " | openEMS " << bits << " -- version " << GIT_VERSION << endl; cout << " | (C) 2010-2023 Thorsten Liebig GPL license" << endl; cout << " ---------------------------------------------------------------------- " << endl; cout << openEMS::GetExtLibsInfo("\t") << endl; @@ -995,6 +996,9 @@ int openEMS::SetupFDTD() 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())); + if (m_CSX->GetQtyPropertyType(CSProperties::LUMPED_ELEMENT)>0) + FDTD_Op->AddExtension(new Operator_Ext_LumpedRLC(FDTD_Op)); + //check all properties to request material storage during operator creation... SetupMaterialStorages();