Specifying P_problem_coef in StaggeredStokesOperator class

[amneetb]

Currently we have setVelocityPoissonSpecifications() routine in StaggeredStokesOperator class to specify U problem coefficients; see

https://github.com/IBAMR/IBAMR/blob/master/include/ibamr/StaggeredStokesOperator.h#L99

For acoustic streaming flows, we compute $\nabla \cdot (\mathcal{D} \vec{u})$ as the low Mach equation instead of the usual $\nabla \cdot \vec{u}$ incompressibility/low-Mach equation. Here, $\mathcal{D}$ is a spatially varying diffusion-like coefficient.

It will be great to have setPressurePoissonSpecifications() routine in this class to specify P problem coefficients.

[boyceg]

Ultimately, are these specifications for a different set of equations than the Stokes equations? If so, it might be better to have a different class altogether for the problem coefficients.

A mistake that I made long ago is to piggyback off of the Poisson specification object provided by SAMRAI. I could have avoided a lot of confusion over the years if the specification objects actually corresponded to the equations that we are using...

I am starting to think that the PETSc approach of passing specifications around using void* contexts might be the best way to handle this kind of thing without being overly prescriptive in the interface.

[amneetb]

Ultimately, are these specifications for a different set of equations than the Stokes equations? If so, it might be better to have a different class altogether for the problem coefficients.

It is a generalized version of the Stokes equation where the difference arises in the (2,1) block when we have variable density. In that case, there is a constraint on the momentum $\nabla \cdot (\mathcal{\rho} \vec{u})$ rather than on velocity $\nabla \cdot \vec{u}$. The momentum equation is exactly the same. In other words, (1,1) and (1,2) blocks remain the same. The (2,2) block is also the same (0 in both cases).

A mistake that I made long ago is to piggyback off of the Poisson specification object provided by SAMRAI. I could have avoided a lot of confusion over the years if the specification objects actually corresponded to the equations that we are using...

I am starting to think that the PETSc approach of passing specifications around using void* contexts might be the best way to handle this kind of thing without being overly prescriptive in the interface.

Yes, that is a good idea. It avoids naming these objects as PoissonSpecifications, which may or may not be correct, depending upon the context. Can you or someone put a PR to show what you have in mind (in which routine is the void* passed) , and I can follow that style?

[boyceg]

Ultimately, are these specifications for a different set of equations than the Stokes equations? If so, it might be better to have a different class altogether for the problem coefficients.

It is a generalized version of the Stokes equation where the difference arises in the (2,1) block when we have variable density. In that case, there is a constraint on the momentum ∇⋅(ρu→) rather than on velocity ∇⋅u→. The momentum equation is exactly the same. In other words, (1,1) and (1,2) blocks remain the same. The (2,2) block is also the same (0 in both cases).

Is there a name attached to these equations? I guess they come from low Mach, but I think that the low Mach equations include other terms beyond the linear ones. I need to think about whether it makes sense to make Stokes the base class, a specialization, or just make this a separate class hierarchy.

A mistake that I made long ago is to piggyback off of the Poisson specification object provided by SAMRAI. I could have avoided a lot of confusion over the years if the specification objects actually corresponded to the equations that we are using... I am starting to think that the PETSc approach of passing specifications around using void* contexts might be the best way to handle this kind of thing without being overly prescriptive in the interface.

Yes, that is a good idea. It avoids naming these objects as PoissonSpecifications, which may or may not be correct, depending upon the context. Can you or someone put a PR to show what you have in mind (in which routine is the void* passed) , and I can follow that style?

@drwells & @abarret, what do you all think?

[abarret]

I don't see why we need to force a uniform problem specification interface for each type of operator. These operators are typically not interchangeable, and you'll have to look up how to specify each operator regardless. If we are going to change how we create operators, my vote is to rewrite GeneralOperator with RAII in mind, so that most parameters are set in the constructor and removing initializeOperatorState() and deallocateOperatorState(). We could then get rid of a lot of member functions which may or may not be used in child classes: GeneralOperator::setSolutionTime(), GeneralOperator::setTimeInterval(), etc. When setting up operators, it is often unclear which of these are actually used. However, this is a very invasive change.

If you want to enforce a uniform interface, one options is to add an empty ProblemSpecification struct. You could add functionality to GeneralOperator:

void GeneralOperator::setProblemSpecification(ProblemSpecification* spec);

Then each type of operator can have an associated struct, e.g. for the StaggeredStokesOperator you can have

struct StaggeredStokesSpec : public ProblemSpecification
{
double d_mu, d_rho;
}

and for VCStaggeredStokesOperator

struct VCStaggeredStokesSpec : public ProblemSpecification
{
int d_mu_idx, d_rho_idx;
}

[amneetb]

Where should the struct ProblemSpecification be defined? ibtk_utilities.h?

[boyceg]

@abarret what do you think? I think it seems like something that should get its own header.

[amneetb]

Now I have three coefficients, C for density, D for shear viscosity and E for bulk viscosity. PoissonSpecifications is not cutting out anymore. We need to have an empty struct ProblemSpecification for defining a uniform interface in the operators. Then specific operator classes can derive from struct ProblemSpecification and use static_cast<> to convert it to the actual specification object. Since struct ProblemSpecification is empty, it can go inside ibtk_utlities.h. That's how I am thinking about struct ProblemSpecification and its usage.

[boyceg]

Even though it is empty, I would still be inclined to put it in its own header.

[amneetb]

Hmm, sounds good. It should be in the IBTK namespace though?

[boyceg]

I think so.

[abarret]

I agree about putting it in it's own header. While now it's only being used for specifications for GeneralOperator, this could eventually be used elsewhere in the library.

[amneetb]

I added and made use of these objects in d8faa3bffe7302f9e0aca514dbfd7a6c654833c5