1
$\begingroup$

I am trying to accurately calculate $\vec{j} \cdot \vec{E}$ for an electron energy equation on a finite-volume mesh. $\vec{j}$ is the electron current density, and $\vec{E}$ is the electric field. On the vertex centered mesh, I know $\phi$ at the center of each cell (where $\vec{E} = \nabla \phi$) and I also have a function that calculates $\vec{j} \cdot \hat{n}$ at each face of each cell.

Using the $\phi$ values, I can easily calculate $\vec{E} \cdot \hat{n}$. The problem is that multiplying these two normal components leaves out the tangential components. Taking the maximum over the faces leads to 10-20% error and volume averaging them is completely incorrect (by several factors). The current solution I am using, is that I also happen to know the direction of $\vec{E}$ although I don't fully trust it. From the physics, it is a fair assumption in most cases that $\vec{j}$ is close in direction to $\vec{E}$, so I assume that they are in the same direction, and divide $\left(\vec{E} \cdot \hat{n}\right) \left(\vec{j} \cdot \hat{n}\right)$ by $\cos^2{\theta}$ which gets me close although I have to ignore faces where the normal is within 30 degrees of perpendicular to $\vec{E}$.

I have also found a solution that works but utilizes a large number of trig functions, which is likely too slow to be done at every cell in every iteration of a unsteady FV solve. The solution works as:

Given $\vec{v}\cdot \hat{n}_1$, $\vec{v} \cdot \hat{n}_2$, $\hat{n}_1$, and $\hat{n}_2$

v = 2*rand(1,2) - 1;    % Original vector we want to reproduce

n_1 = 2*rand(1,2) - 1;  % unit normal vectors
n_1 = n_1/norm(n_1);
n_2 = 2*rand(1,2) - 1;
n_2 = n_2/norm(n_2);

vn_1 = v.'*n_1;         % dot products
vn_2 = v.'*n_2;

thetaN1 = atan2(n_1(2), n_1(1));      % Calculate angle of first normal to x-axis
phi = thetaN1 - atan2(n_2(2),n_2(1)); % Difference in angle between normals

theta1 = atan((vn_1*cos(phi) - vn_2)/(vn_1*sin(phi)));  %Use both dot products to back out angle between first normal and v
v_mag = vn_1/cos(theta1);  % Back out magnitude of v

v_out(1) = v_mag*cos(thetaN1 + theta1);  % Calculate components of v
v_out(2) = v_mag*sin(thetaN1 + theta1);

The above reproduces v accurately down to almost double precision. Even if I introduce small amounts of noise into the original v of about 1% for each of the dot products (representing that the field $\vec{v}$ is changing in space from the center of one face to the center of another), it is still accurate to a few decimal places most cases. As I said, the cost of this is higher than I would like. I would need to do the above calculation a few million times per time step.

I feel like I have somewhat reinvented the wheel and was wondering if there was a more efficient and common way to calculate the dot product of two fields when only the normal components at faces are known.

$\endgroup$
1
$\begingroup$

I would have expected a geometric discretization to place $E$ and $j$ along edges of the mesh (more precisely, $j$ would be fluxes through facets dual to the mesh edges), in which case $\phi$ is represented by scalar values at the vertices of the mesh. In this way, $E$ and $J$ always "live" at the same mesh elements and are collinear, allowing the dot product to be trivially evaluated. In your form, it seems like the best thing to do is just to refine the mesh.

Your snippet of code can be optimized substantially. I common trick in graphics is to notice that $\sin \theta$ is simply the magnitude of the cross product of two unit vectors that span the angle, and $\cos \theta$ is their dot product. This allows you to completely avoid ever calculating trig functions, and reduces everything to vector algebra. Also, you frequently don't need to normalize vectors since if you end up taking ratios the magnitudes cancel out.

$\endgroup$
  • $\begingroup$ Yes, J and E are co-located at the center of each face, and their normal components are obviously colinear, but (J.n)*(E.n) does not equal J.E. On a hexagonal mesh (the dual of an equilateral triangle mesh), no matter the refinement level, if J and E both happen to be directly between the normals of adjacent faces, the error would be equal to cos(30)^2. $\endgroup$ – Godric Seer Jun 3 '14 at 15:43
  • $\begingroup$ This assumes that the degrees of freedom you have for $j$ and $E$ are point-wise samples of an underlying continuous field. Alternatively, if you view your degrees of freedom as coefficients of a basis set localized to each facet, then a simple dot product of normal components does make sense since "only" those normal components are represented by the basis. (The basis captures only the normal components, or rather, it's a basis set for a function space that is everywhere normal to the facets). $\endgroup$ – Victor Liu Jun 3 '14 at 17:38
  • $\begingroup$ I understand what you are saying, but I would argue that (at least for my application) I have to interpret it as sampling of a continuous field. A 10% error in heating in the simulation from what the real continuous fields describe leads to certain chemical process rates to be off by as much as 4-6 orders of magnitude. I see that the basis only describing part of the space may be a larger problem than I realized. $\endgroup$ – Godric Seer Jun 3 '14 at 17:45

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.