Bringing back VOLP for ghost cell (TUM) #124

swenczowski · 2024-11-24T14:29:59Z

--- please, consider this PR a draft for discussion, improvements to follow, not tested ---

Dear all,

as encountered by @javierebb, mglet-base does not provide the field "VOLP", which provides the volume of intersected cells within the ghost cell immersed boundary framework.

This PR outlines an idea to recover this field, which is primarily used for postprocessing on our side. The following considerations have been made, such that I am slightly optimistic that all relevant locations in the code have been touched. As discussed beforehand, both areas and volume are given absolute in [L^2] or [L^3] and not in terms of fractions.

I would much appreciate your ideas and feedback on this very rough first draft.

Best regards,

Simon

hakostra

The significant comments are on some suggestions on how to improve the volp-computations.

A last small request:
Can you simultaneously rename the fields AU, AV, AW - that are now real areas with unit length squared - to AREAU, AREAV, AREAW? This will apply to the fields that are created in gc_mod.F90 and accessed in gc_flowstencils_mod.F90 - not gc_blockbp_mod.F90 - this should be left untouched!

This establish the practice that field AU = non-dimensional field with fraction of open face area - only existing locally in blocking and AREAU = real, physical area with unit length squared used during time-integration.

It is enough to rename the fields them selves, you do not need to rename the pointer arrays unless you see particular needs, i.e. change:

CALL set_field("AU")
CALL get_fieldptr(au, "AU")

to

CALL set_field("AREAU", istag=1, units=[0, 2, 0, 0, 0, 0, 0])
CALL get_fieldptr(au, "AREAU")

I think you get this point as well.

hakostra · 2024-11-25T08:22:24Z

src/ib/gc_stencils_mod.F90


        ! Local variables
-        INTEGER(intk) :: imygrid, idx, igrid, ind, ip3
+        TYPE(field_t), POINTER :: ddx_f, ddy_f, ddz_f
+        INTEGER(intk) :: imygrid, idx, igrid, ind
        INTEGER(intk) :: aucells
        INTEGER(intk) :: kk, jj, ii, k, j, i
+        REAL(realk) :: area
        REAL(realk), POINTER, CONTIGUOUS :: au(:, :, :)
-
-        ! First set AU to value of BU
-        au_f%arr = bu%arr
-        av_f%arr = bv%arr
-        aw_f%arr = bw%arr
-
-        ! Now set AU, AV, AW
+        REAL(realk), POINTER, CONTIGUOUS :: bu(:, :, :)
+        REAL(realk), POINTER, CONTIGUOUS :: ddx(:), ddy(:), ddz(:)
+
+        CALL get_field(ddx_f, "DDX")
+        CALL get_field(ddy_f, "DDY")
+        CALL get_field(ddz_f, "DDZ")
+
+        ! Now set AU, AV, AW ( absolute areas in unit [L^2] )
        all_grids: DO imygrid = 1, nmygrids
+
            igrid = mygrids(imygrid)
            CALL get_mgdims(kk, jj, ii, igrid)
-            CALL get_ip3(ip3, igrid)

+            ! Getting the pressure cell side lengths
+            CALL ddx_f%get_ptr(ddx, igrid)
+            CALL ddy_f%get_ptr(ddy, igrid)
+            CALL ddz_f%get_ptr(ddz, igrid)
+
+            ! Area AU from BU
            CALL au_f%get_ptr(au, igrid)
+            CALL bu_f%get_ptr(bu, igrid)
+            DO i = 1, ii
+                DO j = 1, jj
+                    DO k = 1, kk
+                        area = ddz(k) * ddy(j)
+                        au(k, j, i) = bu(k, j, i) * area
+                    END DO
+                END DO
+            END DO
+            ! Refine for intersected cells
            aucells = SIZE(this%auind(imygrid)%arr)
            DO idx = 1, aucells
                ind = this%auind(imygrid)%arr(idx)
                CALL ind2sub(ind, k, j, i, kk, jj, ii)
