Correct the weights for dipole sources in cylindrical coordinates #2476
Conversation
| amps_array[idx_vol] = | ||
| IVEC_LOOP_WEIGHT(s0, s1, e0, e1, dV0 + dV1 * loop_i2) * amp * data->A(rel_loc); | ||
| else | ||
| amps_array[idx_vol] = IVEC_LOOP_WEIGHT(s0, s1, e0, e1, 1) * amp * data->A(rel_loc); |
There was a problem hiding this comment.
This can block can be more cleanly expressed as:
double interp_weight = fc->gv.dim == Dcyl ? dV0 + dV1 * loop_i2 : 1.0;
amps_array[idx_vol] = IVEC_LOOP_WEIGHT(s0, s1, e0, e1, interp_weight) * amp * data->A(rel_loc);There was a problem hiding this comment.
Sure, but let's make sure things work before we golf the code.
| return sum_to_all(total_sources); | ||
| } | ||
|
|
||
| std::complex<double> fields_chunk::sum_sources(field_type ft) { |
There was a problem hiding this comment.
Is this function just used for debugging and was forgotten to be removed from the commit?
There was a problem hiding this comment.
No I think we should keep it, as mentioned in the code comment above and the blurb above:
Finally, I added some convenience functions that sum up all the source amplitudes (i.e. perform a sort of "integration") to check if what we're doing is right.
| if (where.in_direction(d) == 0.0 && !nosize_direction(d)) { // delta-fun | ||
| if (d == R && where.center().in_direction(d) > 0) | ||
| // in cylindrical coordinates, we need to scale by $r$ for a dipole. | ||
| data.amp *= 1.0 / where.center().in_direction(d); |
There was a problem hiding this comment.
Why was the scaling by
There was a problem hiding this comment.
See the code comment below:
// we handle the resolution scaling at the chunk level for cylindrical coordinates.
|
It might be worthwhile to think this through more carefully from a finite-element perspective (thinking of finite differences as equivalent to finite elements with first-order elements and a uniform grid, which is typically the case). Then the normalization is an integral of the currents "tent functions" times r. This also gives you a better prescription for what to do with an r=0 gridpoint, which in FEM would be "half a tent" function. |
|
It seems this PR does not produce the expected result for the test in #2108. The radiated flux normalized by the dipole position
The results do not change with increasing resolution. Results below for resolution of 100. This is similar to the master branch as demonstrated in #2470 (comment).
|
As discussed, it's challenging to get the weights just right for a dipole source in cylindrical coordinates. This PR attempts to address that.
First, I modified
sources.cppto include therfactor from thedV0anddV1weight factors, just like what's implemented when processing DFT fields. Previously, this factor was completely omitted.Second, for dipole sources, I include the
1/rterm that is necessary for a proper delta function in cylindrical coordinates.Finally, I added some convenience functions that sum up all the source amplitudes (i.e. perform a sort of "integration") to check if what we're doing is right.
Using these three changes, we should be able to confirm that for a dipole source, it's integral should always be 1 (modulo the resolution scale factor) even for dipoles that lie between grid points. Previously, that wasn't the case. Now, it is. I ran a simple experiment that sweeps a dipole between two grid points and sum up the corresponding weights. Indeed, it sums to 1 for all values of$r$ (there's some ambiguity here at $r=0$ ).
Currently, some of the tests are failing. This could be due to a few different factors. Before we tackle that, I want to ensure we get the weight formulation down first.