flaport / fdtd

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

Updating fields #11

Closed tiagovla closed 3 years ago

tiagovla commented 3 years ago
# grid.py
def update_E(self):
        """ update the electric field by using the curl of the magnetic field """

        curl = curl_H(self.H)
        self.E += self.courant_number * self.inverse_permittivity * curl

        # update objects
        for obj in self.objects:
            obj.update_E(curl)

Here you update E for the whole grid. Then you update E inside the objects, changing E in-place, locally.

#objects.py
def update_E(self, curl_H):
        """ custom update equations for inside the anisotropic object
        Args:
            curl_H: the curl of magnetic field in the grid.
        """
        loc = (self.x, self.y, self.z)
        self.grid.E[loc] *= (1 - self.absorption_factor) / (1 + self.absorption_factor)
        self.grid.E[loc] += (
            self.grid.courant_number
            * self.inverse_permittivity
            * curl_H[loc]
            / (1 + self.absorption_factor)
        )

Aren't you doing it twice?

# first time
    E += c*dt*inv(ε)*curl(H)
# second time
    f = 0.5*dt*σ
    E *= inv(1 + f) * (1 - f)
    E += inv(1 + f)*sc*inv(ε)*curl_H
flaport commented 3 years ago

No, when adding an object to the grid, the inverse_permittivity at the location of the object is set to zero:

https://github.com/flaport/fdtd/blob/73afcb8403aa64193a042680aa0e5dd13d62dae1/fdtd/objects.py#L86

This means that effectively, only the second update equation is used at the location of the object.

estebvac commented 2 years ago

Hi flaport, In that sense, wasn't it better if the method _register_grid in object.py set the permittivity of the object directly in the grid. Repeating calculations can be critical when the simulation size is bigger.

flaport commented 2 years ago

That's a fair point and in fact nothing is stopping you to do exactly that... Just do an in-place index replacement within inverse_epsilon...

The reason that I chose the Object route is that back in the day I was working on simulating photo-refractive crystals in which I was also modeling the electro-optic effect. Moreover those objects required non-diagonal epsilon elements which is a much more expensive update. Going the object-route allowed me to just have those expensive update equations in a localized component and still use the quick update equations everywhere else in the grid.

That said, for normal use cases a replacement in the grid is a much better approach and I agree we should probably make that the default.