[IB method] Implemented new ghost point method #9

ktoepperwien · 2024-09-03T18:46:27Z

Implemented new ghost point method following Ghias et al., JCP 225(1), 2007, with image point interpolation based on R. Franke, Mathematics of Computation 38(157), 1982. Image points currently have the same distance from IB surface as ghost points. In this commit, only Dirichlet boundary conditions supported! Includes test case for non-reacting flow over 2D hill ("Witch of Agnesi"/"Agnesi hill") based on Lundquist et al., Month. Meterol. Rev. 140(12), 2012.

Notes:

common_ops.py scatter/gather is incompatible with lists, therefore commented out dim-check.
To run Witch of Agnesi test case, be sure to initialize the Coriolis force in fire.py as follows:

   _DEFAULT_LAT = 0.5497607357 # rad == 31.5 deg 
   # Coriolis force coeff of 1e-4 rad/s from Lundquist et al., Month. Meterol.
   # Rev. 2010 for validation of IB methods corresponds to a latitude 
   # defined by _REF_LAT 
   _REF_LAT = 0.7555285998 # rad == 43.289 deg 
   my_lat = _REF_LAT 
   logging.warn(f'Using reference latitude of {my_lat} rad ' 
                f'corresonding to {np.rad2deg(my_lat):.3f} deg.') 
   self.coriolis_force_fn = cloud_utils.coriolis_force(my_lat, { 
       'u': self.fire_utils.u_init, 
       'v': self.fire_utils.v_init, 
       'w': 0.0
   }, 2)

Ghias et al., JCP 225(1), 2007, with image point interpolation based on R. Franke, Mathematics of Computation 38(157), 1982. Image points currently have the same distance from IB surface as ghost points. In this commit, only Dirichlet boundary conditions supported! Includes test case for flow over 2D hill ("Witch of Agnesi"/"Agnesi hill") based on Lundquist et al., Month. Meterol. Rev. 140(12), 2012, http://journals.ametsoc.org/doi/10.1175/MWR-D-11-00311.1 Notes: - common_ops scatter/gather is incompatible with lists, therefore commented out dim-check. - To run Witch of Agnesi test case, be sure to initialize the Coriolis force as follows: _DEFAULT_LAT = 0.5497607357 # rad == 31.5 deg # Coriolis force coeff of 1e-4 rad/s from Lundquist et al., Month. Meterol. # Rev. 2010 for validation of IB methods corresponds to a latitude # defined by _REF_LAT _REF_LAT = 0.7555285998 # rad == 43.289 deg my_lat = _REF_LAT logging.warn(f'@@@ Using reference latitude of {my_lat} rad ' f'corresonding to {np.rad2deg(my_lat):.3f} deg.') self.coriolis_force_fn = cloud_utils.coriolis_force(my_lat, { 'u': self.fire_utils.u_init, 'v': self.fire_utils.v_init, 'w': 0.0 }, 2)

john-qingwang · 2024-09-04T00:04:58Z

There are quite a bit of lint issues that are not conforming with the Google programming style. I can go ahead and correct them, but I think it would be a good practice for you to work on it. You may found https://github.com/google/yapf useful for this purpose. Thanks!

added Neumann boundary condition for scalars, and reformatted code.

ktoepperwien · 2024-09-11T00:28:09Z

Note: To preserve second order accuracy with Neumann boundary conditions, we need 2 additional points in the fluid domain, see Eq. 14 in Luo et al., International Journal of Heat and Mass Transfer 92, 2016 (Link, DOI: 10.1016/j.ijheatmasstransfer.2015.09.024). This is not yet implemented -- we are currently using a single image point for Neumann boundary conditions.
For non-isothermal cases, we will also need to update the density at the ghost point by invoking the equation of state (which is only first-order accurate). For second-order accuracy, see Eqs. 17 and 18 in above reference.

john-qingwang · 2024-09-12T16:37:10Z

swirl_lm/utility/common_ops.py

@@ -1574,8 +1574,8 @@ def strip_halos(
  # Handles the case of single 3D tensor.
  if isinstance(f, tf.Tensor):
    nz, nx, ny = f.get_shape().as_list()
-    return f[halos[2]:nz - halos[2], halos[0]:nx - halos[0],
-             halos[1]:ny - halos[1]]
+    return f[halos[2]: - halos[2], halos[0]: - halos[0],


Why do we need to change these lines? The previous implementation is supposed to be more robust than the updated one. Is it because f is a dynamic tensor whose shape is None at this point? Can we try tf.shape and see it works?

john-qingwang · 2024-09-12T16:38:48Z

swirl_lm/utility/common_ops.py

@@ -1668,8 +1668,8 @@ def gather(x: tf.Tensor, indices: tf.Tensor) -> tf.Tensor:
    matching the first dimension of `indices`.  If `indices` is empty, return
    a vector with length 0.
  """
-  if  tf.shape(indices)[0] == 0:
-    return tf.constant([])
+  # if  tf.shape(indices)[0] == 0: # <-- Unsupported!


Let's remove these comments. Thanks!

john-qingwang · 2024-09-12T16:39:04Z

swirl_lm/utility/common_ops.py

@@ -1721,8 +1721,8 @@ def scatter(
    with everywhere else being 0. If `indices` is empty, a 3D tensor with all
    zeros will be returned.
  """
-  if tf.shape(indices)[0] == 0:
-    return tf.zeros(shape, dtype=dtype)
+  # if tf.shape(indices)[0] == 0: # <-- Unsupported!


john-qingwang · 2024-09-12T16:45:12Z

swirl_lm/example/shared/wildfire_utils.py

@@ -563,7 +563,8 @@ def __init__(

    if (self.config.combustion is None or
        not self.config.combustion.HasField('wood')):
-      raise ValueError('Wood model is not defined as a combustion model.')
+      # raise ValueError('Wood model is not defined as a combustion model.')


As this is a library for fire utilities, we have to make sure that a combustion model is defined. We may make the instance of fire_utils that's created by the caller None as opposed to having it passing through silently here.

john-qingwang · 2024-09-12T16:46:57Z

swirl_lm/example/shared/wildfire_utils.py

-    self.drag_force_fn = self.vegetation_drag_update_fn(
-        _C_D.value, self.config.combustion.wood.a_v
-    )
+    if not self.config.combustion is None:


Let's keep the conditional check for combustion model at the top level, and assume that it's never a None inside the fire utils.

john-qingwang · 2024-09-13T20:34:55Z