Symmetry breaks quickly in simple flow

[mcgratta]

Run this case below. It only takes a few seconds. The resulting flow should be symmetric and isothermal, but it is neither. It is worst when CONSTANT_SPECIFIC_HEAT_RATIO=T, but it still fails if is F. If I remove one of the VENTs injecting the gas, the case remains iso-thermal, but obviously not symmetric.

&HEAD CHID='sym_test' /

&MESH IJK=64,1,16, XB=0.0,3.2,-0.001,0.001,0.0,0.8 /

&TIME T_END=8.0 /

&PRES CHECK_POISSON=T /

&MISC STRATIFICATION=F, NOISE=F, CONSTANT_SPECIFIC_HEAT_RATIO=T /
&RADI RADIATION=F /

&SPEC ID='GAS BG', MW=29, BACKGROUND=T /
&SPEC ID='GAS A', MW=29 /
&SPEC ID='GAS B', MW=29 /

&SURF ID='WALL', FREE_SLIP=T, DEFAULT=.TRUE. /

&SURF ID='BLOW A', VEL=-0.5, MASS_FRACTION(1)=1., SPEC_ID(1)='GAS A' /
&SURF ID='BLOW B', VEL=-0.5, MASS_FRACTION(1)=1., SPEC_ID(1)='GAS B' /

&VENT XB=0.0,0.0,-0.001,0.001,0.0,0.4, SURF_ID='BLOW A' /
&VENT XB=0.0,0.0,-0.001,0.001,0.4,0.8, SURF_ID='BLOW B' /
&VENT PBX=3.2, SURF_ID='OPEN' /

&SLCF PBY=0.,QUANTITY='MASS FRACTION', SPEC_ID='GAS A',  CELL_CENTERED=T, VECTOR=T /
&SLCF PBY=0.,QUANTITY='MASS FRACTION', SPEC_ID='GAS B',  CELL_CENTERED=T /
&SLCF PBY=0.,QUANTITY='MASS FRACTION', SPEC_ID='GAS BG', CELL_CENTERED=T /
&SLCF PBY=0.,QUANTITY='TEMPERATURE', CELL_CENTERED=T /
&SLCF PBY=0.,QUANTITY='DIVERGENCE',  CELL_CENTERED=T /

&DEVC ID='T_max', XB=0.0,3.2,-0.001,0.001,0.0,0.8, QUANTITY='TEMPERATURE', SPATIAL_STATISTIC='MAX', TEMPORAL_STATISTIC='INSTANT VALUE' /
&DEVC ID='T_min', XB=0.0,3.2,-0.001,0.001,0.0,0.8, QUANTITY='TEMPERATURE', SPATIAL_STATISTIC='MIN', TEMPORAL_STATISTIC='INSTANT VALUE' /

&DUMP DT_DEVC=0.00001 /

&TAIL /

[drjfloyd]

With a quick look. The .out file show these all have the same properties and if I add slices for specific heat, conductivity, etc. Those show identical values until temperature starts to change in which case conductivity increases with temperature. This doesn't seem to be something wrong with the functions for getting properties.

[mcgratta]

I don't think so either. I'm just surprised at how quickly the symmetry and iso-thermal break down. Like in 3 time steps. No time step changes, no changes in geom, no obsts, no mesh boundaries.

[drjfloyd]

It also doesn't happen if one of the VENTs is the background species or both vents are the same species. Only when the two vents are A and B

[mcgratta]

Nice clue.

[mcgratta]

How did you apply the BG species at the VENT. I get

ERROR(330): SURF BLOW B cannot use background species for MASS_FRACTION. (CHID: sym_test)

[mcgratta]

Nevermind. I know how.

[mcgratta]

Setting GVEC=0,0,0 maintains symmetry even though I would think we would not have to set that to do so.

[mcgratta]

When I run the case longer, the initial slop goes away. There seems to be something about the ramp up or initial/boundary conditions.

[drjfloyd]

It is still seen with just a 2 cell wide domain but not as strong

&HEAD CHID='test' /

&MESH IJK=8,1,4, XB=0.0,0.8,-0.001,0.001,0.2,0.6 /

&TIME T_END=2,DT=0.01 /

&PRES CHECK_POISSON=T /

&MISC STRATIFICATION=F, NOISE=F, CONSTANT_SPECIFIC_HEAT_RATIO=T /
&RADI RADIATION=F /

&SPEC ID='GAS BG', MW=29, BACKGROUND=T /
&SPEC ID='GAS A', MW=29 /
&SPEC ID='GAS B', MW=29 /