-                au(k, j, i) = this%auvalue(imygrid)%arr(idx)
+                area = ddz(k) * ddy(j)
+                au(k, j, i) = this%auvalue(imygrid)%arr(idx) * area
            END DO

+            ! Area AV from BV
            CALL av_f%get_ptr(au, igrid)
+            CALL bv_f%get_ptr(bu, igrid)
+            DO i = 1, ii
+                DO j = 1, jj
+                    DO k = 1, kk
+                        area = ddx(i) * ddz(k)
+                        au(k, j, i) = bu(k, j, i) * area
+                    END DO
+                END DO
+            END DO
+            ! Refine for intersected cells
            aucells = SIZE(this%avind(imygrid)%arr)
            DO idx = 1, aucells
                ind = this%avind(imygrid)%arr(idx)
                CALL ind2sub(ind, k, j, i, kk, jj, ii)
-                au(k, j, i) = this%avvalue(imygrid)%arr(idx)
+                area = ddx(i) * ddz(k)
+                au(k, j, i) = this%avvalue(imygrid)%arr(idx) * area
            END DO

+            ! Area AW from BW
            CALL aw_f%get_ptr(au, igrid)
+            CALL bw_f%get_ptr(bu, igrid)
+            DO i = 1, ii
+                DO j = 1, jj
+                    DO k = 1, kk
+                        area = ddx(i) * ddy(j)
+                        au(k, j, i) = bu(k, j, i) * area
+                    END DO
+                END DO
+            END DO
+            ! Refine for intersected cells
            aucells = SIZE(this%awind(imygrid)%arr)
            DO idx = 1, aucells
                ind = this%awind(imygrid)%arr(idx)
                CALL ind2sub(ind, k, j, i, kk, jj, ii)
-                au(k, j, i) = this%awvalue(imygrid)%arr(idx)
+                area = ddx(i) * ddy(j) 
+                au(k, j, i) = this%awvalue(imygrid)%arr(idx) * area
            END DO

        END DO all_grids
    END SUBROUTINE get_auavaw


I wonder if it would be simpler and cleaner to keep the original code, constructing au, av aw as fractions in the range 0..1, and then just in the end have another loop:

DO i = 1, ii DO j = 1, jj DO k = 1, kk au(k, j, i) = ddy(j)*ddz(k)*au(k, j, i) av(k, j, i) = ddx(i)*ddz(k)*av(k, j, i) aw(k, j, i) = ddx(i)*ddy(j)*aw(k, j, i) END DO END DO END DO

What do you think?

Ok, can be done. Admittedly, I do not dislike the current form, as it shows the systematics more clearly. However, I have no strong preferences...

The you make a call - I will accept both.

hakostra · 2024-11-25T08:23:59Z

src/ib/gc_stencils_mod.F90

+            CALL get_fieldptr(x, "X", igrid)
+            CALL get_fieldptr(y, "Y", igrid)
+            CALL get_fieldptr(z, "Z", igrid)


Tip: XSTAG, YSTAG, ZSTAG also exist as fields - you no longer need to construct these on the fly.

hakostra · 2024-11-25T08:24:44Z

src/ib/gc_stencils_mod.F90

+        REAL(realk), INTENT(in) :: ddx(ii), ddy(jj), ddz(kk)
+        REAL(realk), INTENT(in) :: dx(ii), dy(jj), dz(kk)
+        REAL(realk), INTENT(in) :: x(ii), y(jj), z(kk)


Pass in XSTAG, YSTAG, ZSTAG, DDX, DDY, DDZ

Ohh, ok. Will be changed.

hakostra · 2024-11-25T08:26:16Z

src/ib/gc_stencils_mod.F90

+                        areax = au(k, j, i-1) - au(k, j, i)
+                        areay = av(k, j-1, i) - av(k, j, i)
+                        areaz = aw(k-1, j, i) - aw(k, j, i)


In most other stencil operations we have the higher index first and subtract the lower index

Suggested change

areax = au(k, j, i-1) - au(k, j, i)

areay = av(k, j-1, i) - av(k, j, i)

areaz = aw(k-1, j, i) - aw(k, j, i)

