Finite difference time domain method Puns

I'm a PhD physicist, working in the field of optics and photonics.

Many of our problems and simulations are hugely complex, run on background servers which are expensive to maintain, and which aren't running 100% of the time. Upgrading these servers in the lab happens every few years, but again - at a huge cost.

Id be interested in offloading these tasks onto a decentralized computational engine which is "pay to play" - in that I pay small amounts of ADA tokens for people to solve a set of paralleled equations.

In particular, I'd be interested in solving the finite difference time domain problem as per Maxwell's equations.

There already exists a fairly substantial set of code for these solvers - such as lumerical, etc... I really just want to know if I can produce a dApp which solves the problem instead of doing it on my own machine.

for a better idea of exactly what type of problem I'm trying to solve, read this comment I posted : https://www.reddit.com/r/CardanoDevelopers/comments/lpuytp/limits_of_cardano_decentralized_physics_computing/godyk8x?utm_source=share&utm_medium=web2x&context=3 .

Hello, I am in the midst of writing my dis, and I need help with deriving the finite difference method, for time.

I understand how to derive the finite difference approx (forward, backward, central and symmetric) is formed by manipulating the Taylor series -

however I would like to know how the finite difference method is derived via a similar method to above, for time. (for forward, backward, central and symmetric)

Any explanations, links to videos or textbooks would be greatly appreciated

thank you !

Can I know if there is any website/resource that explains how to solve fin numericals using Gauss-Seidel iteration method to obtain temperature distribution curve? I searched the web and couldn't find any. It'll be useful if the source contains how to solve the problem in MATLAB to obtain the TD curve

sparsity structure of coefficient matrix

I have a coefficient matrix which is obtained through finite difference discretization of a 4d dirac equation in terms of the 12 unknown(4 spins and three colours) with 4 points in each dimension of the space time.Overall this sparse matrix has a size of 3072 x 3072 and its sparsity pattern is shown in the image "sparsity structure of coefficient matrix".I am interested to obtain some specific mass dependent blocks in the matrics(these blocks are obtained when one fixes each of the four spins and write the equation of the three colours for each spin.)

If I dig deep into the structure of the matrics ,then I observe the following pattern repeated every 48 diagonal points as shown in thr imagestructure.I can not deduce any equation ordering from this structure.

My question is how can i find out how the 12 equations for each grid point is ordered and which permutation can help me to achieve the desired structure where colours are arranged per spin.Has anyone worked with a 4d finite difference ?

Hey guys, I'm trying to implement a 2D FDTD algorithm to solve Maxwell's Equations for a gaussian pulse source; but it's not working. The pulse explodes to ridiculous values and doesn't really propagate through the grid... I've gone through everything with a fine-toothed comb, but it's still broken :(

If there are any FDTD afficienardos out there who can scan my code to spot any obvious mistakes I'd be very grateful! Thanks :)

from __future__ import division
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

# SET UNITS
kilohertz = 1e3
megahertz = 1e6
gigahertz = 1e9

# SET PHYSICAL CONSTANTS
c0 = 299792458 
eps0 = 8.85418782e-12
mu0 = 1.25663706e-6

# GRID PARAMETERS
nmax = 1   # Air
NLAM = 50  # Want to resolve the shortest wavelength (highest freq) with 50 points

# SOURCE PARAMETERS
fmax = 5.0 * gigahertz 
tau = 0.5 / fmax
t0 = 5 * tau           

# SET COMPUTATIONAL CONSTANTS
lam0 = c0 / fmax 
dx = lam0 / nmax / NLAM
Nx = 200
dy = dx
Ny = 200
steps = 1000

# CREATE GRID AXIS
xa = np.linspace(0, Nx, Nx, endpoint = False) * dx
ya = np.linspace(0, Ny, Ny, endpoint = False) * dy
xa, ya = np.meshgrid(xa, ya)

# INITIALISE MATERIALS TO FREE SPACE
epsR = np.ones((Nx, Ny))
muR = np.ones((Nx, Ny))

# COMPUTE TIME STEP
dt = min(dx, dy) / (2 * c0 )

# COMPUTE THE SOURCE
time = np.linspace(0, steps, steps, endpoint = False) * dt 
source = np.exp(-((time - t0)/tau)**2) 

# COMPUTE UPDATE COEFFICIENTS
mHx = - c0 * dt / muR
mHy = c0 * dt / muR
mDz = c0 * dt
mEz = 1 / epsR

# INITIALIZE FIELDS TO ZERO
Hx = np.zeros((Nx, Ny))
Hy = np.zeros((Nx, Ny))
Dz = np.zeros((Nx, Ny))
Ez = np.zeros((Nx, Ny))

# MAIN FDTD LOOP
for t in range(steps):

    # COMPUTE DZ
    for i in range(Nx):
        for j in range(Ny):
            if i == 0 and j == 0:
                Dz[0][0] = Dz[0][0] + mDz * ((Hy[0][0] - 0) / dx - (Hx[0][0] - 0) / dy) 
            elif i == 0:
                Dz[0][j] = Dz[0][j] + mDz * ((Hy[0][j] - 0) / dx - (Hx[0][j] - Hx[0][j-1]) / dy) 
            elif j == 0:
                Dz[i][0] = Dz[i][0] + mDz * ((Hy[i][0] - Hy[i-1][0]) / dx - (Hx[i][0] - 00) / dy)
            else:
                Dz[i

I searched a lot and found methods such as Gelerkin and collocation methods and can anyone confirm if these are similar to finite difference analysis ? My condition is they have to be discretization methods.

Apparently I mixed these two up in my FEA class. We were solving a problem for homework where we had to calculate the stresses and deflections in an I-beam, and divide it into 3 sections (the two flanges and the web). The question explicitly asked to use the direct method of solving, and I got the correct answers, but my prof told me that the direct method is the same as the finite element method...

To be clear, he wanted me to solve it this way:

{stiffness matrix}*{deflections} = {reaction forces}

As opposed to what I thought was the "direct method":

deflection = (PL)/(AE)

Any clarification is greatly appreciated

I'm writing a 2D FDTD program for the TMz case. I'd like to keep the field updating equations in the PML and the medium of the same form, that is with no split in the Ez field. Which means I guess I'd like to just alter the conductivities of the PML layers.

If anyone has a minute to step away from the ~stonks~ and help me with this homework question I would be super grateful! The question is as follows:

Write a matlab routine that solves the two-dimensional Poisson's equation, (del)u = f(x; y); in the domain x = [0; Lx] and y = [0; Ly] subject to the following Dirichlet boundary conditions on the boundaries:

u(x, 0) = sin(x); u(x, Lx) = sin(x)*cos(Ly); u(0, y) = 0; u(Lx, y) = sin(Lx)*cos(y)

Assume a uniform grid with Nx and Ny equally spaced points in the x and y directions, respectively. Build the full matrix (i.e. don't use a sparse matrix), and show that you get second-order accuracy using the test solution:

u(x, y) = sin(x)*cos(y)

for Lx = Ly = 8: Use an equal number of points in each direction, Nx = Ny = N: Note that you will need to build the required forcing term f(x, y).

I have a code that is working, but it does not match up with the analytical solution and I have no idea why. I've been working on it all day and I just can't look at it anymore. I appreciate any help, thank you in advance! I have pasted what I have below.

clear all; close all; 

% length of the area to be solved
Lx = 8; Ly = 8;

% number of grid points
Nx = 100; Ny = 100;

% grid points in the x-direction
dx = Lx/(Nx-1);

for i = 1:Nx
    xx(i,1) = (i-1)*dx;
end

% grid points in the y-direction
dy = Ly/(Ny-1);

for i = 1:Ny
    yy(i,1) = (i-1)*dy;
end

% initialize the size of matrix AA and vector cc
AA = zeros(Nx*Ny,Nx*Ny);
cc = zeros(Nx*Ny,1);

% apply BCs on lower/upper boundaries
for j = 1:Nx
    % lower wall (y = 0)
    row = j;
    AA(row,row) = 1.0;
    cc(row) = sin(xx(j))*cos(yy(1));
    % upper wall (y = Ly)
    row = Nx*(Ny-1)+j;
    AA(row,row) = 1.0;
    cc(row) = sin(xx(j))*cos(yy(Ny));
end

% apply BCs on left/right boundaries
for i = 2:Ny-1
    % left wall (x = 0)
    row = i;
    AA(row,row) = 1.0;
    cc(row) = sin(xx(1))*cos(yy(i));
    % right wall (x = Lx)
    row = Nx*(Ny-1)+i;
    AA(row,row) = 1.0;
    cc(row) = sin(xx(Nx))*cos(yy(i));
end

% interior grid points
for i = 2:Ny-1
    for j = 2:Nx-1
        row = (i-1)*Nx+j;
        AA(row,row) = -2.0*(1.0/(dx*dx) + 1.0/(dy*dy));
        AA(row-1,row) = 1.0/(dx*dx) + 1.0/(dy*dy); % west stencil point, x(j-1)
        AA(row+

As per the discussion topic vote, September's monthly topic is "Finite Element Method vs Finite Volume Method vs Finite Difference Method vs Spectral Element Method vs Hybrid Methods".

Previous discussions: https://www.reddit.com/r/CFD/wiki/index

Hello again! I have finally written a code for a system of 3 second order differential equations (theta_r, theta_b, and theta_c) with Dirichlet (static) boundary conditions. (Code #1)

I am trying to impose Neumann (differential) boundary conditions for theta_b(x=1) and theta_c(x=0) (code #2). I'm unsure if I have implemented this incorrectly or if my values for the heat transfer parameters are off, so that it makes the new graphs flat. I have posted both of my codes below and I appreciate any help! THANKS!

Code #1

import numpy as np
from scipy import linalg as lin
from matplotlib import pyplot as plt

# dimensional boundary conditions
T0 = 30.
Tf = 500.
Tr0 = 30.
Trf = 500.
Tb0 = 30.
Tbf = 500.
Tc0 = 500.
Tcf = 30.

# nondimensional boundary conditions
thetar0 = (Tr0 - T0)/(Tf - T0)
thetar1 = (Trf - T0)/(Tf - T0)
thetab0 = (Tb0 - T0)/(Tf - T0)
thetab1 = (Tbf - T0)/(Tf - T0)
thetac0 = (Tc0 - T0)/(Tf - T0)
thetac1 = (Tcf - T0)/(Tf - T0)

# geometric parameters
LL = 1.
r_r = 1.        # radius of the screw rod
r_b = 4.        # radius of the biomass
r_c = 8.        # radius of the cooling jacket
Acr = np.pi*r_r**2          # cross sectional area of the screw rod
Asr = 2*np.pi*r_r*LL        # surface area of the screw rod
Acb = np.pi*(r_b-r_r)**2    # cross sectional area of the biomass
Asb = 2*np.pi*r_r*LL        # surface area of the biomass--UNSURE WHICH RADIUS TO USE, ASK
Acc = np.pi*(r_c-r_r)**2    # cross sectional area of the cooling jacket
Asc = 2*np.pi*r_c*LL        # surface area of the cooling jacket

# nondimensional geometric parameters
Ar = Acr/Asr
Ab = Acb/Asb
Arb = Asr/Acb
Ac = Asc/Acc
A_cb = Asb/Acc

# dimensionless heat transfer parameters
Bi_b = 0.1
Bi_r = 1.
Bi_c = 0.01
Bi_rb = 0.1
Bi_cb = 0.01
Pe_b = 10.
Pe_c = 1.


# finite difference parameters
h = 0.01                # step size
nx = int(LL/h) + 1      # number of space segments
#nn = nx - 2             # number of nodes to solve

# matrix parameters
d = h**2*Bi_r*Ar        
lam = h**2*(Bi_b*Ab+Bi_rb*Arb)      
eta = h**2*Bi_rb*Arb    
phi = h*Pe_b/2
psi = h*Pe_c/2
bet = h**2*(Bi_c*Ac - Bi_cb*A_cb)
sig = h**2*Bi_cb*A_cb
mu = h**2*Bi_b*Ab

XX = np.linspace(0, LL, num=nx, endpoint=True)
#x_x = len(XX)
#print(

how can I write these nodes which Cartesian and cylindrical coord. junction with finite difference method. I solve equation nodes in Cartesian coord. but I cant other nodes. Anybody help?

https://preview.redd.it/k4yjfa3km0641.jpg?width=961&format=pjpg&auto=webp&s=101199450d8f2df15562508d565ff9bcf993d484

I'm trying to make a graph for a 1D heat transfer equation using the finite difference method. Python keeps giving me the error:

ValueError: matmul: Input operand 1 does not have enough dimensions (has 0, gufunc core with signature (n?,k),(k,m?)->(n?,m?) requires 1)

I am trying to multiply a 99x99 matrix with a 99 vector in line 62.

I have copy/pasted my code below. I appreciate any help thanks!

# finite difference method
# 8/17/20

import numpy as np
from matplotlib import pyplot as plt

# dimensional boundary conditions
T0 = 30.
Tf = 500.
Tr0 = 30.
Trf = 500.

# nondimensional boundary conditions
thetar0 = (Tr0 - T0)/(Tf - T0)
thetar1 = (Trf - T0)/(Tf - T0)

# geometric parameters
LL = 1.
r_r = 1.
Acr = np.pi*r_r**2 
Asr = 2*np.pi*r_r*LL 

Ar = Acr/Asr

# dimensional heat transfer parameters
h_r = 25.
Cp_r = 490.
rho_r = 7976.667
k_r = 15.1
alp_r = k_r/(rho_r*Cp_r)
gam_r = h_r/(Cp_r*rho_r)

# nondimensional heat transfer parameters
Bi_r = gam_r*LL**2/alp_r

# finite difference parameters
h = 0.01            # step size
nx = int(LL/h) + 1  # number of space segments
d = h**2*Bi_r*Ar    # matrix parameter
nn = nx - 2         # number of nodes to solve

XX = np.linspace(0, LL, num=nx, endpoint=True)
#x_x = len(XX)
#print('xx = % f \n' %(x_x))

# create the tridiagonal matrix
Ea = (d - 2)*np.ones(nn)    # main diagonal
Eb = np.ones(nn-1)          # secondary diagonals
M = np.diagflat(Eb,k=-1) + np.diagflat(Ea,k=0) + np.diagflat(Eb,k=1)

# create the constants vector
A = np.zeros(nn)
A[0] = thetar0      # first point in vector is first BC
A[-1] = thetar1     # last point in vector is second BC

# create array to store the solution
theta_r = np.zeros(nn)
theta_r[0] = thetar0
theta_r[-1] = thetar1

for s in range(0,nx-1):
    theta_r[s+1] = np.matmul(M[s,:],theta_r[s])

# make temperature dimensional
T_r = theta_r*(Tf-T0)+T0

# create plot
plt.plot(XX, T_r, '-')
plt.xlabel('X (x/L)')
plt.ylabel('T (C)')