Scalar field conservation issue during AMR regridding (AdvDiff + INSStaggered)

[Teo Lara]

Hello everyone!

I’m using ConstrainedIBMethod to model swimming bodies and would like to simultaneously advect a scalar field with the fluid flow. Based on navier_stokes/ex6, I have coupled AdvDiffSemiImplicitHierarchyIntegrator with INSStaggeredHierarchyIntegrator for this purpose.

However, I’ve noticed that during regriddings the total integral of the scalar field is not conserved, leading to unphysical jumps in the total scalar “mass.” My setup has no inflow/outflow of the scalar field, and the small amount of diffusion I include is not sufficient to explain these jumps. For phyiscal intuition, this scalar field represents some nutrient "food" field, whose total mass should be conserved. 

In the attached figure, I show the total amount of “food” (the scalar field) obtained by integrating over the finest patches of my entire simulation domain. Some variation is expected, but I am particularly concerned about the large jumps, which I believe are caused by regridding.

Question: Is there a way to enforce conservation of the scalar field during AMR regridding when coupling AdvDiffSemiImplicitHierarchyIntegrator to INSStaggeredHierarchyIntegrator?

Best,

Teo Lara

[Amneet Bhalla]

I think if you use equation for scalar phi of the form:

\partial (phi)/\partial t + div (u phi) = div (D grad phi) 

with no slip BC for velocity and homogeneous Neumann for phi, then the total phi will be conserved. This can be shown by taking a volume integral of the above PDE over the domain and using Gauss-Divergence theorem. 

--Amneet 

--
You received this message because you are subscribed to the Google Groups "IBAMR Users" group.
To unsubscribe from this group and stop receiving emails from it, send an email to ibamr-users...@googlegroups.com.
To view this discussion visit https://groups.google.com/d/msgid/ibamr-users/65b071cc-f2b5-471d-8bab-02521b9a1155n%40googlegroups.com.

[Aaron Barrett]

As Amneet has mentioned, depending on the boundary conditions, the total amount will be conserved across the entire domain. The regridding process uses conservative interpolation and coarsening operators, and I can confirm that regridding is conservative.

In the attached figure, I show the total amount of “food” (the scalar field) obtained by integrating over the finest patches of my entire simulation domain. Some variation is expected, but I am particularly concerned about the large jumps, which I believe are caused by regridding.

Are you integrating over the finest patches or the entire domain? The integral over the finest patches will not be conserved, so you may just be measuring the volume changes.

This snippet of code will print the total amount across the entire hierarchy that you can put in the main time loop.

        auto var_db = VariableDatabase<NDIM>::getDatabase();

        const int T_idx = var_db->mapVariableAndContextToIndex(T_var, adv_diff_integrator->getCurrentContext());

        HierarchyMathOps hier_math_ops("hier_math_ops", patch_hierarchy);

        HierarchyCellDataOpsReal<NDIM, double> hier_cc_ops(patch_hierarchy, 0, patch_hierarchy->getFinestLevelNumber());

        const double val = hier_cc_ops.integral(T_idx, hier_math_ops.getCellWeightPatchDescriptorIndex());

        pout << "At time " << loop_time << " the total value is " << val << "\n";

Aaron

---

From: ibamr...@googlegroups.com <ibamr...@googlegroups.com> on behalf of Amneet Bhalla <mail2...@gmail.com>
Sent: Friday, September 5, 2025 5:54 PM
To: ibamr...@googlegroups.com <ibamr...@googlegroups.com>
Subject: Re: [ibamr-users] Scalar field conservation issue during AMR regridding (AdvDiff + INSStaggered)

 

To view this discussion visit https://groups.google.com/d/msgid/ibamr-users/CAMETWJ00qAnALuOOt%3D17AYczrKzZXsM_no%3DbOh2MT52ZKXW4zw%40mail.gmail.com.

[Teo Lara]

Hi!

Thanks for the answers!

Yes, this is the exact equation I am using to govern the phi field: \partial (phi)/\partial t + div (u phi) = div (D grad phi) 
There is are no source/sink terms, I am using homogenous neumann bcs, there is no scalar flux outside of the simulation domain. I agree, in theory everything should be conserved. 

The "food" quantity that I plotted is obtained using hier_cc_ops.integral, like so:

            int C_idx = var_db->mapVariableAndContextToIndex(concentration_var, adv_diff_integrator->getCurrentContext());
            double total_food_in_simulation = hier_cc_ops.integral(C_idx, hier_math_ops.getCellWeightPatchDescriptorIndex())  ;

I thought this was an integral over the finest patches, but in any case, I'm not calculating the total mass of the scalar field in some other strange way, I am using HierarchyCellDataOpsReal as suggested in the previous response.
What led me to believe this was a regridding issue was that everything works fine without jumps when running on one level, but when running on 4 levels, things go sideways. I don't know what is going on, perhaps it has something to do with how I set up the advection integrator?

Any ideas would be very helpful!

Thanks!

Teo

[Amneet Bhalla]

Also make sure you are using the conservative form of advection in the input file, i.e., div (u phi) and not the non-conservative form, u. grad (phi). 

To view this discussion visit https://groups.google.com/d/msgid/ibamr-users/47da7669-3f24-495d-b2a1-9ba2ecbb2847n%40googlegroups.com.

--

--Amneet 

[Boyce Griffith]

On Sep 5, 2025, at 7:15 PM, Amneet Bhalla <mail2...@gmail.com> wrote:

Also make sure you are using the conservative form of advection in the input file, i.e., div (u phi) and not the non-conservative form, u. grad (phi). 

I agree that is important, but that shouldn’t impact regridding, should it…?

To view this discussion visit https://groups.google.com/d/msgid/ibamr-users/CAMETWJ1uHb_SGJE1FO6n-MiYxfa%3DktZE7cSFpjn9sPuMGNKk5g%40mail.gmail.com.

[Aaron Barrett]

I can't recreate this on navier_stokes/ex5 with forced regridding. Using advective form may result in some amount changes over time, but it shouldn't lead to any jumps.

A couple of suggestions:

The data ops classes can be created to run on different sets of level numbers. Make sure they are operating on the entire hierarchy, and not a subset.

Does this occur on 2 levels, or is it only 4?

Can you confirm that the jumps in concentrations occur after a regrid? If you are running with logging enabled, the log file should note when a regridding operation occurs.

---

From: ibamr...@googlegroups.com <ibamr...@googlegroups.com> on behalf of Boyce Griffith <boy...@gmail.com>
Sent: Friday, September 5, 2025 8:17 PM
To: IBAMR Users <ibamr...@googlegroups.com>

To view this discussion visit https://groups.google.com/d/msgid/ibamr-users/45AB0807-A68D-43FE-AE1F-826500D79C1D%40gmail.com.

[Teo Lara]

Update: This is fixed!!

turns out I was doing the integral incorrectly, I wasn't calling  

HierarchyCellDataOpsReal<NDIM, double> hier_cc_ops(patch_hierarchy, 0, patch_hierarchy->getFinestLevelNumber());

at everytimestep, I just did it once at the beginning. This caused what looked like to be problems during the regriddings, although in reality everything was always fine. 

Thanks everyone for the help!!