&SURF ID='WALL', FREE_SLIP=T, DEFAULT=.TRUE. /

&SURF ID='BLOW A', VEL=-0.5, MASS_FRACTION(1)=1., SPEC_ID(1)='GAS A' /
&SURF ID='BLOW B', VEL=-0.5, MASS_FRACTION(1)=1., SPEC_ID(1)='GAS B' /

&VENT XB=0.0,0.0,-0.001,0.001,0.3,0.4, SURF_ID='BLOW A' /
&VENT XB=0.0,0.0,-0.001,0.001,0.4,0.5, SURF_ID='BLOW B' /

&OBST XB=0,0.8,-0.001,0.001,0.2,0.3/
&OBST XB=0,0.8,-0.001,0.001,0.5,0.6/

&VENT PBX=0.8, SURF_ID='OPEN' /

&SLCF PBY=0.,QUANTITY='MASS FRACTION', SPEC_ID='GAS A', CELL_CENTERED=T, VECTOR=T /
&SLCF PBY=0.,QUANTITY='MASS FRACTION', SPEC_ID='GAS B', CELL_CENTERED=T /
&SLCF PBY=0.,QUANTITY='MASS FRACTION', SPEC_ID='GAS BG', CELL_CENTERED=T /
&SLCF PBY=0.,QUANTITY='TEMPERATURE', CELL_CENTERED=T /
&SLCF PBY=0.,QUANTITY='DIVERGENCE', CELL_CENTERED=T /
&SLCF PBY=0.,QUANTITY='CONDUCTIVITY', CELL_CENTERED=T /
&SLCF PBY=0.,QUANTITY='SPECIFIC HEAT', CELL_CENTERED=T /
&SLCF PBY=0.,QUANTITY='DENSITY', CELL_CENTERED=T /
&SLCF PBY=0.,QUANTITY='VISCOSITY', CELL_CENTERED=T /
&SLCF PBY=0.,QUANTITY='SPECIFIC ENTHALPY', CELL_CENTERED=T /

&TAIL/

[mcgratta]

If you set FLUX_LIMITER='GODUNOV' on the MISC line, the temperatures will be exactly 20 C. I worked my case above every which way, and it seems to me that the flux limiters other than CENTRAL and GODUNOV all create a small error. CENTRAL will eventually go "unstable" in the sense that tiny errors will grow and throw off the isothermal flow. I assume that the two simple limiters work because they act linearly when one sums the species equations, but I haven't nailed down a slam dunk reason.

[mcgratta]

On a related note, I just added a case called Pressure_Solver/obst_activation_default_gases.fds which is what this issue is based on. In this case, there are multiple meshes and obstructions that are crated and removed at mesh boundaries. I set the flux limiter to GODUNOV because, for the moment, I am trying to fix a small discrepancy when obstructions are created or removed at mesh boundaries in the presence of mixing gases. I think that this case be dealt with separately form the flux limiter issue.

[rmcdermo]

OK, I should be able to track this down. The other limiters should work too, so something is amiss.

[mcgratta]

In the case where we have three gases, two on the inlet and one background, and all have the same MW, when I sum the limited fluxes, I do not get what I expect with SUPERBEE, i.e. $\sum_\alpha \overline{\rho Z}_\alpha^{\rm FL} = \rho$, but with GODUNOV I do. I have not worked out the algebra to demonstrate that I should or should not expect this.

[rmcdermo]

I have a guess at what is going on. I sent some notes on Overleaf (they were not rendering correctly here on GitHub).

[rmcdermo]

I added a correction in this PR #13789 . For now, it is hidden behind the MISC flag TEST_FLUX_LIMITER_FACE_CORRECTION. The following input file now gives perfectly isothermal flow.

sym_test.fds.txt

Once we are happy with the coding, we can remove the flag. I think we just need to invoke this if N_TRACKED_SPECIES>2.

[mcgratta]

I'll run this locally. I'd rather not commit a potential grenade because it causes problems when doing a git bisect.

[mcgratta]

I am running the verification cases to see if something will blow up. I set the new parameter to T in the test. But I just ran the sym_test case that is posted above, and I gave Gas A a MW of 25 and Gas B a MW of 35. The temperature does not remain isothermal.

We can work on this more tomorrow. I think we need to run the idea through a series of increasingly complex cases -- multi-mesh, different gases, flux limiters, obstructions, etc.

[rmcdermo]

