Skip to content
Snippets Groups Projects
Commit 73bc50ac authored by Boman Romain's avatar Boman Romain
Browse files

move fdtd to https://github.com/rboman/progs

parent af4a1d88
No related branches found
No related tags found
1 merge request!68Remove old modules (sph, Metafor, fdtd)
Pipeline #4173 passed
Hide whitespace changes
Inline Side-by-side
Showing
with 0 additions and 802 deletions
CMakeLists.txt
+ 0
− 5
View file @ 73bc50ac
......@@ -43,7 +43,6 @@ MESSAGE(STATUS "CMAKE_SYSTEM_NAME=\"${CMAKE_SYSTEM_NAME}\"")
MESSAGE(STATUS "CMAKE_CXX_COMPILER_ID=${CMAKE_CXX_COMPILER_ID}")
# (de)activate modules
OPTION(WAVES_USE_FDTD "Compile fdtd module" ON)
OPTION(WAVES_USE_FLOW "Compile flow module" ON)
OPTION(WAVES_USE_HEAT "Compile heat module" ON)
OPTION(WAVES_USE_MIRRORS "Compile mirrors module" ON)
......@@ -162,10 +161,6 @@ IF(WAVES_USE_TBOXVTK)
ADD_SUBDIRECTORY( tboxVtk )
ENDIF()
IF(WAVES_USE_FDTD)
ADD_SUBDIRECTORY( fdtd )
ENDIF()
IF(WAVES_USE_FLOW)
ADD_SUBDIRECTORY( flow )
ENDIF()
......
fdtd/CMakeLists.txt deleted 100644 → 0
+ 0
− 8
View file @ af4a1d88
# Add source dir
ADD_SUBDIRECTORY( src )
ADD_SUBDIRECTORY( _src )
# Add test dir
MACRO_AddTest(${CMAKE_CURRENT_SOURCE_DIR}/tests)
fdtd/__init__.py deleted 100644 → 0
+ 0
− 6
View file @ af4a1d88
# -*- coding: utf-8 -*-
# fdtd MODULE initialization file
import fwk
import tbox
from fdtdw import *
fdtd/_src/CMakeLists.txt deleted 100644 → 0
+ 0
− 49
View file @ af4a1d88
# Copyright 2020 University of Liège
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
# http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.
# CMake input file of the SWIG wrapper around "fdtdw.so"
INCLUDE(${SWIG_USE_FILE})
FILE(GLOB SRCS *.h *.cpp *.inl *.swg)
FILE(GLOB ISRCS *.i)
SET(CMAKE_SWIG_FLAGS "")
SET_SOURCE_FILES_PROPERTIES(${ISRCS} PROPERTIES CPLUSPLUS ON)
SET(SWINCFLAGS
-I${PROJECT_SOURCE_DIR}/fdtd/src
-I${PROJECT_SOURCE_DIR}/tbox/src
-I${PROJECT_SOURCE_DIR}/tbox/_src
-I${PROJECT_SOURCE_DIR}/fwk/src
-I${PROJECT_SOURCE_DIR}/fwk/_src
)
SET_SOURCE_FILES_PROPERTIES(${ISRCS} PROPERTIES SWIG_FLAGS "${SWINCFLAGS}")
if (${CMAKE_VERSION} VERSION_LESS "3.8.0")
SWIG_ADD_MODULE(fdtdw python ${ISRCS} ${SRCS})
else()
SWIG_ADD_LIBRARY(fdtdw LANGUAGE python SOURCES ${ISRCS} ${SRCS})
endif()
MACRO_DebugPostfix(_fdtdw)
TARGET_INCLUDE_DIRECTORIES(_fdtdw PRIVATE ${PROJECT_SOURCE_DIR}/fwk/_src
${PROJECT_SOURCE_DIR}/tbox/_src
${PYTHON_INCLUDE_PATH}
)
SWIG_LINK_LIBRARIES(fdtdw
fdtd tbox fwk ${PYTHON_LIBRARIES}
)
fdtd/_src/fdtdw.i deleted 100644 → 0
+ 0
− 60
View file @ af4a1d88
/*
* Copyright 2020 University of Liège
*
* Licensed under the Apache License, Version 2.0 (the "License");
* you may not use this file except in compliance with the License.
* You may obtain a copy of the License at
*
* http://www.apache.org/licenses/LICENSE-2.0
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*/
// SWIG input file of the 'fdtd' module
%feature("autodoc","1");
%module(docstring=
"'fdtdw' module: projet MP 2017/2018
(c) ULg - A&M",
directors="1",
threads="1"
) fdtdw
%{
#include <string>
#include <sstream>
#include <typeinfo>
#include "fdtd.h"
#include "wProblem.h"
#include "fwkw.h"
#include "tboxw.h"
%}
%include "fwkw.swg"
// ----------- MODULES UTILISES ------------
%import "tboxw.i"
// ----------- FDTD CLASSES ----------------
%include "fdtd.h"
%shared_ptr(fdtd::Problem);
%feature("director:except") {
if ($error != NULL) {
std::cout << "[in director:except]\n";
//throw Swig::DirectorMethodException();
throw std::runtime_error("Director problem");
}
}
%include "wProblem.h"
fdtd/src/CMakeLists.txt deleted 100644 → 0
+ 0
− 28
View file @ af4a1d88
# Copyright 2020 University of Liège
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
# http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.
# CMake input file of fdtd.so
FILE(GLOB SRCS *.h *.cpp *.inl *.hpp)
ADD_LIBRARY(fdtd SHARED ${SRCS})
MACRO_DebugPostfix(fdtd)
TARGET_INCLUDE_DIRECTORIES(fdtd PUBLIC ${PROJECT_SOURCE_DIR}/fdtd/src)
# ---
TARGET_LINK_LIBRARIES(fdtd tbox fwk)
SOURCE_GROUP(base REGULAR_EXPRESSION ".*\\.(cpp|inl|hpp|h)")
fdtd/src/fdtd.h deleted 100644 → 0
+ 0
− 43
View file @ af4a1d88
/*
* Copyright 2020 University of Liège
*
* Licensed under the Apache License, Version 2.0 (the "License");
* you may not use this file except in compliance with the License.
* You may obtain a copy of the License at
*
* http://www.apache.org/licenses/LICENSE-2.0
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*/
// global header of the "fdtd" module
#ifndef FDTD_H
#define FDTD_H
#if defined(WIN32)
#ifdef fdtd_EXPORTS
#define FDTD_API __declspec(dllexport)
#else
#define FDTD_API __declspec(dllimport)
#endif
#else
#define FDTD_API
#endif
#include "tbox.h"
/**
* @brief this namespace avoids conflicts with similar names in 'waves'
*/
namespace fdtd
{
class Problem;
}
#endif //FDTD_H
fdtd/src/wProblem.cpp deleted 100644 → 0
+ 0
− 32
View file @ af4a1d88
/*
* Copyright 2020 University of Liège
*
* Licensed under the Apache License, Version 2.0 (the "License");
* you may not use this file except in compliance with the License.
* You may obtain a copy of the License at
*
* http://www.apache.org/licenses/LICENSE-2.0
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*/
#include "wProblem.h"
using namespace fdtd;
Problem::Problem()
{
}
void Problem::write(std::ostream &out) const
{
out << "fdtd::Problem:\n";
}
Problem::~Problem()
{
std::cout << "~Problem()\n";
}
fdtd/src/wProblem.h deleted 100644 → 0
+ 0
− 45
View file @ af4a1d88
/*
* Copyright 2020 University of Liège
*
* Licensed under the Apache License, Version 2.0 (the "License");
* you may not use this file except in compliance with the License.
* You may obtain a copy of the License at
*
* http://www.apache.org/licenses/LICENSE-2.0
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*/
#ifndef WPROBLEM_H
#define WPROBLEM_H
#include "fdtd.h"
#include "wObject.h"
#include <memory>
#include <iostream>
#include <vector>
namespace fdtd
{
/**
* @brief a class containing the problem parameters
*/
class FDTD_API Problem : public fwk::wSharedObject
{
public:
Problem();
virtual ~Problem();
#ifndef SWIG
virtual void write(std::ostream &out) const override;
#endif
};
} // namespace fdtd
#endif //WPROBLEM_H
fdtd/tests/fdtd.py deleted 100755 → 0
+ 0
− 28
View file @ af4a1d88
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
# Copyright 2020 University of Liège
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
# http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.
import fdtd as fd
def main():
p = fd.Problem()
print(p)
if __name__ == "__main__":
main()
fdtd/tests/fdtd1d_1.py deleted 100755 → 0
+ 0
− 54
View file @ af4a1d88
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
# Copyright 2020 University of Liège
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
# http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.
import math
import numpy as np
from fwk.wutils import parseargs
args = parseargs()
if not args.nogui:
import matplotlib.pyplot as plt
SIZE=200
ez = np.zeros(SIZE, dtype=float)
hy = np.zeros(SIZE, dtype=float)
imp0 = 377.0
maxTime = 1000
result = np.zeros(maxTime, dtype=float)
for qTime in range(maxTime):
for mm in range(0,SIZE-1):
hy[mm] = hy[mm] + (ez[mm+1] - ez[mm]) / imp0
for mm in range(1,SIZE):
ez[mm] = ez[mm] + (hy[mm] - hy[mm-1]) * imp0
ez[0] = math.exp(-(qTime-30.)*(qTime-30.)/100.)
#print(ez[50])
result[qTime] = ez[50]
#
if not args.nogui:
plt.plot(result)
plt.title('Example p42')
plt.xlabel('time step')
plt.ylabel('Ez[50]')
plt.ylim(-1,1)
plt.grid(True)
plt.show()
fdtd/tests/fdtd1d_1b.py deleted 100755 → 0
+ 0
− 63
View file @ af4a1d88
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
# Copyright 2020 University of Liège
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
# http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.
import math
import numpy as np
from fwk.wutils import parseargs
args = parseargs()
if not args.nogui:
import matplotlib.pyplot as plt
SIZE=200
ez = np.zeros(SIZE, dtype=float)
hy = np.zeros(SIZE, dtype=float)
imp0 = 377.0
maxTime = 100
#plt.ion()
#
if not args.nogui:
line, = plt.plot(ez)
plt.title('Example p42')
plt.xlabel('time step')
plt.ylabel('Ez[50]')
plt.ylim(-1,1)
plt.grid(True)
plt.draw()
for qTime in range(maxTime):
for mm in range(0,SIZE-1):
hy[mm] = hy[mm] + (ez[mm+1] - ez[mm]) / imp0
for mm in range(1,SIZE):
ez[mm] = ez[mm] + (hy[mm] - hy[mm-1]) * imp0
ez[0] = math.exp(-(qTime-30.)*(qTime-30.)/100.)
#args = parseargs()
if not args.nogui:
line.set_ydata(ez)
plt.draw()
plt.pause(0.001)
#plt.draw()
#plt.ioff()
#plt.show()
#raw_input()
\ No newline at end of file
fdtd/tests/fdtd1d_1c.py deleted 100755 → 0
+ 0
− 92
View file @ af4a1d88
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
# Copyright 2020 University of Liège
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
# http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.
import math
import numpy as np
from fwk.wutils import parseargs
args = parseargs()
if not args.nogui:
import matplotlib.pyplot as plt
import matplotlib.animation as animation
SIZE = 200
ez = np.zeros(SIZE, dtype=float)
hy = np.zeros(SIZE, dtype=float)
imp0 = 377.0
maxTime = 300
nloop = 1
if not args.nogui:
fig, ax = plt.subplots()
line, = plt.plot(ez)
plt.title('Example p42')
plt.xlabel('X')
plt.ylabel('Ez')
#plt.ylim(-2,2)
plt.ylim(-1,1)
plt.grid(True)
plt.draw()
else:
line = None
qTime=0
def animate(i):
global qTime
for l in range(nloop):
qTime=qTime+1
# ----- Hy -------
# absorbing boundary condition (ABC)
hy[SIZE-1] = hy[SIZE-2]
# perfect magnetic conductor (PMC)
# hy[SIZE-1] = 0
for mm in range(0,SIZE-1):
hy[mm] = hy[mm] + (ez[mm+1] - ez[mm]) / imp0
# correction for TFSF boundary
hy[49] -= np.exp(-(qTime - 30.)*(qTime - 30.)/100.) / imp0
# ----- Ez -------
# absorbing boundary condition (ABC)
ez[0] = ez[1]
# prescribed source (or perfect electric conductor if =0 (PEC) )
#ez[0] = np.exp(-(qTime-30.)*(qTime-30.)/100.)
for mm in range(1,SIZE):
ez[mm] = ez[mm] + (hy[mm] - hy[mm-1]) * imp0
# additive source
#ez[50] += np.exp(-(qTime-30.)*(qTime-30.)/100.)
# correction for TFSF boundary
ez[50] += np.exp(-(qTime+0.5- (-0.5) - 30.)*(qTime+0.5- (-0.5) - 30.)/100.)
if not args.nogui:
line.set_ydata(ez)
ax.set_title('time = %d' % qTime)
return line,
if not args.nogui:
ani = animation.FuncAnimation(fig, animate, maxTime // nloop,
interval=25, repeat=False)
plt.show()
else:
# calls the loops without matplotlib
for i in range(maxTime // nloop):
animate(i)
\ No newline at end of file
fdtd/tests/fdtd1d_1d.py deleted 100755 → 0
+ 0
− 152
View file @ af4a1d88
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
# Copyright 2020 University of Liège
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
# http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.
import math
import numpy as np
from fwk.wutils import parseargs
args = parseargs()
if not args.nogui:
import matplotlib.pyplot as plt
import matplotlib.animation as animation
class FDTD:
def __init__(self):
self.imp0 = 377.0 # impedance
self.SIZE = 200 # nb of spatial grid points
self.maxTime = 450 # nb of time steps
self.nloop = 1 # nb of loops / refresh of the display
self.MAT_LAYER = 100
self.LOSS_LAYER = 180 # position of the lossy layer
self.LOSS = 0.02
self.EPSR = 9.0
self.TFSF = 50 # position of the TFSF boundary
self.Cdtds = 1.0 # Courant
self.width = 10. # width of the Gaussian input
self.delay = 30. # delay of the Gaussian input
def execute(self):
if not args.nogui:
self.fig, self.ax = plt.subplots() #nrows=2)
else:
self.fig = self.ax = None
self.qTime = 0
self.ez = np.zeros(self.SIZE)
self.hy = np.zeros(self.SIZE-1)
# E eq
self.ceze = np.zeros(self.SIZE)
self.cezh = np.zeros(self.SIZE)
self.ceze[:self.MAT_LAYER] = 1.0
self.ceze[self.MAT_LAYER:self.LOSS_LAYER] = 1.0
self.ceze[self.LOSS_LAYER:] = (1.0-self.LOSS) / (1.0+self.LOSS)
self.cezh[:self.MAT_LAYER] = self.imp0
self.cezh[self.MAT_LAYER:self.LOSS_LAYER] = self.imp0 / self.EPSR
self.cezh[self.LOSS_LAYER:] = self.imp0 / self.EPSR / (1.0+self.LOSS)
# H eq
self.chyh = np.zeros(self.SIZE-1)
self.chye = np.zeros(self.SIZE-1)
self.chyh[:self.LOSS_LAYER] = 1.0
self.chyh[self.LOSS_LAYER:] = (1.0-self.LOSS) / (1.0+self.LOSS)
self.chye[:self.LOSS_LAYER] = 1.0 / self.imp0
self.chye[self.LOSS_LAYER:] = 1.0 / self.imp0 / (1.0+self.LOSS)
# initialise matplotlib
xe = np.linspace(0.0, self.SIZE-1, self.SIZE)
xh = np.linspace(0.0, self.SIZE-2, self.SIZE-1) + 0.5
if not args.nogui:
self.line1, = plt.plot(xe, self.ez, label=r'$E_z$')
self.line2, = plt.plot(xh, self.hy*self.imp0, label=r'$H_y \eta_0$')
plt.title('---')
plt.xlabel(r'$X$')
plt.ylabel(r'$E_z & H_y\eta_0$')
#plt.ylim(-2,2)
plt.ylim(-1.1, 1.5)
plt.plot([self.MAT_LAYER, self.MAT_LAYER], plt.ylim(), 'r--', linewidth=2.0)
plt.plot([self.LOSS_LAYER, self.LOSS_LAYER], plt.ylim(), 'b--', linewidth=2.0)
plt.grid(True)
plt.legend(loc='upper left')
plt.draw()
ani = animation.FuncAnimation(self.fig, self.animate, (self.maxTime-1) // self.nloop,
interval=25, repeat=False)
plt.show()
else:
# calls the loops without matplotlib
for i in range((self.maxTime-1) // self.nloop):
self.animate(i)
def animate(self, i):
for l in range(self.nloop):
self.qTime = self.qTime + 1
self.updateH()
self.updateTFSF()
self.abc()
self.updateE()
# additive source
#ez[50] += np.exp(-(qTime-30.)*(qTime-30.)/100.)
if not args.nogui:
self.line1.set_ydata(self.ez)
self.line2.set_ydata(self.hy*self.imp0)
self.ax.set_title('time = %d' % self.qTime)
return self.line1, self.line2
def abc(self):
"""absorbing boundary condition"""
self.ez[0] = self.ez[1]
#self.ez[self.SIZE-1] = self.ez[self.SIZE-2]
def updateH(self):
for mm in range(0, self.SIZE-1):
self.hy[mm] = self.chyh[mm] * self.hy[mm] + self.chye[mm] * (self.ez[mm+1] - self.ez[mm])
def updateE(self):
for mm in range(1, self.SIZE-1):
self.ez[mm] = self.ceze[mm] * self.ez[mm] + self.cezh[mm] * (self.hy[mm] - self.hy[mm-1])
def ezInc(self, time, location):
return np.exp(-math.pow((time - self.delay - location/self.Cdtds) / self.width, 2))
def updateTFSF(self):
# correction for TFSF boundary
self.hy[self.TFSF] -= self.ezInc(self.qTime, 0.0) * self.chye[self.TFSF]
self.ez[self.TFSF+1] += self.ezInc(self.qTime+0.5, -0.5)
if __name__ == "__main__":
sim = FDTD()
sim.execute()
fdtd/tests/projethpc.py deleted 100755 → 0
+ 0
− 137
View file @ af4a1d88
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
# Copyright 2020 University of Liège
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
# http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.
import math
import numpy as np
from fwk.wutils import parseargs
args = parseargs()
if not args.nogui:
import matplotlib.pyplot as plt
import matplotlib.animation as animation
from matplotlib import cm
from mpl_toolkits.mplot3d import Axes3D
class FDTD:
def __init__(self):
self.L = 300.0e-3
self.Ls = 10.0e-3
self.f = 2.4e9 # Hz
self.mu = 4e-7*math.pi # H/m
self.eps = 8.854e-12 # F/m
self.c = 1 / math.sqrt(self.mu*self.eps)
print('c =', self.c)
print('wavelength =', self.c / self.f) # ~12cm pour 2.4GHz
self.dx = 5e-3
self.dt = self.dx/self.c/2
#self.dt = 1.0/self.f/20 #math.sqrt(3)*self.dx * 0.5 / self.c
print('dt =', self.dt)
self.nstep = 50
def run(self):
print('running...')
# calculate the grid size
npx = int(math.floor(self.L/2.0/self.dx))*2+1
print('npx =', npx)
nps = int(math.floor(self.Ls/2.0/self.dx))*2+1
print('nps =', nps)
is1 = (npx-nps) // 2
Z = math.sqrt(self.mu / self.eps)
print('Z =', Z)
CFL = self.c*self.dt/self.dx
print('CFL =', CFL)
# init fields
self.ez = np.zeros( (npx,npx) )
self.hx = np.zeros( (npx,npx) )
self.hy = np.zeros( (npx,npx) )
X = np.arange(0, npx, 1.0)
Y = np.arange(0, npx, 1.0)
X, Y = np.meshgrid(X, Y)
#print(X)
args = parseargs()
if not args.nogui:
fig = plt.figure()
ax = Axes3D(fig)
msh = ax.plot_surface(X, Y, self.ez, rstride=1, cstride=1, cmap=cm.viridis)
ax.set_zlim(-1, 1)
ax.set_xlabel(r'$X$')
ax.set_ylabel(r'$Y$')
ax.set_zlabel(r'$E_z$')
plt.draw()
# time integration
time = 0.0
for i in range(self.nstep):
time += self.dt
for i in range(npx):
for j in range(npx-1):
self.hx[i,j] += - CFL/Z * (self.ez[i,j+1] - self.ez[i,j])
for i in range(npx-1):
for j in range(npx):
self.hy[i,j] += CFL/Z * (self.ez[i+1,j] - self.ez[i,j])
for i in range(1,npx):
for j in range(1,npx):
self.ez[i,j] += Z*CFL * (self.hy[i,j] - self.hy[i-1,j]) - Z*CFL * (self.hx[i,j] - self.hx[i,j-1])
# apply source
for i in range(is1,is1+nps+1):
for j in range(is1,is1+nps+1):
self.ez[i,j] = math.sin(2*math.pi*self.f*time)
if not args.nogui:
if msh:
ax.collections.remove(msh)
msh = ax.plot_surface(X, Y, self.ez, rstride=1, cstride=1, cmap=cm.viridis)
plt.pause(.001)
if not args.nogui:
plt.show()
if __name__=="__main__":
"""
X = np.arange(-5, 5, 0.25)
Y = np.arange(-5, 5, 0.25)
X, Y = np.meshgrid(X, Y)
R = np.sqrt(X**2 + Y**2)
Z = np.sin(R)
fig = plt.figure()
ax = Axes3D(fig)
ax.plot_surface(X, Y, Z, rstride=1, cstride=1, cmap=cm.viridis)
plt.show()
"""
sim = FDTD()
sim.run()
\ No newline at end of file
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment