flaport / fdtd

A 3D electromagnetic FDTD simulator written in Python with optional GPU support
https://fdtd.readthedocs.io
MIT License
455 stars 116 forks source link

Objects of arbitrary shape #3

Closed BrandonHanx closed 2 years ago

BrandonHanx commented 4 years ago

Thanks for your code! And I wonder how to define a structure like circle? More generally, how can I define an arbitrary shape? Thanks.

flaport commented 4 years ago

You're scratching the surface of something that has not explicitly been implemented (yet).

On top of that, while preparing an example to help you further I found a bug (typo) which would have prevented you to implement it in the first place. The bug is now removed, but bear in mind that the following example will only work with the development version of this library (pip install git+https://github.com/flaport/fdtd)

The goal is definitely to have an explicit implementation at some point, but for now you can use the following workaround (the full python file is attached in the zip):

# the following only works for a numpy backend. Change the numpy arrays
# to pytorch tensors to make it work for the torch backend.

refractive_index = 1.7
x = y = np.linspace(-1,1,100)
X, Y = np.meshgrid(x, y)
circle_mask = X**2 + Y**2 < 1

# 3D permittivity for each voxel in the object
permittivity = np.ones((100,100,1))
permittivity += circle_mask[:,:,None]*(refractive_index**2 - 1)
grid[170:270, 100:200, 0] = fdtd.Object(permittivity=permittivity, name="object")

Obviously, this works for any kind of object (2D, 3D, ...) as long as you can define a mask for it.

However, as you can see, the visualize function does not know how to handle a tensor with locally changing permittivity and shows for both the cases (higher index in the full square vs higher index in the enclosed circle) a pink rectangle. However, as we can see from the fields, the behavior is different.

circle: (notice the focusing) circle_mask

square: rectangle_mask

The full python file if you want to play around with this: circular_object.zip

BrandonHanx commented 4 years ago

Thanks for your answer. Now I understand it.

Kiranalu commented 4 years ago

Thanks for your answer. Now I understand it.

Man,do you try to make circle look like circle?

sofroniewn commented 3 years ago

Hello! First off I'd like to say that I'm loving this package! I've only discovered it today, but I'm already having loads of fun with it and I think it's going to work really well for my needs, so thanks!

I have a question about this topic though. If I run the example from here https://github.com/flaport/fdtd/issues/3#issuecomment-523802210

grid = fdtd.Grid(shape=(300, 300, 1))

refractive_index = 1.7
x = y = np.linspace(-1,1,100)
X, Y = np.meshgrid(x, y)
circle_mask = X**2 + Y**2 < 1

# 3D permittivity for each voxel in the object
permittivity = np.ones((100,100,1))
permittivity += circle_mask[:,:,None]*(refractive_index**2 - 1)
grid[170:270, 100:200, 0] = fdtd.Object(permittivity=permittivity, name="object")

plt.imshow(np.linalg.norm(grid.inverse_permittivity[:, :, 0, :], axis=2))

and then plot the inverse_permittivity I just get a "hole" where my object is

I was naively (and I mean very naively as I am new to these sorts of simulations) expecting to see the circle in side the grid.inverse_permittivity but it appears that is not the case, instead there is just 0 there. I guess that is because you don't actually set the values back into the grid. I'm curious why, one guess I had is it allows you to possibly have objects of different "resolutions", but maybe I'm missing something.

I'm planning to write a function right now which would look through all the objects, extract there permittivity values and "stitch" them into a copy of the grid, so that I could do something like imshow and see the full environment that I had created, but before doing that I wanted to ask if that seems like a reasonable thing to do.

If you thought it was useful to others I'd happily contribute it to this repo as a PR

sofroniewn commented 3 years ago

Just a quick followup - I'm able to achieve what I want with something like

stiched_inverse_permittivity = grid.inverse_permittivity.copy()

for obj in grid.objects:
    stiched_inverse_permittivity[obj.x, obj.y, obj.z] = obj.inverse_permittivity

plt.imshow(np.linalg.norm(stiched_inverse_permittivity[:, :, 0, :], axis=2))

But I'm curious if that's the recommended way to achieve what I want, if that will always work - for example what happens if there are multiple objects in the same location in the grid? I'm not sure if that's possible right now, or if things are additive etc, but one would probably want to consider situations like that before adding a general method.

flaport commented 3 years ago

Hi @sofroniewn ,

Thanks for liking the library 🙂

Indeed, placing an object in the grid sets the inverse_permittivity to zero within the location of the object. In those locations, the inverse_permittivity of the Object is used (I'm not a 100% happy with this solution, see for example #4 and #6).

The reason I'm doing this however is to indeed allow for different update equations within the object as opposed to update equations within the grid. For example, AnisotropicObject has a different inverse_permittivity in x, y and z direction and AbsorbingObject includes absorption. An object that uses a smaller internal grid size would be a great idea to implement as well!

At the moment, stitching the inverse_permittivity together like you did is currently the only way to accurately visualize the index contrast in the grid.