Terrain-Aware Immersed Forcing Method #1115

hgopalan · 2024-07-02T18:05:45Z

Summary

This draft PR introduces an immersed forcing method for representing terrain within AMR-Wind. The terrain is generated outside amr-wind and this method does not rely on the existing immersed boundary method. The terrain runs are performed in steps:

Run precursor
Generate the terrain using the included utility from the mmctools project and SRTM data
Generate roughness using NLCD database (under development and not fully integrated)
Run the terrain simulation with the terrain-related forcing terms included in the amp-wind input file and a terrain file created using the utilities.

The method can also be used to run cases without terrain such as an array of bluff bodies. The user has to create a XYZ file with the relevant information to do perform these simulations. The user can also use the stl2XYZ.py in tools/terrain folder to convert a STL file to XYZ file for AMR-Wind. It is recommended to have a smooth boundary near the inflow and outflow to reduce the number of iterations of the nodal and Mac projections.

The terrain modules have been designed to run with the Kosovic turbulence model. To use other turbulence models, the terrainBlanking information should be used to set the turbulent viscosity to zero in the solid region.

Pull request type

This is a draft PR and will be updated over the next weeks.

Please check the type of change introduced:

Checklist

The following is included:

new unit-test(s)
new regression test(s)
documentation for new capability

This PR was tested by running:

the unit tests
- on GPU
- on CPU (Kestrel and MacOS)
the regression tests
- on GPU
- on CPU (Kestrel and MacOS)

Additional background

Issue Number:

… changes.

moprak-nrel

Here's some quick comments from a first pass

moprak-nrel · 2024-07-08T14:53:48Z

amr-wind/equation_systems/icns/source_terms/DragForcing.cpp

+    const amrex::Real gpu_startY = (m_spongeDistanceY > 0)
+                                       ? prob_hi[1] - m_spongeDistanceY
+                                       : prob_lo[1] - m_spongeDistanceY;
+    const amrex::Real gpu_spongeX = (m_spongeDistanceX > 0) ? 1 : 0;


Do these need to be amrex::Real? Looks like the value is just 0 or 1.

Does not have to be real. Can change it to int.

moprak-nrel · 2024-07-08T15:01:01Z

amr-wind/equation_systems/icns/source_terms/DragForcing.cpp

