Lumped RLC parallel & series implementation (openEMS) (#121)

pull/129/head
G. L 2023-11-18 13:23:15 +02:00 committed by GitHub
parent 5f36e7f3a2
commit ee3f2b7d80
No known key found for this signature in database
GPG Key ID: 4AEE18F83AFDEB23
8 changed files with 889 additions and 5 deletions

View File

@ -24,6 +24,8 @@ set(SOURCES
${CMAKE_CURRENT_SOURCE_DIR}/engine_ext_tfsf.cpp ${CMAKE_CURRENT_SOURCE_DIR}/engine_ext_tfsf.cpp
${CMAKE_CURRENT_SOURCE_DIR}/operator_ext_steadystate.cpp ${CMAKE_CURRENT_SOURCE_DIR}/operator_ext_steadystate.cpp
${CMAKE_CURRENT_SOURCE_DIR}/engine_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 PARENT_SCOPE
) )

View File

@ -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 <http://www.gnu.org/licenses/>.
*/
#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;
}

View File

@ -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 <http://www.gnu.org/licenses/>.
*/
#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

View File

@ -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 <http://www.gnu.org/licenses/>.
*/
#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<Operator_Ext_LumpedRLC*>(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<CSProperties*> cs_props;
int dir;
CSPropLumpedElement::LEtype lumpedType;
vector<uint> v_pos[3];
vector<int> v_dir;
vector<double> v_ilv;
vector<double> v_i2v;
vector<double> v_vv2;
vector<double> v_vj1;
vector<double> v_vj2;
vector<double> v_vvd;
vector<double> v_ib0;
vector<double> v_b1;
vector<double> 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<CSPropLumpedElement*>(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<CSPrimitives*> cs_RLC_prims = cs_RLC_props->GetAllPrimitives();
for (size_t boxIdx = 0 ; boxIdx < cs_RLC_prims.size() ; ++boxIdx)
{
CSPrimBox* cBox = dynamic_cast<CSPrimBox*>(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]<uiStop[dir_p1])
{
m_Op->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]<uiStop[dir_p2])
{
m_Op->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]<uiStop[dir_p1])
{
m_Op->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]<uiStop[dir_p2])
{
m_Op->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;
}

View File

@ -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 <http://www.gnu.org/licenses/>.
*/
#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_ */

View File

