Suggestions for further SGS improvements:

* filter radius issues: replace delta_target; include dx,dz 
* decide on VD damping or delta shape for bottom boundary. delta works OK,
  but parameter c_mild may be tunable, as specified.
* check behavior of routines in 3d (not tested yet)
* Inclusion of Ri-dependent depedent behavior in the vertical direction
* Inclusion of Ri-dependent depedent behavior in turbulent Prandtl number 

Suggestion for scientific qeustion:

* What is the behavior of the bottom gravity current as a function of Ra 
(the only remaining physical parameter), how do turbulent structures, entrainment, 
salt flux change over  a wide range, e.g., 10^4 <= Ra <= 10^8 ?

Best parameters for homog. Smag. SGS (no Ri-dependence):

(a) 0.03 <= cs <= 0.05 (or "normal" Cs=0.17 to 0.22, consistent with published results)

(b) Also some hidden parameters for delta, to take care of the bottom boundary:
    c_mild=300 (or I may suggest even larger values, 500, 1000?)
    delta_target=0.01 (for this mesh, should be a dynamic function)
    no VD damping (replaced by Traian's delta distributions)

#####################################################################################


Saturday - 7 August 2004:

###################################################################################
 
ref exp:  Ra=5.e6; r=0.02; c_s=0

step (x1000) |   1   |    2   |   3   |   4   |   5   |   6   |  7   |
----------------------------------------------------------------------
Tmax         | 1.023 |  1.138 1.173 | 1.158 |  1.133 | 1.066 | 1.046 |

runs fine until the end of the domain, 12000 steps.

It is important that we find a way to reduce Tmax if we wish to attain
an advantage over constant r case. So, lets' focus on this issue...

****************************************************************************
exp smari39a:

Ra=5.e6
r=0.
c_s=10.,let's start with super high value to notice the difference  
factor=1, remove Ri-dependence to see how much vertical mixing takes place
          with homog. smag.
c_mild=20, Traian's default
delta_target=0.01, Traian's default, seems OK for deltay but small for deltax?
delta=delta_y, Traian's default, x direction not considered for the time being

step (x1000) |   1   |    2   |   3   |   4   |   5   |   6   |  7   |
----------------------------------------------------------------------
Tmax         | 1.093 |  1.307 | 1.309 | 1.299 |  1.253| 1.171 | 1.128 |

runs fine until the end of the domain, 12000 steps. 

Super smooth in general, but still some KH formation.

Overshoots are typically points under the nose. Given that Ri-dependence
is removed, it must be related to delta, (a) distribution of delta - cmild value
or (b) magnitude - delta_target.

****************************************************************************

exp smari39b:

same as smari39a, but
delta_target=0.02

step (x1000) |   1   |    2   |   3   
--------------------------------------
Tmax         | 1.085 |  1.285 | 1.329 

blew up at 3442 steps. 

Much more diffused than smari30a, but small overshoots remain at the bottom.
Model is very sensitive to delta_target!

****************************************************************************

exp smari39c:

same as smari39a, but
c_mild=100 (c_aggresive...)

step (x1000) |   1   |    2   |   3   |   4   
----------------------------------------------
Tmax         | 1.084 |  1.285 | 1.297 | 1.306 

It is somewhat smoother than smari39a, but generally the same.
The bottom overshoots remain elusive, they must be related
to Van Priest damping?

****************************************************************************

exp smari39d:

from smari - subroutine compute_sma:

       y_tilde = 0.06*(10.0-xm1(kx,ky,1,ie)) : 

                bottom level, which is 0.6 at x=0 and 0 at x=10, OK.

       VD_const = 1.0 - exp(ym1(kx,ky,1,ie)-y_tilde)

                which is 0 at the bottom (OK) and -0.5 at the top at x=0
                and -1.7 at the top at x=10... what is this? square of VD_const is
                used, so (-) is not a problem, but exp() is not consistent?  

Eq. (50) from Hughes et al. (2001) on Van Driest damping is:

vT=(cs*delta*(1-exp(-y'/A)))^2*|S|

such that for large y', vT=(cs*delta)^2*|S|.

This needs to be modified, change to:

       VD_const = 1.0 - exp(-(ym1(kx,ky,1,ie)-y_tilde)/Aplus) 

Now I see why Traian has Aplus floating around, with a default of
Aplus=25, which will not work here (needs to be <1).
Let's start with 
Aplus=1.

Essentially the same as smari39c, which is expected.

****************************************************************************

exp smari39e:

Why do we need both delta and Van Driest function to go to zero at the wall?
Let's keep it simple, just use one of them... Let's eliminate Van Priest damping,
and play with Traian's delta functions.

VD out, delta_target=0.01, c_mild=100: blew up at  40 steps
  "           "            c_mild= 20: blew up at 156 steps

Probably super large c_s is the problem now?

   "          "       "  , c_s=0.1   : 

OK, works

step (x1000) |   1   |    2   |   3   |   4   |   5   
------------------------------------------------------
Tmax         | 1.058 |  1.179 | 1.189 | 1.224 |  1.066

this case is very similar, in terms of snapshots, to ref exp
with r=0.02. It is only slightly more diffused. This is not-so-bad
news given that we do not manipulate vertical viscosity, and it is
homog. Smag model. However, bottom overshoots remain... 
Lets try  one more case with the aggressive profile, which I like, 
based on Traian's delta profles from his website document):

VD out, delta_target=0.01, c_s=0.1, c_mild=100

step (x1000) |   1   |    2   |   3   |   4   |   5   
----------------------------------------------------
Tmax         |   1   |  1.006 | 1.002 |   1   |   1

Snapshots are somewhat more diffused than ref case with r=0.02,
but no more overshoots; the model is quite healthy and
both cs and delta are in reasonable range, It seems like
a successful experiment with SGS...

Let's keep in mind that conventional range for
Smagorinsky constant is:
0.1  <=  CS  <= 0.24 which corresponds to our implementation
0.01 <=  c_s <= 0.06 so that c_s is still little high.

****************************************************************************

exp smari39f:

- Let's have a real test: 

Ra=20.e6 (4x increase wrt ref)
dt=0.5e-6 (however same number of steps necessary to travel trough)

crashed and burned at 2300 steps... Looked good until then though.
Again the bottom problem? Maybe sharpening delta profile further
will help:

- c_mild=300 (very aggressive):

OK, it works, and Tmax=1 at all times.
However, the results are nearly indentical to smari39e. It seems
that the increase in shear via higher buoyancy forcing produces 
enough diffusion to suppress emergence of further modes. 
Maybe c_s is too high, let's crank it
down to the middle of the conventional range:

- c_s=0.03 

Seems to work quite well, generating sharper features and more turbulence.

Let's push it further:
-Ra=200.e6 (100x increase wrt ref); dt=0.1e-6; c_s=0.05 (little increase
just to be conservative)

results are nearly identical to Ra=20.e6 case!

Let's check what ref exp would do:

- c_s=0; r=0.02, Ra=20.e6

it runs and generates more structure than Smag-SGS, but it is not healty,
Tmax>1.4, and increase in r would probably lead to something similar.


****************************************************************************

Let's experiment with Ri-stuff:

Activate Ri-dependent factor: crashed and 200 steps.

Zero vertical viscosity is not received well or what is happening?
More debugging necesary, i.e. is Ri calculated right? 

Modify as:

            if (j.eq.2) then
               factor = vert_diff_ratio
            else
               factor = 1.0
             endif

where vert_diff_ratio=0.1 
to try it. In this case, it is the ratio
of vertical and horizontal Smagorinsky viscosities; no Ri dependence...

-> sharper features are generated, routines seem to work OK.

Now modify as:

             Ri = -Ra*w_T(kx,ky,kz,1,2)  
! there seems to be a s ign mistake and abs not neccesary
     $          / (w(kx,ky,kz,1,2)**2 + w(kx,ky,kz,3,2)**2)
           ...

           if (j.eq.2) then
               factor = factor_Ri
            else
               factor = 1.0
            endif

The code hangs... something is wrong.. The code hates factor_Ri.
Nek500 is also 2-3x slower, U-Press requires much more iterations.

Let's resort to our old friend,
vert_diff_ratio=1.e-5 in subroutine reset_vert_diff
as something very small as "background" vertical viscosity.

No, doesn't work, blew up quite early.
Not much success with Ri-stuff, enough for today.

#############################################################################


Sunday - 8 August 2004:


#############################################################################


Yesterday was pretty good. Let's summarize the results of the day:

1) After taking care of a few things regarding the bottom boundary,
   homog. Smag. model is now in good working condition in nek5000,
   with all routines working well in 2d.

2) Hypothesis:
Initial results seem to indicate that homog. Smag. gives very similar results
to anisotropic Laplacian, with the added advantage of requiring no
justification of the famous parameter, r=vert_diff_ratio. 
Noting that we used as a minimum parameter:

r=0.02-0.05 for Ra=5e6  (2004 papers)
r=1.e-5     for Ra=1.e5 (2002 paper)
r=1         for Ra=5.e8 (Pearl experiments)

to keep the model stable while generating max. KH rolls, anisotropic
Laplacian requires adjustment of this parameter manually, depending on
forcing (Ra) and probably resolution as well. 
In contrast, homog. Smag SGS may provide an automatic way to adjust
dissipation depending on forcing (and probably resolution as well - needs to
be tested).

It would be surprisingly good news. If that is the case, we can attack immediately
the scientific question of, for instance:

* What is the behavior of the bottom gravity current as a function of Ra 
(the only remaining physical parameter), how do turbulent structures, entrainment, 
salt flux change over  a wide range, e.g., 10^4 <= Ra <= 10^8 ?

So, let's focus on this issue today.

********************************************************************************

* smari40a: Ra=5.e6, dt=1.e-6, r=0
           cs=0.05, cmild=300, delta_target=0.01, no VD, no factor_Ri 

(let's change the name, finally, and rerun ref exp with sgs to start).

step (x1000) |   1   |    2   |   3   |   4   |   5   |   6  |   7  
----------------------------------------------------------------------
Tmax         |   1   |  1.005 | 1.001 |   1   |   1   !   1  !   1

Nearly identical results to ref case with r=0.02; only sligthly
more diffused as cs=0.05 is conservatively high (very well behaved experiment)
and r=0.02 is little to brave (leading to high Tmax, see above).

********************************************************************************

* smari40b: Ra=312500 (16x smaller wrt smari40a), dt=4.e-6 (4x larger)... cs=0.05 ...

Let's crank down the forcing. 
Can it still generate KH by reducing appropriately the dissipation ? 

-> Formation of head, but not much else, no KH, probably because cs is too high
but it is supposed to be a constant. Maybe there is still some room for 
Ri-dependent coefficients, but one must expect less turbulence anyway.

********************************************************************************

* smari40c: Ra=80.e6 (16x larger wrt smari40a), dt=0.25e-6 (4x smaller); cs=0.05 

Let's go the other way, crank up the forcing. We expect more turbulence,
which is of course much more testing on the sgs model:
will SGS hold up without touching any parameter?

step (x1000) |   1   |    2   |   3   |   4   |   5   |   6  |   7  
----------------------------------------------------------------------
Tmax         | 1.017 |  1.045 | 1.061 | 1.042 |   1   !   1  !   1

Runs fine (hard stop at 12000 steps).

This is very impressive, just crank up Ra, and SGS supplies
necessary dissipation where needed. No wonder Smag. is so popular!

Dynamically, nice KH rolls appear to be replaced by a bunch of
hydraulic jumps, with rolls detaching off their tips. Interesting because
hydraulic jumps circulate the opposite direction (clockwise) 
as KH rolls (ccw). Also, "head" disintegrates after 9000 steps.

********************************************************************************

* smari40d: Ra=320.e6 (64x larger wrt smari40a), dt=0.125e-6 (8x smaller), cs=0.05

Push further...

step (x1000) |   1   |    2   |   3   |   4   |   5   |   6  |   7  
----------------------------------------------------------------------
Tmax         | 1.026 |  1.115 | 1.055 | 1.078 | 1.040 !   !     !   

It runs (stpped manually at 6100 steps, no problems)
This is impressive: SGS model automatically adjust
dissipation via strain tensor, and keeps the model from blowing up.
HOWEVER, there is virtually NO increase in the level of turbulence
wrt case with Ra=80.e6 (snapshots are identical), 
hence I now understand the argument that Smag. model is too dissipative.
It also does not employ any information about temperature, only velocity
strain. Maybe there is room for Ri-dependent formulation.