Add active scalar in ConstraintIB model

[Shuang Shan]

Hi all, 

I am a new IBAMR user and have been working with the code for a few months. My PhD project is on simulation of particles in stratified fluids. In order to achieve this within IBAMR, I started from the ConstraintIB falling_sphere example (ibamr/IBAMR/examples/ConstraintIB/falling_sphere) and the navier_stokes/ex6 to incorporate an active scalar to the code. I then tried to combine these two codes, by registering the AdvDiffSemiImplicitHierarchyIntegrator to the Navier-Stokes solver in the ConstraintIB falling_sphere example following the approach used in navier_stokes/ex6. I have included the code link for clarity (please see https://drive.google.com/file/d/1dW2qYA-80oiTt568KYh2Glx2Acpz9YM0/view?usp=drive_link; https://drive.google.com/file/d/1Pv7iABfC6kM0BljqSKGAK2q9EcCGzF2t/view?usp=drive_link).

However, the temperature field does not react to motion of the particle (see images), which suggests that somehow the advection velocity is not passed to the scalar field integrator. Does anyone have an explanation of what's going on, and more importantly on how to solve this?

Many thanks and have a great week.

Shuang Shan

[Boyce Griffith]

Hmmm, it looks like your code is set up similarly to navier_stokes/ex6. Have you confirmed that coupling between the fluid solver and the concentration field is working correctly in ex6? There are also some IB examples that include this kind of coupling.

Can you try adding a few flags: “-adv_diff_ksp_rtol 1.0e-12 -adv_diff_ksp_min_it 1 -adv_diff_ksp_converged_reason” and send results from a time step or two towards the middle of the computation?

— Boyce

<temperature field 0s.png><temperature field 2.18s.png><velocity field 2.18s.png>

--
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/35e82c32-eb61-4424-8c43-76f19e713ea1n%40googlegroups.com.
<temperature field 2.18s.png><temperature field 0s.png><velocity field 2.18s.png>

[Shuang Shan]

Dear Prof. Boyce Griffith

Many thanks for your reply.

Yes, I confirm that coupling between fluid solver and concentration field works correctly by simulating a lock-exchange flow.

I tried to add the flags: “-adv_diff_ksp_rtol 1.0e-12 -adv_diff_ksp_min_it 1 -adv_diff_ksp_converged_reason” and some results in the middle of simulation are below:

The number of Linear adv_diff_0_ solve converged due to CONVERGED_RTOL iterations keeps 2 during whole simulation.

Warm wishes,

Shuang Shan

[Shuang Shan]

Dear Prof. Boyce Griffith

Many thanks for your useful suggestions, I have managed to add scalar in ConstraintIB/falling_sphere by referring IBAMR/examples/multiphase_flow/ex8.

Through reading the article on ConstraintIB method by Prof. Bhalla in 2013, I got that this falling_sphere case is a one-way coupling simulation, where the sphere is subject to 

a net body force: −(ρ_sphere​−ρ_fluid​)/ρ_fluid​⋅9.81, and fluid is not subject to any body force. Currently, I aim to implement a two-way coupling simulation, where the sphere is subject to 

a body force: −ρ_sphere​​/ ρ_fluid⋅9.81, while the fluid is subject to a Boussinesq body force (GAMMA*T, T is the scalar concentration) like navier_stokes/ex6. Given these objectives, I would like to ask whether the ConstraintIB/falling_sphere is still a suitable case for such modifications? If so, would it be feasible to directly modify ConstraintIB/falling_sphere/ForceProjector.cpp by using a levelset function to define a mask distinguishing between fluid and solid regions—applying a body force of GAMMA * T in the fluid, while setting the force projector for the solid region in input3d as follows:

ForceProjector {
  rho_fluid = RHO
  rho_solid = rho_fluid
  gravitational_constant = 0.0, -9.81*ρ_sphere​​/ ρ_fluid, 0.0
}

Would this approach be appropriate for implementing the desired two-way coupling? 

Besides, I would also like to impose a homogeneous Neumann boundary condition for the scalar field at the fluid-solid interface. Would it be sufficient to define a level set variable and apply the following command as in adv_diff/ex8:

brinkman_adv_diff->registerHomogeneousBC(Temperature, levelset, "Neumann", indicator_function_type, num_of_interface_cells, eta);

to achieve this?  Could the coupling with ConstraintIB pose any challenges or complications in setting this boundary condition?

As always, any help, advice, or recommendations would be very much appreciated.

Warm wishes,

Shuang Shan

[Boyce Griffith]

On May 6, 2025, at 11:46 AM, Shuang Shan <shanshu...@gmail.com> wrote:

Dear Prof. Boyce Griffith

Many thanks for your useful suggestions, I have managed to add scalar in ConstraintIB/falling_sphere by referring IBAMR/examples/multiphase_flow/ex8.

Great. What did you need to do to fix this? The code looked OK to me, and so I was clearly missing something.

Through reading the article on ConstraintIB method by Prof. Bhalla in 2013, I got that this falling_sphere case is a one-way coupling simulation, where the sphere is subject to 

a net body force: −(ρ_sphere​−ρ_fluid​)/ρ_fluid​⋅9.81, and fluid is not subject to any body force. Currently, I aim to implement a two-way coupling simulation, where the sphere is subject to 

a body force: −ρ_sphere​​/ ρ_fluid⋅9.81, while the fluid is subject to a Boussinesq body force (GAMMA*T, T is the scalar concentration) like navier_stokes/ex6.

No, unless I am misremembering, all of the cases described in Amneet’s 2013 paper are two-way coupled simulations. They are not using a Boussinesq-like formulation. The reason that they are able to use a constant-coefficient solver for the incompressible Navier-Stokes equations instead of a variable-density solver is because we are tracking the “excess” mass density associated with the solid in Lagrangian form. (If rho_s is the solid density and rho_f is the fluid density, then rho_s~ = rho_s - rho_f is the excess solid mass density. This quantity needs to be positive for the ConstraintIB formulation implemented in IBAMR to work.)

Given these objectives, I would like to ask whether the ConstraintIB/falling_sphere is still a suitable case for such modifications? If so, would it be feasible to directly modify ConstraintIB/falling_sphere/ForceProjector.cpp by using a levelset function to define a mask distinguishing between fluid and solid regions—applying a body force of GAMMA * T in the fluid, while setting the force projector for the solid region in input3d as follows:

ForceProjector {
  rho_fluid = RHO
  rho_solid = rho_fluid
  gravitational_constant = 0.0, -9.81*ρ_sphere​​/ ρ_fluid, 0.0
}

Would this approach be appropriate for implementing the desired two-way coupling? 

It can work — and in fact, Amneet and others have already implemented something like this in IBAMR as IBLevelSetMethod.

Besides, I would also like to impose a homogeneous Neumann boundary condition for the scalar field at the fluid-solid interface. Would it be sufficient to define a level set variable and apply the following command as in adv_diff/ex8:

brinkman_adv_diff->registerHomogeneousBC(Temperature, levelset, "Neumann", indicator_function_type, num_of_interface_cells, eta);

to achieve this?  Could the coupling with ConstraintIB pose any challenges or complications in setting this boundary condition?

I think this will do what you want — and I would probably use it with IBLevelSetMethod to have a consistent representation of the fluid-solid interface.

— Boyce

To view this discussion visit https://groups.google.com/d/msgid/ibamr-users/ef1976e1-feb9-40b7-af3a-a3387404d07dn%40googlegroups.com.

[Shuang Shan]

Dear Prof. Griffith

Many thanks for your suggestions. Sorry for the late reply.

The problem was with the input file. The initial condition for Temperature should be TemperatureInitialConditions{function = "(X_0 - X_com)^2 + (X_1 - Y_com)^2 + (X_2 - Z_com)^2 <= R^2 ? 0.0 : 8+625*X_1" // X_com, Y_com and Z_com are the particle center x, y and z coordinate} instead of function = "8+625*X_1". When I changed the initial condition, the problem was solved. I had followed the method of adding adv_diff in IBFE before, but later I realized that the solid part of IBFE is a mesh of a specific shape, while ConstraintIB is only composed of Lagrangian points.

When I try to study the physics of the falling sphere in the passive scalar field, I ran into two problems: (I have included the code link for clarity, please see https://drive.google.com/file/d/1NAhBpUxFwPArIBzaSvdvTLvHR73UszXk/view?usp=sharing   https://drive.google.com/file/d/1dW2qYA-80oiTt568KYh2Glx2Acpz9YM0/view?usp=drive_link, it is almost the same as the initial post in this thread just modify the TemperatureInitialConditions in the input file):

1. Scalar field error when increasing the resolution

I increase the resolution from CartesianGeometry {domain_boxes = [ (0,0,0) , (7, 31, 7) ]} to CartesianGeometry {domain_boxes = [ (0,0,0) , (15, 63, 15) ]}. Other parameters keeps unchanged (MAX_LEVELS = 3, level_1 = 4, 4, 4 // vector ratio to next coarser level, level_2 = 2, 2, 2). Plots 1 and 2 attached are velocity and scalar field at time = 0.4s in low and high resolution. In low resolution, both velocity and scalar fields overlap well with the Lagrangian markers (Plot 1).However, the scalar field obviously doesn't overlap with Lagrangian markers in high resolution (Plot 2). I have increased the resolution of the Lagrangian markers accordingly (sphere3d.vertex). I guess the problem may lie in the intermittent fragmented and disconnected finest level meshes, as shown in Plot 3. Does anyone have an explanation of what's going on, and more importantly on how to solve this? In addition, can I manually expand area of the Subset level 1 grid (the green region in Plot 3)? I hope the moving ball to always be at the finest level. Thank you so much.

2. HDF5 error during the parallel simulations:

When I run the built-in examples (e.g., IBFE, navier_stokes, ConstraintIB/falling_sphere, eel2d) using mpirun -np <processor number> ./main3d input3d command, everything works fine. However, the parallel simulation for modified falling_sphere with adv_diff integrator failed. It works well when <processor number> = 1. Here is the error: 

IBStandardInitializer:  Reading from input files.
  base filename: sphere3d
  assigned to level 1 of the Cartesian grid patch hierarchy

 ++++++++++++++++  STRUCTURE NO. 0  ++++++++++++++++++++++++++


 VOLUME OF THE MATERIAL ELEMENT           = 1.5625e-11
 VOLUME OF THE STRUCTURE                  = 6.49062e-08
 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

IBStandardInitializer:  Deallocating initialization data.
IBStandardInitializer:  Deallocating initialization data.


Writing visualization files...

HDF5-DIAG: Error detected in HDF5 (1.10.6) thread 0:
  #000: H5F.c line 509 in H5Fopen(): unable to open file
    major: File accessibilty
    minor: Unable to open file
  #001: H5Fint.c line 1498 in H5F_open(): unable to open file: time = Wed May 21 13:27:31 2025
, name = 'viz_FS/visit_dump.00000/processor_cluster.00001.samrai', tent_flags = 1
    major: File accessibilty
    minor: Unable to open file
  #002: H5FD.c line 734 in H5FD_open(): open failed
    major: Virtual File Layer
    minor: Unable to initialize object
  #003: H5FDsec2.c line 346 in H5FD_sec2_open(): unable to open file: name = 'viz_FS/visit_dump.00000/processor_cluster.00001.samrai', errno = 2, error message = 'No such file or directory', flags = 1, o_flags =                           2
    major: File accessibilty
    minor: Unable to open file
P=00003:Program abort called in file ``../../../../SAMRAI-2.4.4/source/toolbox/restartdb/HDFDatabase.C'' at line 2427
P=00003:ERROR MESSAGE:
P=00003:Unable to open HDF5 file viz_FS/visit_dump.00000/processor_cluster.00001.samrai
P=00003:
HDF5-DIAG: Error detected in HDF5 (1.10.6) thread 0:
  #000: H5F.c line 509 in H5Fopen(): unable to open file
    major: File accessibilty
    minor: Unable to open file
  #001: H5Fint.c line 1498 in H5F_open(): unable to open file: time = Wed May 21 13:27:31 2025
, name = 'viz_FS/visit_dump.00000/processor_cluster.00000.samrai', tent_flags = 1
    major: File accessibilty
    minor: Unable to open file
  #002: H5FD.c line 734 in H5FD_open(): open failed
    major: HDF5-DIAG: Error detecVirtual File Layer
    minor: Unable to initialize object
  #003: H5FDsec2.c line 346 in H5FD_sec2_open(): unable to open file: name = 'viz_FS/visit_dump.00000/processor_cluster.00000.samrai', errno = 2, error message = 'No such file or directory', flags = 1, o_flags =                           2
    major: HDF5-DIAG: Error detected in HDF5 (1.10.6) thread 0:
  #ted in HDF5 (1.10.6) thread 0:
  #000: H5F.c line 509 in H5Fopen(): unable to open fileFile accessibilty
    minor: Unable to open file

    major: File accessibilty
    minor: Unable to open file
  #001: H5Fint.c line 1498 in H5F_open(): unable to open file: time = Wed May 21 13:27:31 2025
, name = 'viz_FS/visit_dump.00000/processor_cluster.00003.samrai', tent_flags = 1
    major: File accessibilty
    minor:000: H5F.c line 509 in H5Fopen(): unable to open file
    major: File accessibilty
    minor: Unable to open file
  #001: H5Fint.c line 1498 P=00001:Program abort called in file ``../../../../SAMRAI-2.4.4/source/toolbox/restartdb/HDFDatabase.C'' at line 2427
P=00001:ERROR MESSAGE:
P=00001:Unable to open HDF5 file viz_FS/visit_dump.00000/processor_cluster.00000.samrai
P=00001:
in H5F_open(): unable to open file: time = Wed May 21 13:27:31 2025
, name = 'viz_FS/visit_dump.00000/processor_cluster.00002.samrai', tent_flags = 1
    major: File accessibilty
    minor: Unable to open file
  #002: H5FD.c line 734 in H5FD_open(): open failed
    major: Virtual File Layer
    minor: Unable to initialize object
 Unable to open file
  #002: H5FD.c line 734 in H5FD_open(): open failed
    major: Virtual File Layer
    minor: Unable to initialize object
  #003: H5FDsec2.c line 346 in H5FD_sec2_open(): unable to open file: name = 'viz_FS/visit_dump.00000/processor_cluster.00003.samrai', errno = 2, error message = 'No such file or directory', flags = 1, o_flags =                           2
    major: File accessibilty
    minor:   #003: H5FDsec2.c line 346 in H5FD_sec2_open(): unable to open file: name = 'viz_FS/visit_dump.00000/processor_cluster.00002.samrai', errno = 2, error message = 'No such file or directory', flags = 1                          , o_flags = 2
    major: File accessibilty
    minor: Unable to open file
Unable to open file
P=00007:Program abort called in file ``../../../../SAMRAI-2.4.4/source/toolbox/restartdb/HDFDatabase.C'' at line 2427
P=00007:ERROR MESSAGE:
P=00007:Unable to open HDF5 file viz_FS/visit_dump.00000/processor_cluster.00003.samrai
P=00007:
P=00005:Program abort called in file ``../../../../SAMRAI-2.4.4/source/toolbox/restartdb/HDFDatabase.C'' at line 2427
P=00005:ERROR MESSAGE:
P=00005:Unable to open HDF5 file viz_FS/visit_dump.00000/processor_cluster.00002.samrai
P=00005:
--------------------------------------------------------------------------
MPI_ABORT was invoked on rank 7 in communicator MPI_COMM_WORLD
with errorcode -1.

NOTE: invoking MPI_ABORT causes Open MPI to kill all MPI processes.
You may or may not see output from other processes, depending on
exactly when Open MPI kills them.
--------------------------------------------------------------------------
[login-e:3093925] 3 more processes have sent help message help-mpi-api.txt / mpi-abort
[login-e:3093925] Set MCA parameter "orte_base_help_aggregate" to 0 to see all help / error messages

As always, any help, advice, or recommendations would be very much appreciated. Have a nice day!

Warm wishes,

Shuang Shan

[Boyce Griffith]

On May 21, 2025, at 7:45 AM, Shuang Shan <shanshu...@gmail.com> wrote:

Dear Prof. Griffith

Many thanks for your suggestions. Sorry for the late reply.

The problem was with the input file. The initial condition for Temperature should be TemperatureInitialConditions{function = "(X_0 - X_com)^2 + (X_1 - Y_com)^2 + (X_2 - Z_com)^2 <= R^2 ? 0.0 : 8+625*X_1" // X_com, Y_com and Z_com are the particle center x, y and z coordinate} instead of function = "8+625*X_1". When I changed the initial condition, the problem was solved. I had followed the method of adding adv_diff in IBFE before, but later I realized that the solid part of IBFE is a mesh of a specific shape, while ConstraintIB is only composed of Lagrangian points.

When I try to study the physics of the falling sphere in the passive scalar field, I ran into two problems: (I have included the code link for clarity, please see https://drive.google.com/file/d/1NAhBpUxFwPArIBzaSvdvTLvHR73UszXk/view?usp=sharing   https://drive.google.com/file/d/1dW2qYA-80oiTt568KYh2Glx2Acpz9YM0/view?usp=drive_link, it is almost the same as the initial post in this thread just modify the TemperatureInitialConditions in the input file):

1. Scalar field error when increasing the resolution

I increase the resolution from CartesianGeometry {domain_boxes = [ (0,0,0) , (7, 31, 7) ]} to CartesianGeometry {domain_boxes = [ (0,0,0) , (15, 63, 15) ]}. Other parameters keeps unchanged (MAX_LEVELS = 3, level_1 = 4, 4, 4 // vector ratio to next coarser level, level_2 = 2, 2, 2). Plots 1 and 2 attached are velocity and scalar field at time = 0.4s in low and high resolution. In low resolution, both velocity and scalar fields overlap well with the Lagrangian markers (Plot 1).However, the scalar field obviously doesn't overlap with Lagrangian markers in high resolution (Plot 2). I have increased the resolution of the Lagrangian markers accordingly (sphere3d.vertex). I guess the problem may lie in the intermittent fragmented and disconnected finest level meshes, as shown in Plot 3. Does anyone have an explanation of what's going on, and more importantly on how to solve this? In addition, can I manually expand area of the Subset level 1 grid (the green region in Plot 3)? I hope the moving ball to always be at the finest level. Thank you so much.<Plot 1.png><Plot 2.png><Plot 3.png>

What does the field behind the structure look like? Are those values “leaking out” through the fluid-structure interface? Is this just passive transport (advection), or is it an advection-diffusion problem where you may need a boundary condition along the fluid-solid interface?

You might try reducing the time step size, which will reduce the amount of slip between the fluid and the structure in ConstraintIBMethod.

2. HDF5 error during the parallel simulations:

When I run the built-in examples (e.g., IBFE, navier_stokes, ConstraintIB/falling_sphere, eel2d) using mpirun -np <processor number> ./main3d input3d command, everything works fine. However, the parallel simulation for modified falling_sphere with adv_diff integrator failed. It works well when <processor number> = 1. Here is the error: 

These kinds of errors are often related to an underlying issue with the file system — permissions are not set correctly, quota exceeded, or something like that.

To view this discussion visit https://groups.google.com/d/msgid/ibamr-users/51b78f05-c56c-471d-8b51-acb28b332445n%40googlegroups.com.
<Plot 2.png><Plot 1.png><Plot 3.png>