+        amrex::Real gpu_spongeVelZ = 0.0;
+        amrex::Real residual = 1000;
+        amrex::Real height_error = 0.0;
+        for (unsigned ii = 0; ii < vsize; ++ii) {


There might be a better way to find the first closest index to the reference height x3 using std::lower_bound and std::upper_bound or some such instead of this loop?

amr-wind/physics/TerrainDrag.cpp

moprak-nrel · 2024-07-08T15:18:42Z

amr-wind/equation_systems/icns/source_terms/DragForcing.cpp

+        const amrex::Real uy1 = vel(i, j, k, 1);
+        const amrex::Real uz1 = vel(i, j, k, 2);
+        const amrex::Real m = std::sqrt(ux1 * ux1 + uy1 * uy1 + uz1 * uz1);
+        if (drag(i, j, k) == 1) {


Similar to Michael's comments, is drag an integer/bool field? Also, perhaps call this something more descriptive than drag to not confuse this with m_drag?

amr-wind/physics/TerrainDrag.cpp

moprak-nrel · 2024-07-08T15:33:30Z

amr-wind/turbulence/LES/Kosovic.cpp

@@ -18,7 +18,7 @@ Kosovic<Transport>::Kosovic(CFDSim& sim)
    , m_vel(sim.repo().get_field("velocity"))
    , m_rho(sim.repo().get_field("density"))
    , m_Nij(sim.repo().declare_field("Nij", 9, 1, 1))
-    , m_divNij(sim.repo().declare_field("divNij", 3))
+    , m_divNij(sim.repo().declare_field("divNij", 3, 1, 1))


I thought @marchdf changed this line as a bugfix in #1030, is there a reason to change this back?

I updated to latest repo updated and included the changes. Not sure why it changed back.

marchdf

Made a first pass at the non code stuff. My main comment has to do with how we want to present this work.

marchdf · 2024-07-08T15:27:02Z

test/CMakeLists.txt

@@ -242,6 +242,9 @@ add_test_re(burggraf_flow)
 add_test_re(abl_godunov_rayleigh_damping)
 add_test_re(rankine)
 add_test_re(rankine-sym)
+add_test_re(terrainbox)
+add_test_re(nrel_precursor)


If you want to test the whole chain you can add test dependencies. The input-output runs do this.

marchdf · 2024-07-08T15:27:42Z

test/test_files/nrel_terrain/README

@@ -0,0 +1,6 @@
+This is the sample terrain case. You need to run "python3.XX backendinterface.py NREL.yaml" in the tools/terrain folder.


Let's remove this readme (kinda strange to have that in a reg test folder) and move it to proper documentation.

Also, remove/move the png from here.

marchdf · 2024-07-08T15:27:56Z

test/test_files/terrainbox/README

@@ -0,0 +1,6 @@
+This is a simple test case to verify the terrain setup. A box is created and a freestream velocity


Same comment about the readme here.

marchdf · 2024-07-08T15:28:59Z

test/test_files/terrainbox/computeError.py

@@ -0,0 +1,13 @@
+from netCDF4 import Dataset


Let's remove this file or follow the MMS test template (which has a python call). Still my preference would be to remove/move to documentation.

marchdf · 2024-07-08T15:29:40Z

test/test_files/terrainbox/createTerrain.py

@@ -0,0 +1,15 @@
+import numpy as np


Same thing here. I guess these regtests are really "tutorials". I like the regtest and think the tutorial part should be moved to documenation.

marchdf · 2024-07-08T15:31:37Z

tools/terrain/SRTM_to_STL_example.py

@@ -0,0 +1,414 @@
+def SRTM_Converter(outputDir,refLat,refLon,refHeight,left,right,bottom,top):


this looks copy-pasted from jupyter. Can we make this a proper script?

marchdf · 2024-07-08T15:32:07Z

tools/terrain/SRTM_to_STL_example.py

+    # In[26]:
+
+
+    if getWRFdata is True:


shouldn't need to test like this

marchdf · 2024-07-08T15:32:57Z

unit_tests/wind_energy/test_abl_terrain.cpp

+    pp.addarr("physics", physics);
+    amr_wind::terraindrag::TerrainDrag terrain_drag(sim());
+    terrain_drag.post_init_actions();
+    int value = 100;


marchdf · 2024-07-08T16:03:10Z

amr-wind/equation_systems/icns/source_terms/DragForcing.H

+    amrex::Gpu::DeviceVector<amrex::Real> gpu_vel_ht;
+    amrex::Gpu::DeviceVector<amrex::Real> gpu_vel_vals;
+    amrex::Real m_drag{100.0};
+    amrex::Real m_spongeStrength{1.0};


let's maintain the consistent naming case. These should be m_sponge_strength etc.

marchdf · 2024-07-08T16:03:55Z

amr-wind/equation_systems/icns/source_terms/DragForcing.cpp

+{
+    amrex::ParmParse pp("DragForcing");
+    pp.query("dragCoefficient", m_drag);
+    pp.query("spongeStrength", m_spongeStrength);


same remark about case here. I think we typically use snake case in the input files

marchdf · 2024-07-08T16:04:24Z

amr-wind/equation_systems/icns/source_terms/DragForcing.cpp

+    const auto& phy_mgr = m_sim.physics_manager();
+    if (phy_mgr.contains("ABL")) {
+        const auto& abl = m_sim.physics_manager().get<amr_wind::ABL>();
+        const VelPlaneAveraging& fa_velocity =


const auto&

marchdf · 2024-07-08T16:05:00Z

amr-wind/equation_systems/icns/source_terms/DragForcing.H

+    const CFDSim& m_sim;
+    const amrex::AmrCore& m_mesh;
+    const Field& m_velocity;
+    amrex::Gpu::DeviceVector<amrex::Real> gpu_vel_ht;


rename gpu_ to device_

marchdf · 2024-07-08T16:05:29Z

amr-wind/equation_systems/icns/source_terms/DragForcing.cpp

+{
+    const auto& vel =
+        m_velocity.state(field_impl::dof_state(fstate))(lev).const_array(mfi);
+    bool is_terrain = this->m_sim.repo().field_exists("terrainBlank");


marchdf · 2024-07-08T16:09:36Z

amr-wind/physics/TerrainDrag.H

+    // Terrain Drag Force Term
+    Field& m_terrainDrag;
+    // Reading the Terrain Coordinates from  file
+    amrex::Vector<amrex::Real> m_xterrain, m_yterrain, m_zterrain;


just to be consistent with other things, can you declare on 1 line? Even better an amrex::Array of size AMREX_SPACEDIM. Same comment for line 56

marchdf · 2024-07-08T16:10:10Z

amr-wind/physics/TerrainDrag.cpp

+    }
+    amrex::Real value1, value2, value3;
+    while (file >> value1 >> value2 >> value3) {
+        m_xterrain.push_back(value1);


so if you use Array, these become nice loops.

marchdf · 2024-07-08T16:10:38Z

amr-wind/physics/TerrainDrag.cpp

+        auto& terrainz0 = m_terrainz0(level);
+        auto& drag = m_terrainDrag(level);
+        // copy terrain data to gpu
+        amrex::Gpu::DeviceVector<amrex::Real> gpu_xterrain(m_xterrain.size());


gpu -> device

marchdf · 2024-07-08T16:11:44Z

amr-wind/physics/TerrainDrag.cpp

+            amrex::ParallelFor(
+                vbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept {
+                    // compute the source term
+                    const amrex::Real x1 = prob_lo[0] + (i + 0.5) * dx[0];


call these x/y/z? or number from 0? Or stuff into gpuarray.

marchdf · 2024-07-08T16:12:33Z

amr-wind/physics/TerrainDrag.cpp

+                            break;
+                        }
+                    }
+                    int turnOn = (x3 <= terrainHt) ? 1 : 0;


just do levelBlanking(i, j, k, 0) = static_cast(x3 <= terrainHt)

… into hgopalan-terrain

hgopalan · 2024-07-09T19:06:55Z

Have addressed the comments

marchdf · 2024-07-09T20:35:32Z

There's a stray update to the amrex submods that shouldn't be there. Can you remove that?

marchdf · 2024-07-09T20:30:01Z

amr-wind/equation_systems/icns/source_terms/DragForcing.H

+    const CFDSim& m_sim;
+    const amrex::AmrCore& m_mesh;
+    const Field& m_velocity;
+    amrex::Gpu::DeviceVector<amrex::Real> device_vel_ht;


need the m_ prefix on these 2

marchdf · 2024-07-09T20:31:45Z

amr-wind/equation_systems/icns/source_terms/DragForcing.cpp

+    const auto& dx = geom.CellSize();
+    const auto& prob_lo = geom.ProbLoArray();
+    const auto& prob_hi = geom.ProbHiArray();
+    const amrex::Real device_drag = m_drag;


sorry to be a pain. But you can remove all the device_ prefix here. That's the nice thing about using the m_ pfx in the first place. There's a bunch of places like that here.

marchdf · 2024-07-09T20:32:44Z

amr-wind/equation_systems/temperature/source_terms/DragTempForcing.cpp

+    const auto& geom_vec = m_mesh.Geom();
+    const auto& geom = geom_vec[lev];
+    const auto& dx = geom.CellSize();
+    const amrex::Real gpu_drag = m_drag;


similar to above. Just call this drag etc.

marchdf · 2024-07-09T20:33:59Z

amr-wind/physics/TerrainDrag.cpp

+        m_z0rough.push_back(value3);
+    }
+    file1.close();
+    m_sim.io_manager().register_output_int_var("terrainDrag");


I would use snake case for the output variables. It keeps things consistent.

marchdf · 2024-07-09T20:36:20Z

test/test_files/terrainbox/terrainBox.inp

@@ -0,0 +1,63 @@
+#¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨#


I think this case/file should be named terrain_box

Agreed that this should be named terrain_box. But more importantly, the input file must have the same name as the directory in order for it to run with ctest. When you update the directory and input file names, make sure you update the CMakeLists.txt too

marchdf · 2024-07-17T14:22:06Z

tools/terrain/NREL.yaml

@@ -0,0 +1,29 @@
+solver: "amrWind"


We should remove these python scripts. There are no tests for these and we don't support python in the CI so these are liable to drift. We need to think of a better place to put these things. Right now it feels like they would belong to the amr-wind-frontend. Any thoughts on where this could live?

I am not sure of the best way. I can create a separate toolkit GitHub repo for terrain.

The setting up of a terrain case by hand is not feasible and this was the reason I had to create the python codes. The user typically creates the terrain and then zooms in to see that the region of interest is sufficient for him and again reruns the script to create the workflow. Also when we create the refinement boxes, the python script lets us keep it at the terrain height and reduces the marking of unnecessary cells.

This sure sounds useful. I would recommend a bunch of hardening: linting, proper packaging (with clear dependencies, python versions, etc), scripting (users don't edit the scripts, they pass arguments to the scripts).

The amr-wind-frontend could be a good place since that's what a lot of people are using to set up complex cases. I would hesitate to create yet another repo...

… GPU errors. Added a test case to show that IB does not affect the MOST BC.

hgopalan added 24 commits April 3, 2024 19:49

Initial Commit of Terrain Model

d3087d5

Update for clang format consistency

69320fa

Coding modifications for consistency

ec4086f

Clang Formatting

fc3724a

GPU Compatibility

7797847

Clang format

c3ad937

Code update for GPU compatibility

bd87c48

Update for Clang Format

25ffb13

GPU compatiblity for Drag Force

1b1b612

GPU Compatbility for Vertical Location

b735987

GPU compatibility

7fca03d

GPU Compatible Local Copy

903ae8e

Clang

201602c

GPU Updated

6055ba7

Clang

7e74c79

GPU Compatibility

425d455

Merge branch 'Exawind:main' into hgopalan-terrain

cd03664

Merge branch 'Exawind:main' into hgopalan-terrain

5b401df

Merge branch 'Exawind:main' into hgopalan-terrain

9dccb0d

Merge branch 'Exawind:main' into hgopalan-terrain

88dfaf9

Merge branch 'Exawind:main' into hgopalan-terrain

d871f53

This update modifies the code-base to make it more convenient to make…

db46f4a

… changes.

Python Utilities to create the terrain file

b225703

Updates for the clang format

383ce04

hgopalan requested review from mbkuhn, marchdf, moprak-nrel and mchurchf July 2, 2024 18:06

hgopalan added 2 commits July 2, 2024 12:16

Removing the unused variables for Temperature

2e7d011

Formatting

673d25f

Added a terrain case.

63968e3

moprak-nrel requested changes Jul 8, 2024

View reviewed changes

marchdf requested changes Jul 8, 2024

View reviewed changes

hgopalan added 6 commits July 8, 2024 17:26

Commit to address some nomenclature changes and moving python files

2473677

Merge branch 'main' into hgopalan-terrain

0e679d4

Merge branch 'hgopalan-terrain' of https://github.com/hgopalan/amr-wind…

c7f463b

… into hgopalan-terrain

Clang Formatting

4ecbf15

Fixing CMake Issue

02505aa

Coupled testing

f81b41a

hgopalan mentioned this pull request Jul 9, 2024

Calculation of flow field in complex terrain #1114

Closed

hgopalan added 5 commits July 9, 2024 09:04

Merge branch 'main' into hgopalan-terrain

a481179

Merge branch 'main' into hgopalan-terrain

cd73d96

Merge branch 'Exawind:main' into hgopalan-terrain

1693c6a

Modification forc converting blanking to int field

64c5d5c

Merge branch 'hgopalan-terrain' of https://github.com/hgopalan/amr-wind…

7a54646

… into hgopalan-terrain

marchdf requested changes Jul 9, 2024

View reviewed changes

hgopalan added 8 commits July 9, 2024 15:12

Updated for consistency and formatting

cd01d02

Clang formatting

18ee504

Formatting

13557f0

Merge branch 'Exawind:main' into hgopalan-terrain

fbc496e

Merge branch 'Exawind:main' into hgopalan-terrain

7da4b68

Rename terrainBox.inp to terrainbox.inp

64c285b

Merge branch 'Exawind:main' into hgopalan-terrain

ef6194b

Fixing terrain_box

e8e50eb

marchdf requested changes Jul 17, 2024

View reviewed changes

hgopalan and others added 2 commits July 17, 2024 11:14

Removed terrain tools from utilities. Modified forcing terms to avoid…

e005a1d

… GPU errors. Added a test case to show that IB does not affect the MOST BC.

Merge branch 'main' into hgopalan-terrain

0254ceb

		@@ -0,0 +1,6 @@
		This is the sample terrain case. You need to run "python3.XX backendinterface.py NREL.yaml" in the tools/terrain folder.

		@@ -0,0 +1,6 @@
		This is a simple test case to verify the terrain setup. A box is created and a freestream velocity

		@@ -0,0 +1,414 @@
		def SRTM_Converter(outputDir,refLat,refLon,refHeight,left,right,bottom,top):

		@@ -0,0 +1,63 @@
		#¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨#

Terrain-Aware Immersed Forcing Method #1115

Are you sure you want to change the base?

Terrain-Aware Immersed Forcing Method #1115

Conversation

hgopalan commented Jul 2, 2024 • edited Loading

Summary

Pull request type

Checklist

Additional background

moprak-nrel left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

marchdf left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

hgopalan commented Jul 9, 2024

marchdf commented Jul 9, 2024

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

hgopalan commented Jul 2, 2024 •

edited

Loading