About modifying the physical property calculation method for calculating LMN.

[fanxiaojie69-dotcom]

Dear incompact3d developers,
I'd like to ask about a problem you're encountering with Incompact3D, and about modifying the physical property calculation method for calculating LMN.
In the original Icompact3D code, I modified the subroutine for calculating temperature and density to call coolProp to generate a table of different temperatures at a specific pressure, then interpolate and calculate the relevant physical properties. I then made modifications to the temperature and density calls, as well as the stretching code you mentioned on the forum. However, after adding the stretching code, I still experience discontinuous density display in ParaView.
Below is a screenshot of my modified code and a ParaView visualization. I'd appreciate your response. Best wishes.

Code to calculate temperature

Code to calculate density

Code for transeq.f90


Code for time_integrators.f90



Code for Navier.f90

paraview

[pbartholomew08]

Hello, the temperature transport configuration is not as used/debugged as the density transport formulation, can you convert your code to try this?

What boundary conditions are you using at the top/bottom for temperature/density?

Finally, depending on how you are redefining your EOS you might need to make a more fundamental change to ensure the velocity divergence is correctly formulated - see https://arc.aiaa.org/doi/10.2514/3.9371

[fanxiaojie69-dotcom]

This is how I set it up in the input.i3d file：
nclx1 = 0
nclxn = 0
ncly1 = 2
nclyn = 2
nclz1 = 0
nclzn = 0 and
&ScalarParam
!================

Sc(1) = 1 ! Schmidt number

nclxS1 = 0
nclxSn = 0
nclyS1 = 2
nclySn = 2
nclzS1 = 0
nclzSn = 0

/End
This is how I set the upper and lower bounds in case-channel.f90
phi1(:,:,:,:) = 0.49 !change as much as you want
if ((nclyS1 == 2).and.(xstart(2) == 1)) then
!! Generate a hot patch on bottom boundary
phi1(:,1,:,:) = 0
endif
and
if (iscalar /= 0) then
if (iimplicit <= 0) then
if ((nclyS1 == 2).and.(xstart(2) == 1)) then
!! Generate a hot patch on bottom boundary
phi(:,1,:,:) = 0
endif
if ((nclySn == 2).and.(xend(2) == ny)) THEN
phi(:,xsize(2),:,:) = 1
endif
else

[pbartholomew08]

Actually, to try and narrow down the issue - can you try using your equation of state in the mixing layer example without stretching? In non-dimensional terms the "hot" flow temperature is 2x the "cold" flow.

[fanxiaojie69-dotcom]

OK. Thank you for your answer. I'll try to narrow down the scope of the question.