Spatially varying viscosity that is a function of a state variable?

[dcoveney]

Hi,

I am looking to run multiphase incompressible simulations in IAMR.

As per many commonly used formulations, the idea is for the two fluids to be indicated by a passive advected scalar (e.g. the volume fraction) that is used to compute the viscosity:

where the scalar varies between 0 and 1.

As far as I can see from the documentation/comments, this is supported in IAMR, through the calcViscosity function in the NavierStokes class. From what I understand, the entire MultiFab referenced by the pointer visc[dir] is filled with a constant viscosity coefficient ParmParsed in from the inputs file (vel_visc_coef), as below:

for (int dir=0; dir<AMREX_SPACEDIM; dir++)
 {
     visc[dir]->setVal(visc_coef[Xvel], 0, visc[dir]->nComp(), visc[dir]->nGrow());          
 }

What I want to do is rewrite this so that I can dereference the MultiFab pointer and extract the data using the array() functionality. Additionally, I would like the viscosity to be calculated using the Tracer variable in theMultiFab containing the state data. The code I have written for this is shown below (which would be inside calcViscosity):

 ParmParse pp("ns");
pp.query("muG", muG);
pp.query("muL", muL);

MultiFab&   S_new    = get_new_data(State_Type);
auto whichTime = which_time(State_Type,time);
BL_ASSERT(whichTime == AmrOldTime || whichTime == AmrNewTime);
auto visc = (whichTime == AmrOldTime ? viscn : viscnp1);

for(int dir = 0; dir <AMREX_SPACEDIM; dir++)
{

     MultiFab& viscMF = *visc[dir];
     #ifdef _OPENMP
     #pragma omp parallel
     #endif
     for (MFIter viscnewmfi(viscMF,true); viscnewmfi.isValid(); ++viscnewmfi)
     {

         const Box& vbx = viscnewmfi.tilebox();
         FArrayBox& Sfab = S_new[viscnewmfi];
         FArrayBox& viscFab = viscMF[viscnewmfi];
         auto lo = lbound(vbx);
         auto hi = ubound(vbx);
         Array4<Real>const& state = Sfab.array();
         Array4<Real>const& viscArray = viscFab.array();
         for(int i = lo.x; i <= hi.x; i++)
          {
              for(int j = lo.y; j <= hi.y; j++)
              {
                  for(int k = lo.z; k <= hi.z; k++)
                  {
                     viscArray(i,j,k,0) = muG*state(i,j,k,Tracer) + muL*(1.0-state(i,j,k,Tracer));
                  }
               }
           }
        }
}

My questions are:

1. Would the MultiFab referenced by visc[dir] necessarily have the same spatial indices as S_new in this case?
2. Is there native support for variable viscosity or a simpler solution to this than I have presented in the above code that I haven't noticed? I saw in the April 2020 release notes that it says "no longer supporting constant mu" so I was wondering if there was, or if this comment refers to another aspect of the code such as Diffusion.cpp.
3. Assuming my code is correct, is there anything else I would have to change in the code? The only comment is in the code above calcViscosity:

   // Functions for calculating the variable viscosity and diffusivity.
   // These default to setting the variable viscosity and diffusivity arrays
   // to the values in visc_coef and diff_coef.  These functions would
   // need to be replaced in any class derived from NavierStokes that
   // wants variable coefficients.

   So I'm assuming not.

Thank you for your support, I look forward to your response.

[drummerdoc]

Couple of quick comments. Support for variable transport coefficients refers to the fact that they are computed and stored in the code as MultiFab's rather than as a scalar. Your approach overall is correct. In detail, there are a few things to fix. Firstly, the state data (typically) lives at the cell centroid, yet the transport coefficients are stored at the centroid of the faces. In other applications that use state-dependent transport, different strategies are used...such as evaluating on cell centroids and moving the faces. That is likely what you want to do here. Finally, you might step through the Tutorials/GPU/CNS code to for examples to see how to wrap the ijk access code you wrote into a lambda function, as it will make it more compatible across different hardware systems.

This was a quick answer, I realize and maybe a bit cryptic. Ask more questions to clarify.

[dcoveney]

Thank you for your prompt response. I will investigate the wrapper for the ijk loops, thanks for referring me to that tutorial.

As far as moving the viscosity coefficient to the cell faces, is there a simple way of doing this in IAMR? Are there any examples of state-dependent transport in tutorials? Thanks

[drummerdoc]

PeleLM does this. In that code, there's calcViscosity which computes mu at cell centers, and then getViscosity which moves it to face centers. There are functions in AMReX to move CC data to face centroids (look in the Src/EB folder for utility code), but those do not properly account for the location of Dirichlet boundary data. In PeleLM, there is an attempt to do this more correctly, but the code in there now is not properly EB aware yet. So, there are a few different places to look and assemble. Note that if the viscosity is variable, you need to be using the tensor solver - should be obvious, but figured I should underscore that one.

[emotheau]

Just a side note: in IAMR the viscosity is now only solved with the tensor solver. So even a constant viscosity is treated as a variable. It was not the case in the past, when we had two separate solver for either constant or variable viscosity.

Also in IAMR, if you look at the development branch, in the routine calcViscosity I have added a feature to compute an additional modeled turbulent viscosity for LES. Maybe you could follow a similar way to adapt the viscosity as you wish.

[asalmgren]

Emmanuel -- why did we eliminate that option? We believe that div(u) = 0 and constant mu should enable us to solve separately for the different components which might be faster -- is there a reason not to allow that option?

On Fri, Apr 24, 2020 at 3:47 PM Emmanuel Motheau notifications@github.com wrote:

Just a side note: in IAMR the viscosity is now only solved with the tensor solver. So even a constant viscosity is treated as a variable. It was not the case in the past, when we had two separate solver for either constant or variable viscosity.

— You are receiving this because you are subscribed to this thread. Reply to this email directly, view it on GitHub https://github.com/AMReX-Codes/IAMR/issues/39#issuecomment-619269745, or unsubscribe https://github.com/notifications/unsubscribe-auth/ACRE6YTRTWLVEZ6BIEJM3P3ROIJJFANCNFSM4MQLHDRA .

-- Ann Almgren Senior Scientist; CCSE Group Lead

[drummerdoc]

that was probably my fault. I wanted to eliminate as much code as possible and there didn't seem to be a big call for such things...nobody cares about 2D or r-z either.

[emotheau]

I really don't remember why and who decided to remove the constant solver for viscosity. According to git, last year I had to put it back to validate EB implementation, but it was removed again later.

[dcoveney]

Thank you for your comments and guidance - I will look at the various options for development of a state-dependent viscosity, including PeleLM, EB functionality as well as potentially editing the turbulent viscosity feature.