@ -379,7 +379,6 @@ int Operator::SnapLine2Mesh(const double* start, const double* stop, unsigned in
return ret; return ret;
} }
Grid_Path Operator::FindPath(double start[], double stop[]) Grid_Path Operator::FindPath(double start[], double stop[])
{ {
Grid_Path path; Grid_Path path;
@ -790,7 +789,7 @@ void Operator::DumpMaterial2File(string filename)
delete vtk_Writer; delete vtk_Writer;
} }
bool Operator::SetupCSXGrid(CSRectGrid* grid) bool Operator::SetupCSXGrid(CSRectGrid* grid)
{ {
for (int n=0; n<3; ++n) 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() bool Operator::Calc_LumpedElements()
{ {
vector<CSProperties*> props = CSX->GetPropertyByType(CSProperties::LUMPED_ELEMENT); vector<CSProperties*> props = CSX->GetPropertyByType(CSProperties::LUMPED_ELEMENT);
for (size_t i=0;i<props.size();++i) for (size_t i=0;i<props.size();++i)
{ {
CSPropLumpedElement* PLE = dynamic_cast<CSPropLumpedElement*>(props.at(i)); CSPropLumpedElement* PLE = dynamic_cast<CSPropLumpedElement*>(props.at(i));
if (PLE==NULL) if (PLE==NULL)
return false; //sanity check: this should never happen! return false; //sanity check: this should never happen!
vector<CSPrimitives*> prims = PLE->GetAllPrimitives(); vector<CSPrimitives*> prims = PLE->GetAllPrimitives();
for (size_t bn=0;bn<prims.size();++bn) for (size_t bn=0;bn<prims.size();++bn)
{ {
@ -1555,6 +1558,11 @@ bool Operator::Calc_LumpedElements()
if (R<0) if (R<0)
R = NAN; R = NAN;
// If this is not a parallel RC, skip this.
if (!(this->IsLEparRC(PLE)))
continue;
if ((std::isnan(R)) && (std::isnan(C))) if ((std::isnan(R)) && (std::isnan(C)))
{ {
cerr << "Operator::Calc_LumpedElements(): Warning: Lumped Element R or C not specified! skipping. " cerr << "Operator::Calc_LumpedElements(): Warning: Lumped Element R or C not specified! skipping. "
@ -1703,6 +1711,18 @@ bool Operator::Calc_LumpedElements()
return true; 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() void Operator::Init_EC()
{ {
for (int n=0; n<3; ++n) for (int n=0; n<3; ++n)

View File

@ -33,12 +33,16 @@ class Operator : public Operator_Base
{ {
friend class Engine; friend class Engine;
friend class Engine_Interface_FDTD; 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_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_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_PML_SF_Plane;
friend class Operator_Ext_Excitation; friend class Operator_Ext_Excitation;
friend class Operator_Ext_UPML; friend class Operator_Ext_UPML;
friend class Operator_Ext_Cylinder; 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: public:
enum DebugFlags {None=0,debugMaterial=1,debugOperator=2,debugPEC=4}; enum DebugFlags {None=0,debugMaterial=1,debugOperator=2,debugPEC=4};
@ -244,6 +248,9 @@ protected:
//! Calculate and setup lumped elements //! Calculate and setup lumped elements
virtual bool Calc_LumpedElements(); 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 //! Store the size of the applied boundary conditions
int m_BC_Size[6]; int m_BC_Size[6];

View File

@ -30,6 +30,7 @@
#include "FDTD/extensions/operator_ext_mur_abc.h" #include "FDTD/extensions/operator_ext_mur_abc.h"
#include "FDTD/extensions/operator_ext_upml.h" #include "FDTD/extensions/operator_ext_upml.h"
#include "FDTD/extensions/operator_ext_lorentzmaterial.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_conductingsheet.h"
#include "FDTD/extensions/operator_ext_steadystate.h" #include "FDTD/extensions/operator_ext_steadystate.h"
#include "FDTD/extensions/engine_ext_steadystate.h" #include "FDTD/extensions/engine_ext_steadystate.h"
@ -292,7 +293,7 @@ void openEMS::WelcomeScreen()
#endif #endif
cout << " ---------------------------------------------------------------------- " << endl; cout << " ---------------------------------------------------------------------- " << endl;
cout << " | openEMS " << bits << " -- version " GIT_VERSION << endl; cout << " | openEMS " << bits << " -- version " << GIT_VERSION << endl;
cout << " | (C) 2010-2023 Thorsten Liebig <thorsten.liebig@gmx.de> GPL license" << endl; cout << " | (C) 2010-2023 Thorsten Liebig <thorsten.liebig@gmx.de> GPL license" << endl;
cout << " ---------------------------------------------------------------------- " << endl; cout << " ---------------------------------------------------------------------- " << endl;
cout << openEMS::GetExtLibsInfo("\t") << endl; cout << openEMS::GetExtLibsInfo("\t") << endl;
@ -995,6 +996,9 @@ int openEMS::SetupFDTD()
FDTD_Op->AddExtension(new Operator_Ext_LorentzMaterial(FDTD_Op)); FDTD_Op->AddExtension(new Operator_Ext_LorentzMaterial(FDTD_Op));
if (m_CSX->GetQtyPropertyType(CSProperties::CONDUCTINGSHEET)>0) if (m_CSX->GetQtyPropertyType(CSProperties::CONDUCTINGSHEET)>0)
FDTD_Op->AddExtension(new Operator_Ext_ConductingSheet(FDTD_Op, m_Exc->GetMaxFreq())); 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... //check all properties to request material storage during operator creation...
SetupMaterialStorages(); SetupMaterialStorages();