From 06bbec106f522212518f441c8d7eefe97b35c24c Mon Sep 17 00:00:00 2001 From: Thorsten Liebig Date: Mon, 4 Oct 2010 11:52:39 +0200 Subject: [PATCH] new extension: upml Operator extension implementing an uniaxial perfectly matched layer (upml) This new pml implementation is going to replace the old split-field pml --- FDTD/engine_ext_upml.cpp | 173 ++++++++++++++++++++ FDTD/engine_ext_upml.h | 48 ++++++ FDTD/operator.cpp | 47 +++--- FDTD/operator_ext_upml.cpp | 321 +++++++++++++++++++++++++++++++++++++ FDTD/operator_ext_upml.h | 88 ++++++++++ openEMS.pro | 8 +- openems.cpp | 7 +- 7 files changed, 669 insertions(+), 23 deletions(-) create mode 100644 FDTD/engine_ext_upml.cpp create mode 100644 FDTD/engine_ext_upml.h create mode 100644 FDTD/operator_ext_upml.cpp create mode 100644 FDTD/operator_ext_upml.h diff --git a/FDTD/engine_ext_upml.cpp b/FDTD/engine_ext_upml.cpp new file mode 100644 index 0000000..7b5193c --- /dev/null +++ b/FDTD/engine_ext_upml.cpp @@ -0,0 +1,173 @@ +/* +* 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 "engine_ext_upml.h" +#include "operator_ext_upml.h" +#include "engine.h" +#include "engine_sse.h" +#include "tools/array_ops.h" + +Engine_Ext_UPML::Engine_Ext_UPML(Operator_Ext_UPML* op_ext) : Engine_Extension(op_ext) +{ + m_Op_UPML = op_ext; + + //this ABC extension should be executed first! + m_Priority = 1e6; + + volt = Create_N_3DArray(m_Op_UPML->m_numLines); + volt_flux = Create_N_3DArray(m_Op_UPML->m_numLines); + curr = Create_N_3DArray(m_Op_UPML->m_numLines); + curr_flux = Create_N_3DArray(m_Op_UPML->m_numLines); +} + +Engine_Ext_UPML::~Engine_Ext_UPML() +{ + Delete_N_3DArray(volt,m_Op_UPML->m_numLines); + volt=NULL; + Delete_N_3DArray(volt_flux,m_Op_UPML->m_numLines); + volt_flux=NULL; + Delete_N_3DArray(curr,m_Op_UPML->m_numLines); + curr=NULL; + Delete_N_3DArray(curr_flux,m_Op_UPML->m_numLines); + curr_flux=NULL; +} + +void Engine_Ext_UPML::DoPreVoltageUpdates() +{ + if (m_Eng==NULL) + return; + + unsigned int pos[3]; + unsigned int loc_pos[3]; + + for (loc_pos[0]=0;loc_pos[0]m_numLines[0];++loc_pos[0]) + { + pos[0] = loc_pos[0] + m_Op_UPML->m_StartPos[0]; + for (loc_pos[1]=0;loc_pos[1]m_numLines[1];++loc_pos[1]) + { + pos[1] = loc_pos[1] + m_Op_UPML->m_StartPos[1]; + for (loc_pos[2]=0;loc_pos[2]m_numLines[2];++loc_pos[2]) + { + pos[2] = loc_pos[2] + m_Op_UPML->m_StartPos[2]; + + volt[0][loc_pos[0]][loc_pos[1]][loc_pos[2]] = m_Op_UPML->vv[0][loc_pos[0]][loc_pos[1]][loc_pos[2]] * m_Eng->GetVolt(0,pos) + - m_Op_UPML->vvfo[0][loc_pos[0]][loc_pos[1]][loc_pos[2]] * volt_flux[0][loc_pos[0]][loc_pos[1]][loc_pos[2]]; + volt[1][loc_pos[0]][loc_pos[1]][loc_pos[2]] = m_Op_UPML->vv[1][loc_pos[0]][loc_pos[1]][loc_pos[2]] * m_Eng->GetVolt(1,pos) + - m_Op_UPML->vvfo[1][loc_pos[0]][loc_pos[1]][loc_pos[2]] * volt_flux[1][loc_pos[0]][loc_pos[1]][loc_pos[2]]; + volt[2][loc_pos[0]][loc_pos[1]][loc_pos[2]] = m_Op_UPML->vv[2][loc_pos[0]][loc_pos[1]][loc_pos[2]] * m_Eng->GetVolt(2,pos) + - m_Op_UPML->vvfo[2][loc_pos[0]][loc_pos[1]][loc_pos[2]] * volt_flux[2][loc_pos[0]][loc_pos[1]][loc_pos[2]]; + + m_Eng->SetVolt(0,pos, volt_flux[0][loc_pos[0]][loc_pos[1]][loc_pos[2]]); + m_Eng->SetVolt(1,pos, volt_flux[1][loc_pos[0]][loc_pos[1]][loc_pos[2]]); + m_Eng->SetVolt(2,pos, volt_flux[2][loc_pos[0]][loc_pos[1]][loc_pos[2]]); + } + } + } +} + +void Engine_Ext_UPML::DoPostVoltageUpdates() +{ + if (m_Eng==NULL) + return; + + unsigned int pos[3]; + unsigned int loc_pos[3]; + + for (loc_pos[0]=0;loc_pos[0]m_numLines[0];++loc_pos[0]) + { + pos[0] = loc_pos[0] + m_Op_UPML->m_StartPos[0]; + for (loc_pos[1]=0;loc_pos[1]m_numLines[1];++loc_pos[1]) + { + pos[1] = loc_pos[1] + m_Op_UPML->m_StartPos[1]; + for (loc_pos[2]=0;loc_pos[2]m_numLines[2];++loc_pos[2]) + { + pos[2] = loc_pos[2] + m_Op_UPML->m_StartPos[2]; + + volt_flux[0][loc_pos[0]][loc_pos[1]][loc_pos[2]] = m_Eng->GetVolt(0,pos); + volt_flux[1][loc_pos[0]][loc_pos[1]][loc_pos[2]] = m_Eng->GetVolt(1,pos); + volt_flux[2][loc_pos[0]][loc_pos[1]][loc_pos[2]] = m_Eng->GetVolt(2,pos); + + m_Eng->SetVolt(0,pos, volt[0][loc_pos[0]][loc_pos[1]][loc_pos[2]] + m_Op_UPML->vvfn[0][loc_pos[0]][loc_pos[1]][loc_pos[2]] * volt_flux[0][loc_pos[0]][loc_pos[1]][loc_pos[2]]); + m_Eng->SetVolt(1,pos, volt[1][loc_pos[0]][loc_pos[1]][loc_pos[2]] + m_Op_UPML->vvfn[1][loc_pos[0]][loc_pos[1]][loc_pos[2]] * volt_flux[1][loc_pos[0]][loc_pos[1]][loc_pos[2]]); + m_Eng->SetVolt(2,pos, volt[2][loc_pos[0]][loc_pos[1]][loc_pos[2]] + m_Op_UPML->vvfn[2][loc_pos[0]][loc_pos[1]][loc_pos[2]] * volt_flux[2][loc_pos[0]][loc_pos[1]][loc_pos[2]]); + } + } + } +} + +void Engine_Ext_UPML::DoPreCurrentUpdates() +{ + if (m_Eng==NULL) + return; + + unsigned int pos[3]; + unsigned int loc_pos[3]; + + for (loc_pos[0]=0;loc_pos[0]m_numLines[0];++loc_pos[0]) + { + pos[0] = loc_pos[0] + m_Op_UPML->m_StartPos[0]; + for (loc_pos[1]=0;loc_pos[1]m_numLines[1];++loc_pos[1]) + { + pos[1] = loc_pos[1] + m_Op_UPML->m_StartPos[1]; + for (loc_pos[2]=0;loc_pos[2]m_numLines[2];++loc_pos[2]) + { + pos[2] = loc_pos[2] + m_Op_UPML->m_StartPos[2]; + + curr[0][loc_pos[0]][loc_pos[1]][loc_pos[2]] = m_Op_UPML->ii[0][loc_pos[0]][loc_pos[1]][loc_pos[2]] * m_Eng->GetCurr(0,pos) + - m_Op_UPML->iifo[0][loc_pos[0]][loc_pos[1]][loc_pos[2]] * curr_flux[0][loc_pos[0]][loc_pos[1]][loc_pos[2]]; + curr[1][loc_pos[0]][loc_pos[1]][loc_pos[2]] = m_Op_UPML->ii[1][loc_pos[0]][loc_pos[1]][loc_pos[2]] * m_Eng->GetCurr(1,pos) + - m_Op_UPML->iifo[1][loc_pos[0]][loc_pos[1]][loc_pos[2]] * curr_flux[1][loc_pos[0]][loc_pos[1]][loc_pos[2]]; + curr[2][loc_pos[0]][loc_pos[1]][loc_pos[2]] = m_Op_UPML->ii[2][loc_pos[0]][loc_pos[1]][loc_pos[2]] * m_Eng->GetCurr(2,pos) + - m_Op_UPML->iifo[2][loc_pos[0]][loc_pos[1]][loc_pos[2]] * curr_flux[2][loc_pos[0]][loc_pos[1]][loc_pos[2]]; + + m_Eng->SetCurr(0,pos, curr_flux[0][loc_pos[0]][loc_pos[1]][loc_pos[2]]); + m_Eng->SetCurr(1,pos, curr_flux[1][loc_pos[0]][loc_pos[1]][loc_pos[2]]); + m_Eng->SetCurr(2,pos, curr_flux[2][loc_pos[0]][loc_pos[1]][loc_pos[2]]); + } + } + } +} + +void Engine_Ext_UPML::DoPostCurrentUpdates() +{ + if (m_Eng==NULL) + return; + + unsigned int pos[3]; + unsigned int loc_pos[3]; + + for (loc_pos[0]=0;loc_pos[0]m_numLines[0];++loc_pos[0]) + { + pos[0] = loc_pos[0] + m_Op_UPML->m_StartPos[0]; + for (loc_pos[1]=0;loc_pos[1]m_numLines[1];++loc_pos[1]) + { + pos[1] = loc_pos[1] + m_Op_UPML->m_StartPos[1]; + for (loc_pos[2]=0;loc_pos[2]m_numLines[2];++loc_pos[2]) + { + pos[2] = loc_pos[2] + m_Op_UPML->m_StartPos[2]; + + curr_flux[0][loc_pos[0]][loc_pos[1]][loc_pos[2]] = m_Eng->GetCurr(0,pos); + curr_flux[1][loc_pos[0]][loc_pos[1]][loc_pos[2]] = m_Eng->GetCurr(1,pos); + curr_flux[2][loc_pos[0]][loc_pos[1]][loc_pos[2]] = m_Eng->GetCurr(2,pos); + + m_Eng->SetCurr(0,pos, curr[0][loc_pos[0]][loc_pos[1]][loc_pos[2]] + m_Op_UPML->iifn[0][loc_pos[0]][loc_pos[1]][loc_pos[2]] * curr_flux[0][loc_pos[0]][loc_pos[1]][loc_pos[2]]); + m_Eng->SetCurr(1,pos, curr[1][loc_pos[0]][loc_pos[1]][loc_pos[2]] + m_Op_UPML->iifn[1][loc_pos[0]][loc_pos[1]][loc_pos[2]] * curr_flux[1][loc_pos[0]][loc_pos[1]][loc_pos[2]]); + m_Eng->SetCurr(2,pos, curr[2][loc_pos[0]][loc_pos[1]][loc_pos[2]] + m_Op_UPML->iifn[2][loc_pos[0]][loc_pos[1]][loc_pos[2]] * curr_flux[2][loc_pos[0]][loc_pos[1]][loc_pos[2]]); + } + } + } +} diff --git a/FDTD/engine_ext_upml.h b/FDTD/engine_ext_upml.h new file mode 100644 index 0000000..ea0a278 --- /dev/null +++ b/FDTD/engine_ext_upml.h @@ -0,0 +1,48 @@ +/* +* 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 ENGINE_EXT_UPML_H +#define ENGINE_EXT_UPML_H + +#include "engine_extension.h" +#include "engine.h" +#include "operator.h" + +class Operator_Ext_UPML; + +class Engine_Ext_UPML : public Engine_Extension +{ +public: + Engine_Ext_UPML(Operator_Ext_UPML* op_ext); + virtual ~Engine_Ext_UPML(); + + virtual void DoPreVoltageUpdates(); + virtual void DoPostVoltageUpdates(); + + virtual void DoPreCurrentUpdates(); + virtual void DoPostCurrentUpdates(); + +protected: + Operator_Ext_UPML* m_Op_UPML; + + FDTD_FLOAT**** volt; + FDTD_FLOAT**** curr; + FDTD_FLOAT**** volt_flux; + FDTD_FLOAT**** curr_flux; +}; + +#endif // ENGINE_EXT_UPML_H diff --git a/FDTD/operator.cpp b/FDTD/operator.cpp index cd5a2c1..99459de 100644 --- a/FDTD/operator.cpp +++ b/FDTD/operator.cpp @@ -563,6 +563,25 @@ int Operator::CalcECOperator() } } + //Apply PEC to all boundary's + bool PEC[6]={1,1,1,1,1,1}; + //make an exception for BC == -1 + for (int n=0;n<6;++n) + if ((m_BC[n]==-1)) + PEC[n] = false; + ApplyElectricBC(PEC); + + CalcPEC(); + + bool PMC[6]; + for (int n=0;n<6;++n) + PMC[n] = m_BC[n]==1; + ApplyMagneticBC(PMC); + + InitExcitation(); + + if (CalcFieldExcitation()==false) return -1; + //all information available for extension... create now... for (size_t n=0;nBuildExtension(); @@ -576,25 +595,6 @@ int Operator::CalcECOperator() delete[] EC_R[n];EC_R[n]=NULL; } - //Apply PEC to all boundary's - bool PEC[6]={1,1,1,1,1,1}; - //exception for pml boundaries - for (int n=0;n<6;++n) - if ((m_BC[n]==3) || (m_BC[n]==-1)) - PEC[n] = false; - ApplyElectricBC(PEC); - - InitExcitation(); - - if (CalcFieldExcitation()==false) return -1; - - CalcPEC(); - - bool PMC[6]; - for (int n=0;n<6;++n) - PMC[n] = m_BC[n]==1; - ApplyMagneticBC(PMC); - return 0; } @@ -654,6 +654,15 @@ void Operator::ApplyMagneticBC(bool* dirs) GetIV(nP,pos[0],pos[1],pos[2]) *= (FDTD_FLOAT)!dirs[2*n+1]; GetII(nPP,pos[0],pos[1],pos[2]) *= (FDTD_FLOAT)!dirs[2*n+1]; GetIV(nPP,pos[0],pos[1],pos[2]) *= (FDTD_FLOAT)!dirs[2*n+1]; + + //the last current lines are outside the FDTD domain and cannot be iterated by the FDTD engine + pos[n]=numLines[n]-1; + GetII(n,pos[0],pos[1],pos[2]) = 0; + GetIV(n,pos[0],pos[1],pos[2]) = 0; + GetII(nP,pos[0],pos[1],pos[2]) = 0; + GetIV(nP,pos[0],pos[1],pos[2]) = 0; + GetII(nPP,pos[0],pos[1],pos[2]) = 0; + GetIV(nPP,pos[0],pos[1],pos[2]) = 0; } } } diff --git a/FDTD/operator_ext_upml.cpp b/FDTD/operator_ext_upml.cpp new file mode 100644 index 0000000..752f950 --- /dev/null +++ b/FDTD/operator_ext_upml.cpp @@ -0,0 +1,321 @@ +/* +* 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 "operator_ext_upml.h" +#include "operator_cylinder.h" +#include "engine_ext_upml.h" +#include "tools/array_ops.h" +#include "fparser.hh" + +Operator_Ext_UPML::Operator_Ext_UPML(Operator* op) : Operator_Extension(op) +{ + m_GradingFunction = new FunctionParser(); + //default grading function + SetGradingFunction(" -log(1e-6)*log(2.5)/(2*dl*pow(2.5,W/dl)-1) * pow(2.5, D/dl) / Z "); + + for (int n=0;n<6;++n) + { + m_BC[n]=0; + m_Size[n]=0; + } + for (int n=0;n<3;++n) + { + m_StartPos[n]=0; + m_numLines[n]=0; + } + + vv = NULL; + vvfo = NULL; + vvfn = NULL; + ii = NULL; + iifo = NULL; + iifn = NULL; +} + +Operator_Ext_UPML::~Operator_Ext_UPML() +{ + delete m_GradingFunction; + m_GradingFunction = NULL; + DeleteOp(); +} + +bool Operator_Ext_UPML::Create_UPML(Operator* op, int BC[6], unsigned int size[6], string gradFunc) +{ + Operator_Ext_UPML* op_ext_upml=NULL; + unsigned int start[3]={0 ,0 ,0}; + unsigned int stop[3] ={op->GetNumberOfLines(0)-1,op->GetNumberOfLines(1)-1,op->GetNumberOfLines(2)-1}; + + //create a pml in x-direction over the full width of yz-space + if (BC[0]==3) + { + op_ext_upml = new Operator_Ext_UPML(op); + op_ext_upml->SetGradingFunction(gradFunc); + start[0]=0; + stop[0] =size[0]; + op_ext_upml->SetBoundaryCondition(BC, size); + op_ext_upml->SetRange(start,stop); + op->AddExtension(op_ext_upml); + } + if (BC[1]==3) + { + op_ext_upml = new Operator_Ext_UPML(op); + op_ext_upml->SetGradingFunction(gradFunc); + start[0]=op->GetNumberOfLines(0)-1-size[1]; + stop[0] =op->GetNumberOfLines(0)-1; + op_ext_upml->SetBoundaryCondition(BC, size); + op_ext_upml->SetRange(start,stop); + op->AddExtension(op_ext_upml); + } + + //create a pml in y-direction over the xz-space (if a pml in x-direction already exists, skip that corner regions) + start[0]=(size[0]+1)*(BC[0]==3); + stop[0] =op->GetNumberOfLines(0)-1-(size[0]+1)*(BC[1]==3); + + if (BC[2]==3) + { + op_ext_upml = new Operator_Ext_UPML(op); + op_ext_upml->SetGradingFunction(gradFunc); + start[1]=0; + stop[1] =size[2]; + op_ext_upml->SetBoundaryCondition(BC, size); + op_ext_upml->SetRange(start,stop); + op->AddExtension(op_ext_upml); + } + if (BC[3]==3) + { + op_ext_upml = new Operator_Ext_UPML(op); + op_ext_upml->SetGradingFunction(gradFunc); + start[1]=op->GetNumberOfLines(1)-1-size[3]; + stop[1] =op->GetNumberOfLines(1)-1; + op_ext_upml->SetBoundaryCondition(BC, size); + op_ext_upml->SetRange(start,stop); + op->AddExtension(op_ext_upml); + } + + //create a pml in z-direction over the xy-space (if a pml in x- and/or y-direction already exists, skip that corner/edge regions) + start[1]=(size[2]+1)*(BC[2]==3); + stop[1] =op->GetNumberOfLines(1)-1-(size[3]+1)*(BC[3]==3); + + if (BC[4]==3) + { + op_ext_upml = new Operator_Ext_UPML(op); + op_ext_upml->SetGradingFunction(gradFunc); + start[2]=0; + stop[2] =size[4]; + op_ext_upml->SetBoundaryCondition(BC, size); + op_ext_upml->SetRange(start,stop); + op->AddExtension(op_ext_upml); + } + if (BC[5]==3) + { + op_ext_upml = new Operator_Ext_UPML(op); + op_ext_upml->SetGradingFunction(gradFunc); + start[2]=op->GetNumberOfLines(2)-1-size[5]; + stop[2] =op->GetNumberOfLines(2)-1; + op_ext_upml->SetBoundaryCondition(BC, size); + op_ext_upml->SetRange(start,stop); + op->AddExtension(op_ext_upml); + } + + return true; +} + + +void Operator_Ext_UPML::DeleteOp() +{ + Delete_N_3DArray(vv,m_numLines); + vv = NULL; + Delete_N_3DArray(vvfo,m_numLines); + vvfo = NULL; + Delete_N_3DArray(vvfn,m_numLines); + vvfn = NULL; + Delete_N_3DArray(ii,m_numLines); + ii = NULL; + Delete_N_3DArray(iifo,m_numLines); + iifo = NULL; + Delete_N_3DArray(iifn,m_numLines); + iifn = NULL; +} + + +bool Operator_Ext_UPML::SetGradingFunction(string func) +{ + if (func.empty()) + return true; + + m_GradFunc = func; + int res = m_GradingFunction->Parse(m_GradFunc.c_str(), "D,dl,W,Z,N"); + if(res < 0) return true; + + cerr << "Operator_Ext_UPML::SetGradingFunction: Warning, an error occured parsing the pml grading function (see below) ..." << endl; + cerr << func << "\n" << string(res, ' ') << "^\n" << m_GradingFunction->ErrorMsg() << "\n"; + return false; +} + +void Operator_Ext_UPML::CalcGradingKappa(int ny, unsigned int pos[3], double Zm[3], double kappa_v[3], double kappa_i[3]) +{ + double depth=0; + double width=0; + for (int n=0;n<3;++n) + { + if ((pos[n] <= m_Size[2*n]) && (m_BC[2*n]==3)) //lower pml in n-dir + { + width = (m_Op->GetDiscLine(n,m_Size[2*n]) - m_Op->GetDiscLine(n,0))*m_Op->GetGridDelta(); + depth = width - (m_Op->GetDiscLine(n,pos[n]) - m_Op->GetDiscLine(n,0))*m_Op->GetGridDelta(); + + if (n==ny) + depth-=m_Op->GetMeshDelta(n,pos)/2; + double vars[5] = {depth, width/m_Size[2*n], width, Zm[n], m_Size[2*n]}; + if (depth>0) + kappa_v[n] = m_GradingFunction->Eval(vars); + else + kappa_v[n]=0; + if (n==ny) + depth+=m_Op->GetMeshDelta(n,pos)/2; + + if (n!=ny) + depth-=m_Op->GetMeshDelta(n,pos)/2; + if (depth<0) + depth=0; + vars[0]=depth; + if (depth>0) + kappa_i[n] = m_GradingFunction->Eval(vars); + else + kappa_i[n] = 0; + } + else if ((pos[n] >= m_Op->GetNumberOfLines(n) -1 -m_Size[2*n+1]) && (m_BC[2*n+1]==3)) //upper pml in n-dir + { + width = (m_Op->GetDiscLine(n,m_Op->GetNumberOfLines(n)-1) - m_Op->GetDiscLine(n,m_Op->GetNumberOfLines(n)-m_Size[2*n+1]-1))*m_Op->GetGridDelta(); + depth = width - (m_Op->GetDiscLine(n,m_Op->GetNumberOfLines(n)-1) - m_Op->GetDiscLine(n,pos[n]))*m_Op->GetGridDelta(); + + if (n==ny) + depth+=m_Op->GetMeshDelta(n,pos)/2; + double vars[5] = {depth, width/(m_Size[2*n]), width, Zm[n], m_Size[2*n]}; + if (depth>0) + kappa_v[n] = m_GradingFunction->Eval(vars); + else + kappa_v[n]=0; + if (n==ny) + depth-=m_Op->GetMeshDelta(n,pos)/2; + + if (n!=ny) + depth+=m_Op->GetMeshDelta(n,pos)/2; + if (depth>width) + depth=0; + vars[0]=depth; + if (depth>0) + kappa_i[n] = m_GradingFunction->Eval(vars); + else + kappa_i[n]=0; + } + else + { + kappa_v[n] = 0; + kappa_i[n] = 0; + } + } +} + +bool Operator_Ext_UPML::BuildExtension() +{ + /*Calculate the upml coefficients as defined in: + Allen Taflove, computational electrodynamics - the FDTD method, third edition, chapter 7.8, pages 297-300 + - modified by Thorsten Liebig to match the equivalent circuit (EC) FDTD method + - kappa is used for conductivities (instead of sigma) + */ + if (m_Op==NULL) + return false; + + DeleteOp(); + vv = Create_N_3DArray(m_numLines); + vvfo = Create_N_3DArray(m_numLines); + vvfn = Create_N_3DArray(m_numLines); + ii = Create_N_3DArray(m_numLines); + iifo = Create_N_3DArray(m_numLines); + iifn = Create_N_3DArray(m_numLines); + + unsigned int pos[3]; + unsigned int loc_pos[3]; + int nP,nPP; + double kappa_v[3]={0,0,0}; + double kappa_i[3]={0,0,0}; + double eff_Mat[3][4]; + double Zm[3]; + double dT = m_Op->GetTimestep(); + + for (loc_pos[0]=0;loc_pos[0]Calc_EffMatPos(0,pos,eff_Mat[0]); + m_Op->Calc_EffMatPos(1,pos,eff_Mat[1]); + m_Op->Calc_EffMatPos(2,pos,eff_Mat[2]); + Zm[0] = sqrt(eff_Mat[0][2]/eff_Mat[0][0]); + Zm[1] = sqrt(eff_Mat[1][2]/eff_Mat[1][0]); + Zm[2] = sqrt(eff_Mat[2][2]/eff_Mat[2][0]); + for (int n=0;n<3;++n) + { + CalcGradingKappa(n, pos,Zm,kappa_v,kappa_i); + nP = (n+1)%3; + nPP = (n+2)%3; + //check if pos is on PEC + if ( (m_Op->GetVV(n,pos[0],pos[1],pos[2]) + m_Op->GetVI(n,pos[0],pos[1],pos[2])) != 0 ) + { + //modify the original operator to perform eq. (7.85) by the main engine (EC-FDTD: equation is multiplied by delta_n) + //the engine extension will replace the original voltages with the "voltage flux" (volt*eps0) prior to the voltage updates + //after the updates are done the extension will calculate the new voltages (see below) and place them back into the main field domain + m_Op->GetVV(n,pos[0],pos[1],pos[2]) = (2*__EPS0__ - kappa_v[nP]*dT) / (2*__EPS0__ + kappa_v[nP]*dT); + m_Op->GetVI(n,pos[0],pos[1],pos[2]) = (2*__EPS0__*dT) / (2*__EPS0__ + kappa_v[nP]*dT) * m_Op->GetEdgeLength(n,pos) / m_Op->GetEdgeArea(n,pos); + + //operators needed by eq. (7.88) to calculate new voltages from old voltages and old and new "voltage fluxes" + vv [n][loc_pos[0]][loc_pos[1]][loc_pos[2]] = (2*__EPS0__ - kappa_v[nPP]*dT) / (2*__EPS0__ + kappa_v[nPP]*dT); + vvfn[n][loc_pos[0]][loc_pos[1]][loc_pos[2]] = (2*__EPS0__ + kappa_v[n]*dT) / (2*__EPS0__ + kappa_v[nPP]*dT)/eff_Mat[n][0]; + vvfo[n][loc_pos[0]][loc_pos[1]][loc_pos[2]] = (2*__EPS0__ - kappa_v[n]*dT) / (2*__EPS0__ + kappa_v[nPP]*dT)/eff_Mat[n][0]; + } + + //check if pos is on PMC + if ( (m_Op->GetII(n,pos[0],pos[1],pos[2]) + m_Op->GetIV(n,pos[0],pos[1],pos[2])) != 0 ) + { + //modify the original operator to perform eq. (7.89) by the main engine (EC-FDTD: equation is multiplied by delta_n) + //the engine extension will replace the original currents with the "current flux" (curr*mu0) prior to the current updates + //after the updates are done the extension will calculate the new currents (see below) and place them back into the main field domain + m_Op->GetII(n,pos[0],pos[1],pos[2]) = (2*__EPS0__ - kappa_i[nP]*dT) / (2*__EPS0__ + kappa_i[nP]*dT); + m_Op->GetIV(n,pos[0],pos[1],pos[2]) = (2*__EPS0__*dT) / (2*__EPS0__ + kappa_i[nP]*dT) * m_Op->GetEdgeLength(n,pos,true) / m_Op->GetEdgeArea(n,pos,true); + + //operators needed by eq. (7.90) to calculate new currents from old currents and old and new "current fluxes" + ii [n][loc_pos[0]][loc_pos[1]][loc_pos[2]] = (2*__EPS0__ - kappa_i[nPP]*dT) / (2*__EPS0__ + kappa_i[nPP]*dT); + iifn[n][loc_pos[0]][loc_pos[1]][loc_pos[2]] = (2*__EPS0__ + kappa_i[n]*dT) / (2*__EPS0__ + kappa_i[nPP]*dT)/eff_Mat[n][2]; + iifo[n][loc_pos[0]][loc_pos[1]][loc_pos[2]] = (2*__EPS0__ - kappa_i[n]*dT) / (2*__EPS0__ + kappa_i[nPP]*dT)/eff_Mat[n][2]; + } + } + } + } + } + return true; +} + +Engine_Extension* Operator_Ext_UPML::CreateEngineExtention() +{ + Engine_Ext_UPML* eng_ext = new Engine_Ext_UPML(this); + return eng_ext; +} diff --git a/FDTD/operator_ext_upml.h b/FDTD/operator_ext_upml.h new file mode 100644 index 0000000..9d5ca0d --- /dev/null +++ b/FDTD/operator_ext_upml.h @@ -0,0 +1,88 @@ +/* +* 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 OPERATOR_EXT_UPML_H +#define OPERATOR_EXT_UPML_H + +#include "operator.h" +#include "operator_extension.h" + +class FunctionParser; + +//! Operator extension implemention an uniaxial perfectly matched layer (upml) +/* + The priority for this extension should be the highest of all extensions since this operator will use the main engine to perform vital parts in the upml implementation. + Therefore the voltages and currents as well as the operator are replaced during these update process. + This extension is propably incompatible with the most other extensions operating in the same regions. + */ +class Operator_Ext_UPML : public Operator_Extension +{ + friend class Engine_Ext_UPML; +public: + virtual ~Operator_Ext_UPML(); + + void SetBoundaryCondition(int* BCs, unsigned int size[6]) {for (int n=0;n<6;++n) {m_BC[n]=BCs[n];m_Size[n]=size[n];}} + + void SetRange(const unsigned int start[3], const unsigned int stop[3]) {for (int n=0;n<3;++n) {m_StartPos[n]=start[n];m_numLines[n]=stop[n]-start[n]+1;}} + + //! Set the grading function for the pml + /*! + Define the pml grading grading function. + Predefined variables in this grading function are: + D = depth in the pml in meter + dl = mesh delta inside the pml in meter + W = width (length) of the pml in meter + N = number of cells for the pml + Z = wave impedance at the current depth and position + example: SetGradingFunction("-log(1e-6)*log(2.5)/(2*dl*pow(2.5,W/dl)-1) * pow(2.5, D/dl) / Z"); + + An empty function string will be ignored. + */ + virtual bool SetGradingFunction(string func); + + virtual bool BuildExtension(); + + virtual Engine_Extension* CreateEngineExtention(); + + virtual string GetExtensionName() const {return string("Uniaxial PML Extension");} + + //! Create the UPML + static bool Create_UPML(Operator* op, int BC[6], unsigned int size[6], string gradFunc); + +protected: + Operator_Ext_UPML(Operator* op); + int m_BC[6]; + unsigned int m_Size[6]; + + unsigned int m_StartPos[3]; + unsigned int m_numLines[3]; + + string m_GradFunc; + FunctionParser* m_GradingFunction; + + void CalcGradingKappa(int ny, unsigned int pos[3], double Zm[3], double kappa_v[3], double kappa_i[3]); + + void DeleteOp(); + FDTD_FLOAT**** vv; //calc new voltage from old voltage + FDTD_FLOAT**** vvfo; //calc new voltage from old voltage flux + FDTD_FLOAT**** vvfn; //calc new voltage from new voltage flux + FDTD_FLOAT**** ii; //calc new current from old current + FDTD_FLOAT**** iifo; //calc new current from old current flux + FDTD_FLOAT**** iifn; //calc new current from new current flux +}; + +#endif // OPERATOR_EXT_UPML_H diff --git a/openEMS.pro b/openEMS.pro index ea109dd..76b8124 100644 --- a/openEMS.pro +++ b/openEMS.pro @@ -84,7 +84,9 @@ SOURCES += main.cpp \ FDTD/processmodematch.cpp \ FDTD/operator_cylindermultigrid.cpp \ FDTD/engine_cylindermultigrid.cpp \ - FDTD/engine_ext_cylindermultigrid.cpp + FDTD/engine_ext_cylindermultigrid.cpp \ + FDTD/operator_ext_upml.cpp \ + FDTD/engine_ext_upml.cpp HEADERS += tools/ErrorMsg.h \ tools/AdrOp.h \ tools/constants.h \ @@ -126,7 +128,9 @@ HEADERS += tools/ErrorMsg.h \ FDTD/operator_cylindermultigrid.h \ FDTD/engine_cylindermultigrid.h \ FDTD/engine_ext_cylindermultigrid.h \ - tools/aligned_allocator.h + tools/aligned_allocator.h \ + FDTD/operator_ext_upml.h \ + FDTD/engine_ext_upml.h QMAKE_CXXFLAGS_RELEASE = -O3 \ -g \ -march=native diff --git a/openems.cpp b/openems.cpp index 47a3734..994eda6 100644 --- a/openems.cpp +++ b/openems.cpp @@ -24,6 +24,7 @@ #include "FDTD/operator_multithread.h" #include "FDTD/operator_ext_mur_abc.h" #include "FDTD/operator_ext_pml_sf.h" +#include "FDTD/operator_ext_upml.h" #include "FDTD/processvoltage.h" #include "FDTD/processcurrent.h" #include "FDTD/process_efield.h" @@ -191,7 +192,7 @@ bool openEMS::SetupBoundaryConditions(TiXmlElement* BC) { int EC; //error code of tinyxml int bounds[6] = {0,0,0,0,0,0}; //default boundary cond. (PEC) - int pml_size[6] = {8,8,8,8,8,8}; //default pml size + unsigned int pml_size[6] = {8,8,8,8,8,8}; //default pml size string s_bc; const char* tmp = BC->Attribute("PML_Grading"); string pml_gradFunc; @@ -244,7 +245,9 @@ bool openEMS::SetupBoundaryConditions(TiXmlElement* BC) FDTD_Op->AddExtension(op_ext_mur); } } - Build_Split_Field_PML(FDTD_Op,bounds,pml_size,pml_gradFunc); + + //create the upml + Operator_Ext_UPML::Create_UPML(FDTD_Op,bounds,pml_size,pml_gradFunc); return true; }