Use 'quantities' library to handle conversions between simulation units and actual units

[aizvorski]

I'm wondering what the units used for each quantity are, both for parameters which are passed in and for values stored internally or produced by the calculation.

A lot of things seem to be unitless, but other things seem to have assumed units. A quick list:

• grid.E - the docs say "E := √ε0E" - does this mean that grid.E divided by sqrt(numeric value of epsilon0 expressed in farads/meter) (= grid.E/2.98e-6) would produce the E field in volts/meter? is this basically because sqrt(epsilon0)sqrt(mu0) = speed of light?
• grid.H?
• Grid(grid_spacing) in meters?
• grid.time_step in seconds?
• AbsorbingObject(conductivity) in siemens/meter? (really not sure about this one but I need it to match specific materials)
• PointSource(power) in watts?

It would be great to track all quantities with units and add this to the documentation. I can try to write that up as long as someone can help answer questions.

[flaport]

Hi @aizvorski ,

• Indeed, choosing E_sim= √ε0*E and H_sim=õ0*H basically results in nicely symmetric update equations.
• grid_spacing is indeed in meters
• time_step is indeed in seconds (derived from grid_spacing and courant_number)
• The conductivity σ should indeed be given by S/m = (A^2 s^3) / (kg m^3), however in the chosen simulation units we have σ_sim = inv(ε) * σ / ε0, i.e. to obtain the conductivity in simulation units, divide by the permittivity of the medium (I'm open to change this convention).
• The power of a source should ideally be watts or alternatively be a varying current density J (with units A/m^2). Or in stead of power we should just define a electric field amplitude that's being added to the grid, but I think the current implementation is neither of that... I honestly forgot my reasoning in implementing it this way... This probably needs to be standardized and cleaned up.

[aizvorski]

Hi @flaport - Thanks, that's really helpful!

I'm still slightly confused by "simulation units" :) There may be two different interpretations (which may be equivalent, maybe)

(a) the equations are changed to curl(H_sim) = ε*(1/c)*dE_sim/dt and curl(E_sim) = -µ*(1/c)*dH_sim/dt where H_sim and E_sim are new physical quantities which both have the same dimensions of [mass]^1/2 [length]^-1/2 [time]^-1, and c is a velocity as usual?

(b) the equations are not changed, but E is expressed in units of 3.36e+5 volts/meter, and H is expressed in units of 8.92e+2 amperes/meter, where these are just the original quantity multiplied by a dimensionless scaling factor, ie when you say sqrt(epsilon0) you don't mean the sqrt of the physical quantity with units but rather sqrt(epsilon0 / (1 farad/meter)) which is dimensionless. All the scaled values still have the same dimensions, and the equations are the same, but are in units chosen so that the numerical values of epsilon0 and mu0 (both expressed in those units) happen to be equal to each other? (*and I think the only unit that is being scaled is current, everything else can stay in SI units)

Which of these two is more correct? Are they equivalent?

Also is the main advantage of this to make sure both equations can use the same Courant number?

[aizvorski]

@flaport "instead of power we should just define a electric field amplitude" - I like this, it's intuitive. It seems that the point/line sources convert power to amplitude and then just use amplitude internally, might as well specify amplitude as input.

Looking at the code, amplitude is distributed proportionally to inv(epsilon) along lines? (and maybe planes?) That may produce unexpected results if a source is inside an object :)

I'd like to set up a simulation which has a point dipole (short line) with a given amplitude in V/m and measure power density in the far field in W/m^2, and make sure it matches the analytical values and all the units and conversions are as expected.

If I wanted to try that with the current definition of line sources, what are the power units now - are they essentially E_sim units*H_sim units*meter^3/second?

[flaport]

Hi @aizvorski,

(a) is correct. Just be know that in simulation units curl(H) has the same unit as E and H. i.e. the curl does not divide by the grid spacing, in fact we have that ∇×H = curl(H)/du (see this, where I note that the factor 1/du was left out). Adding this missing factor to your equation in (a) gives the following update equations:

E  += (c*dt/du)*inv(ε)*curl_H = sc*inv(ε)*curl_H
H  -= (c*dt/du)*inv(µ)*curl_E = sc*inv(μ)*curl_E

This means that the unit of E and H is [mass]^(+1/2) * [length]^(+1/2) * [time]^(-1) (notice the positive sign in the exponent for the length dimension).

The advantage of this convention is that E and H will be of comparable magnitude, which should make the FDTD algorithm numerically more stable.