areax = au(k, j, i) - au(k, j, i-1)

areay = av(k, j, i) - av(k, j-1, i)

areaz = aw(k, j, i) - aw(k-1, j, i)

This will flip the sign of the areas, see also comment below.

Ok. Will be changed.

hakostra · 2024-11-25T08:29:52Z

src/ib/gc_stencils_mod.F90

+                        vpx = (-au(k, j, i-1) * (xstagm - xstag)) &
+                            + areax * (swx - xstag)


The difference between two xstag values are ddx: xstagm - xstag = -ddx.
(-au)*(-ddx) = au*ddx

Since we flipped the sign of areax, we also flip the sign of (swx - xstag).

So you end up with something like: vpx = (au(k, j, i-1)*ddx(i)) + areax*(xstag-swx)

This will be more numerically accurate (since we take one less difference between numbers of similar magnitude - subject to high round-off errors), easier to read and less code.

The same of course apply for the y- and z-directions as well - I guess you get the point.

Sorry, that also came to my mind yesterday evening -- of course, the difference between two staggered positions is the pressure cell side length...

When coding it, I was too focussed on a pure replication of the previous implementation... Will be changed.

swenczowski · 2024-12-01T18:36:14Z

Dear all,

this PR has been updated based on the preceding discussion. Thank you for the hints!

@hakostra, please, let me know if all points have been addressed.

@javierebb, please, already give the implementation a try.

Best regards,

Simon

hakostra · 2024-12-02T07:07:26Z

Hi Simon, now try ci/lint-fortran.py --fixall

hakostra · 2024-12-02T08:58:45Z

@swenczowski This looks good. Are your intentions now to apply the areau/v/w and volp to the scalar solver as discussed?

swenczowski · 2024-12-02T11:17:39Z

@swenczowski This looks good. Are your intentions now to apply the areau/v/w and volp to the scalar solver as discussed?

Yes, that will come now. Also, we are planning to add a volume source for the scalar field, which would also need the VOLP field.

swenczowski · 2024-12-02T12:42:48Z

If you agree, I suggest merging the content of this PR before we continue with the actual "usage" of VOLP for the scalar.
Would that be OK for you?

swenczowski requested a review from hakostra as a code owner November 24, 2024 14:29

hakostra requested changes Nov 25, 2024

View reviewed changes

swenczowski force-pushed the simon/volp branch from e7f9dd5 to b5bfb61 Compare December 1, 2024 18:33

VOLP and change from AU to AREAU as actual area (TUM)

a162bc6

swenczowski force-pushed the simon/volp branch from b5bfb61 to a162bc6 Compare December 2, 2024 11:15

hakostra approved these changes Dec 2, 2024

View reviewed changes

hakostra merged commit 398a527 into kmturbulenz:master Dec 2, 2024
1 check passed

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Bringing back VOLP for ghost cell (TUM) #124

Bringing back VOLP for ghost cell (TUM) #124

swenczowski commented Nov 24, 2024

hakostra left a comment

hakostra Nov 25, 2024

swenczowski Nov 25, 2024

hakostra Nov 25, 2024

hakostra Nov 25, 2024

hakostra Nov 25, 2024

swenczowski Nov 25, 2024

hakostra Nov 25, 2024

swenczowski Nov 25, 2024

hakostra Nov 25, 2024

swenczowski Nov 25, 2024

swenczowski commented Dec 1, 2024

hakostra commented Dec 2, 2024

hakostra commented Dec 2, 2024

swenczowski commented Dec 2, 2024

swenczowski commented Dec 2, 2024

		vpx = (-au(k, j, i-1) * (xstagm - xstag)) &
		+ areax * (swx - xstag)

Bringing back VOLP for ghost cell (TUM) #124

Bringing back VOLP for ghost cell (TUM) #124

Conversation

swenczowski commented Nov 24, 2024

hakostra 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

swenczowski commented Dec 1, 2024

hakostra commented Dec 2, 2024

hakostra commented Dec 2, 2024

swenczowski commented Dec 2, 2024

swenczowski commented Dec 2, 2024