Created
June 25, 2010 09:43
-
-
Save mizutomo/452649 to your computer and use it in GitHub Desktop.
1次元FDTD Erlang Parallel版
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
%%% File: Fdtd.erl | |
%%% Description: FDTD 1-Dimentional Code | |
-module(fdtd). | |
%-export([main/0]). | |
-compile(export_all). | |
%% Physical Constant | |
-record(physics, {pi, light, permittivity, permeability}). | |
%% Analysis Area Info | |
-record(area, {px, nx, dx}). | |
%% Time Info. | |
-record(time, {cfl, stop_time, dt, last_step}). | |
calc_const() -> | |
PI = 3.14159, | |
LIGHT = 2.99792458e8, | |
PERMITTIVITY = 8.85418782e-12, | |
PERMEABILITY = 4.0 * PI * 1.0e-7, | |
Physics = #physics{pi = PI, light = LIGHT, | |
permittivity = PERMITTIVITY, | |
permeability = PERMEABILITY}, | |
PX = 1.0, | |
NX = 1000, | |
DX = PX / NX, | |
Area = #area{px = PX, nx = NX, dx = DX}, | |
CFL = 0.5, | |
STOP_TIME = 30.0e-9, | |
DT = CFL / (LIGHT * (1.0 / DX)), | |
LAST_STEP = erlang:trunc(STOP_TIME / DT), | |
Time = #time{cfl = CFL, stop_time = STOP_TIME, dt = DT, last_step = LAST_STEP}, | |
{Physics, Area, Time}. | |
print_analysis_info({_, Area, Time}) -> | |
io:format("Area X: ~f [m]~n", [Area#area.px]), | |
io:format("Grid X: ~B~n", [Area#area.nx]), | |
io:format("Delta X: ~f [m]~n", [Area#area.dx]), | |
io:format("CFL : ~f\n", [Time#time.cfl]), | |
io:format("dt : ~e [sec]~n", [Time#time.dt]), | |
io:format("Step : ~B~n", [Time#time.last_step]). | |
create_list(N) -> | |
[0.0 || _X <- lists:seq(1, N)]. | |
set_soft_source(0, [H | T], Value) -> | |
[H + Value | T]; | |
set_soft_source(Index, [H | T], Value) -> | |
[H | set_soft_source(Index - 1, T, Value)]. | |
calc_stimulus(Step, Dt) -> | |
Tc = 10.0e-9, | |
Pw = 0.1e-9, | |
T_current = (erlang:float(Step) - 0.5) * Dt, | |
T_eff = (T_current - Tc) / Pw, | |
1.0 * math:exp(-1.0 * (T_eff * T_eff)). | |
create_diff_list(L) -> | |
[_ | T] = L, | |
[_ | T1] = lists:reverse(L), | |
NewL = lists:reverse(T1), | |
[X - Y || {X, Y} <- lists:zip(T, NewL)]. | |
cut_both_end_element([_|T]) -> | |
[_ | T1] = lists:reverse(T), | |
lists:reverse(T1). | |
calc_ey(Const, Ey, Hz) -> | |
Diff_Hz = create_diff_list(Hz), | |
Cut_Ey = cut_both_end_element(Ey), | |
[send_calc_ey_message({E, DH}, Const) || {E, DH} <- lists:zip(Cut_Ey, Diff_Hz)]. | |
send_calc_ey_message(EH, Const) -> | |
Pid = spawn(fdtd, calc_ey_main, []), | |
Pid ! {self(), {EH, Const}}, | |
receive | |
{Pid, Msg} -> Msg | |
end. | |
calc_ey_main() -> | |
receive | |
{From, {{E, DH}, {Eps, Dt, Dx}}} -> | |
From ! {self(), E - Dt / (Eps * Dx) * DH} | |
end. | |
calc_hz(Const, Hz, Ey) -> | |
Diff_Ey = create_diff_list(Ey), | |
[send_calc_hz_message({H, DE}, Const) || {H, DE} <- lists:zip(Hz, Diff_Ey)]. | |
send_calc_hz_message(EH, Const) -> | |
Pid = spawn(fdtd, calc_hz_main, []), | |
Pid ! {self(), {EH, Const}}, | |
receive | |
{Pid, Msg} -> Msg | |
end. | |
calc_hz_main() -> | |
receive | |
{From, {{H, DE}, {Mu, Dt, Dx}}} -> | |
From ! {self(), H - Dt / (Mu * Dx) * DE} | |
end. | |
ey_boundary({Light, Dt, Dx}, {Ey_at_1, Ey_at_N1}, Hist_Ey) -> | |
C = Light * Dt, | |
[Hist_Ey_0, Hist_Ey_1, Hist_Ey_2, Hist_Ey_3] = Hist_Ey, | |
Ey_0 = Hist_Ey_1 + (C - Dx) / (C + Dx) * (Ey_at_1 - Hist_Ey_0), | |
Ey_N = Hist_Ey_2 + (C - Dx) / (C + Dx) * (Ey_at_N1 - Hist_Ey_3), | |
{Ey_0, Ey_N}. | |
calc(Const, EH_Field, Step) -> | |
{Physics, Area, Time} = Const, | |
{Ey, Hz, Hist_Ey} = EH_Field, | |
Stimulus = calc_stimulus(Step, Time#time.dt), | |
Hz1 = set_soft_source(500, Hz, Stimulus * Time#time.dt), | |
Ey1 = calc_ey({Physics#physics.permittivity, Time#time.dt, Area#area.dx}, Ey, Hz1), | |
{Ey_at_0, Ey_at_N} = ey_boundary({Physics#physics.light, Time#time.dt, Area#area.dx}, | |
{lists:nth(1, Ey1), lists:last(Ey1)}, Hist_Ey), | |
Ey2 = [Ey_at_0] ++ Ey1 ++ [Ey_at_N], | |
Hist_Ey1 = [lists:nth(1, Ey2), lists:nth(2, Ey2), | |
lists:nth(length(Ey2) - 1, Ey2), lists:nth(length(Ey2), Ey2)], | |
Hz2 = calc_hz({Physics#physics.permeability, Time#time.dt, Area#area.dx}, Hz1, Ey2), | |
{Ey2, Hz2, Hist_Ey1, Stimulus}. | |
print_point(Fd, Time, Stim, V1, V2) -> | |
io:fwrite(Fd, "~g, ~g, ~g, ~g~n", [Time, Stim, V1, V2]). | |
print_line(Fd, Step, Hz) -> | |
io:fwrite(Fd, "# Step = ~B~n", [Step]), | |
[io:fwrite(Fd, "~B, ~g~n", [I, Data]) || | |
{I, Data} <- lists:zip(lists:seq(0, length(Hz)-1), Hz)], | |
io:fwrite(Fd, "~n",[]). | |
fdtd_loop(Const, EH_Field, Step, {Fd1, Fd2}) -> | |
{_, _, Time} = Const, | |
Last_Step = Time#time.last_step, | |
_ = if Step rem 500 == 0 -> io:format("Step: ~B~n", [Step]); true -> 0 end, | |
case Step of | |
Last_Step -> file:sync(Fd1), file:close(Fd1), | |
file:sync(Fd2), file:close(Fd2); | |
_ -> {Ey, Hz, Hist_Ey, Stimulus} = calc(Const, EH_Field, Step), | |
print_point(Fd1, Step * Time#time.dt, Stimulus, | |
lists:nth(500, Hz), lists:nth(900, Hz)), | |
_ = if Step rem 100 == 0 -> print_line(Fd2, Step, Hz); true -> 0 end, | |
fdtd_loop(Const, {Ey, Hz, Hist_Ey}, Step+1, {Fd1, Fd2}) | |
end. | |
calc_fdtd(Const) -> | |
{_, Area, _} = Const, | |
Ey = create_list(Area#area.nx + 1), | |
Hz = create_list(Area#area.nx), | |
Hist_Ey = create_list(4), | |
{ok, Fd1} = file:open("fdtd_erlang_point.csv", [write]), | |
{ok, Fd2} = file:open("fdtd_erlang_line.csv", [write]), | |
fdtd_loop(Const, {Ey, Hz, Hist_Ey}, 0, {Fd1, Fd2}). | |
main() -> | |
Const = calc_const(), | |
print_analysis_info(Const), | |
calc_fdtd(Const). | |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment