idaholab / moose

Multiphysics Object Oriented Simulation Environment
https://www.mooseframework.org
GNU Lesser General Public License v2.1
1.77k stars 1.05k forks source link

Add RZ ray tracing support #22499

Open loganharbour opened 2 years ago

loganharbour commented 2 years ago

Reason

Ray tracing is only supported in XYZ and R-spherical geometries. We should support RZ geometry. Example use case: #22500

Design

Add support for ray tracing in RZ geometry.

@snschune has the math for making this work, I have the design for making it work. We need to hijack the tracing routines for QUAD4 elems in RZ geometry. The math can be described later, but essentially this involves moving the intersection points on a per-segment basis. We do a change of variable from the xyz intersection to an rz intersection. Kernels should work the same.

Impact

Added capability.

loganharbour commented 2 years ago

@snschune care to comment on the maths?

snschune commented 1 year ago

We need to carry the three-dimensional direction to distinguish directions that make different angles with the radius. To do that let's call v = (vx, vy, vz) the direction in 3D xyz space, p = (x,y,z) is a point in 3D space, and rho = (r,eta) is a point in RZ space. The mesh and the elements are given in RZ space, in particular the normals on faces are denoted by N=(Nr, Neta).

Rays are started at some point in rho space, let's call that point rho_0. We need to compute the intersection of this ray with the first face in the element it's currently in. To do that, we first compute the representation of rho in xyz space:

p_0(rho_0) = (r_0, 0, eta_0) = (x_0, y_0, z_0)

Then we compute the intersection of the ray with face k characterized by the normal form:

(1) N_k * rho(p) - B = 0

where:

(2) rho(p) = (sqrt((x_0 + s * vx)^2 + (y_0 + s * vy)^2), z_0 + s * vz)

where s is a Real parameter. Substituting (2) into (1) gives a quadratic equation for s. Solve for s and only accept the forward solution. Do that for all faces in the element and select the smallest s.

Then compute the next point p_1 by

p_1 = p_0 + s * v

and rho_1 by:

rho_1 = (sqrt(x_1^2 + y_1^2), z_1)

Then rinse repeat for each element until ray dies.