I think I know what is going on. I updated the notes. To maintain isothermal flow, we need to keep $T = \frac{\overline{W}p}{\rho R}$ constant at all cells and cell faces. The current correction does not consider this. So, I think I need to add a molecular weight factor such that $\sum_\alpha \frac{(\rho Z_\alpha)}{W_\alpha} = \frac{\rho}{\overline{W}}$ is maintained at the face.

[mcgratta]

PR #13789 with the test flag set to T knocked about a dozen cases out of tolerance. Since there is no immediate crisis, I suggest we revert the change and work it over the next few weeks. That will give us time to progress from the simple cases to the complex, like couch.fds which was one of the cases that failed. In the meantime, I want to figure out what is happening at mesh boundaries when obstructions are created or removed This can be done using the simple GODUNOV limiter because I think the issue there is separate. It would be easier not to have the test code complicating things until it's ready.

[rmcdermo]

Please do not revert. I will deal with this. I am working on the correction. Please just leave the flag to false.

…

Sent from my iPhone

On Fri, Nov 22, 2024 at 8:06 AM Kevin McGrattan ***@***.***> wrote: PR #13789 <#13789> with the test flag set to T knocked about a dozen cases out of tolerance. Since there is no immediate crisis, I suggest we revert the change and work it over the next few weeks. That will give us time to progress from the simple cases to the complex, like couch.fds which was one of the cases that failed. In the meantime, I want to figure out what is happening at mesh boundaries when obstructions are created or removed This can be done using the simple GODUNOV limiter because I think the issue there is separate. It would be easier not to have the test code complicating things until it's ready. — Reply to this email directly, view it on GitHub <#13764 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/ABBWXQOXWF6H6OBJOFXNO6T2B5B7PAVCNFSM6AAAAABR33IOEGVHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMZDIOJTHE3DQNZRGM> . You are receiving this because you commented.Message ID: ***@***.***>

[mcgratta]

I transferred this issue to Randy. The best test case for this is in /Verification/Pressure_Solver/obst_activation_default_gases.fds. You would just need to remove the flag for GODUNOV.

[rmcdermo]

I made a bit of progress on this issue today. This PR #13854 is working for this test case:

isothermal_3spec_const_gamma.fds.txt

However, this case is single mesh. When I go to multi-mesh I end up with errors at the INTERPOLATED_BC which is what I'm working on now.

[drjfloyd]

Here is a test case. The attached case defines five species that allow for testing with the same MW and CP or varying MW and CP. Eight meshes have a two cell high gas phase domain contracting down to one cell high. Pairs of species are injected on the left at the same or varied temperature. For varying temperature comments in the input file note what the solution should be.

test.txt

[rmcdermo]

This PR #13863 fixes the temperature issue at interpolated boundaries for CONSTANT_SPECIFIC_HEAT_RATIO=T (still). I will now start working the variable CP aspects using Jason's case.

Another note is that these cases still require CFL_VELOCITY_NORM=1 and CFL_MAX=0.8 to stay perfectly isothermal. Here is the working input file for obst_activation_default_gases.

test.fds.txt

[rmcdermo]

@drjfloyd Have your case working in a branch. Need to run it through firebot then I'll push it up.

[rmcdermo]

This should be solved now with FLUX_LIMITER_MW_CORRECTION=T, which is now default for LES, PR #13891

I have pushed this also the verification case here PR #13893 .

I'll close this, but we can reopen if we see issues. I think that improving efficiency is something that needs its own thread.

[drjfloyd]

I just ran firebot with the current source with and without the flux correction set to true. The symmetry_test_2 case in flowfields failed with it true. The plot shows the difference in left and right temperature. With it true the difference is almost a degree vs. .03 degree with it false.

[rmcdermo]

OK, thanks. I'll take a look.

[rmcdermo]

@drjfloyd I'm not sure this is something we can avoid. If I set FLUX_LIMITER_MW_CORRECTION=F and run the case to 8 seconds, then I get this,

So, this is what we get by default.

If I set,

&MISC NOISE=.FALSE., FLUX_LIMITER_MW_CORRECTION=T, CFL_VELOCITY_NORM=1 /

then I get this,

I attribute the earlier onset of the randomness simply to the fact that the correction requires more flops and so the roundoff error accumulates a bit faster.

What I can do is to automatically set CFL_VELOCITY_NORM=1 if FLUX_LIMITER_MW_CORRECTION=T.

I don't think it makes sense to try to improve this at the moment, as we have another idea we want to pursue that would both eliminate the need for a special correction and hopefully also speed up the whole transport algorithm. But let me know if you disagree.

[drjfloyd]

That all makes sense. I agree.