Working back from the source value being added to the grid, we have that

amplitude ~ kg^(1/2) * m^(1/2) * s^(-1)

Which means that the unit for power must be:

power ~ amplitude^2 ~ kg m / s^2 = N

This of course makes no sense. I am going to change this.

The idea behind scaling the amplitude of the source with n=sqrt(ε^{-1}) is that I wanted people to specify the free-space amplitude/power for the source, but it might indeed might be counterintuitive (and maybe even wrong)... I am going to change this too.

[flaport]

Ok, the power keyword issue was solved in 1c3b98dbc5a43bc5231db0d883760a81d4a1e3d9 (fdtd version 0.2.0).

From now on, you need to specify the field amplitude in simulation units.

[pingpongballz]

Hi,

From what I'm understanding by reading this, I'm assuming one can convert E from simulation units to V/m by simply dividing the simulation units by √ε0 and H_sim to A/m by dividing by √μ0?

Thanks!

[aizvorski]

@flaport Thanks, very cool! I'm going to try specifying a point dipole with known amplitude and see if the detected power vs distance in the far field matches the expected value. Would you like to have an example like that in the standard examples?

Ok, the power keyword issue was solved in 1c3b98d (fdtd version 0.2.0).

From now on, you need to specify the field amplitude in simulation units.

[aizvorski]

@flaport Thanks for explaining, good to know that these are new quantities with new units. I think I'm getting close to actually understanding this :)

This means that the unit of E and H is [mass]^(+1/2) * [length]^(+1/2) * [time]^(-1) (notice the positive sign in the exponent for the length dimension).

I see the formula but I admit I don't understand why that means the units are as described. (c*dt/du) is dimensionless (=velocity*inverse velocity) so that doesn't really narrow down what the H_sim and E_sim units are, as long as they are the same.

Going from the E_sim= √ε0*E and H_sim=õ0*H formulas: E is in volts/meter and ε0 in farads/meter, E_sim is in [mass]^1/2 [length]^-1/2 [time]^(-1); and H is in amperes/meter and µ0 in henry/meter gives H_sim in [mass]^1/2 [length]^-1/2 [time]^(-1) as well (and they do match)

So where does the [length]^+1/2 come from? I think you mean it comes from multiplying by du, but that's not actually multiplying as much as redefining curl. ∇×H_sim would have dimension of H_sim/length; curl(H_sim) (the function in the code which doesn't divide by du) has the same dimension as H_sim. Could you please explain this part?

[aizvorski]

From what I'm understanding by reading this, I'm assuming one can convert E from simulation units to V/m by simply dividing the simulation units by √ε0 and H_sim to A/m by dividing by √μ0?

@pingpongballz I think that's correct, as long as ε0 is in farads/meter and μ0 is in henry/meter. The simulation units for both E and H should be kilogram^(+1/2)*meter^(-1/2)*second^(-1). Taking kilogram^(+1/2)*meter^(-1/2)*second^(-1) / sqrt(farad/meter) gives volt/meter, so that seems to check out.

However @flaport said it should be meter^(+1/2), so I may be totally wrong about that (and the du value may be needed as well). I asked for clarification, so please stay tuned :)

[pingpongballz]

@aizvorski thanks so much! My initial doubt was also in the metre^0.5 issue. But ill just bite the bullet and go ahead and multiply my poynting vectors by 3e8 m/s.

Again, thanks for your help! 😃

[flaport]

Hi @aizvorski , @pingpongballz ,

Indeed, I seem to have made a mistake when I was quickly double checking the units...

Esim = ε0^(1/2) E
~ F^(1/2) m^(-1/2) V m^(-1)
~ C s kg^(-1/2) m^(-3/2) m kg C^(-1) s^(-2)
~ kg^(1/2)  m^(-1/2) s^(-1)

My apologies for the confusion.

Closing this issue as I think this has been resolved. But feel free to continue the discussion if you want any of this changed.

[aizvorski]

@flaport Thank you for the excellent explanation - I hope this is useful to anyone who has the same questions in the future. Indeed mostly resolved.

How do you feel about using a package like quantities to allow passing in values (eg amplitude, conductivity) and retrieving values (eg field strength from detectors) with proper units? It could be an optional dependency - if not installed (of if passing just floats or numpy arrays) assume simulation units, but if passing in a value with units then do the conversion to/from simulation units as needed.

[flaport]

That's indeed a good idea! I will look into this. I'm reopening the issue to track progress on this issue.