[MITgcm-devel] implicit formulation for bottom (and side) drag?

[Jean-Michel Campin]

Hi Dan and others,

Would like to clarify a coupled of things regarding implicit bottom drag
(and back to the list after few other exchanged msg):

1) Non-hydrostatic (using 3-D solver) is more complicated, and even with current
 implementation of implicit vertical viscosity (only applied to horiz. Mom), 
 it's probably not right (3-D solver matrix does not commute with vertical 
 viscosity operator, even if it's only applied to U,V and not to W). 
 Bottom drag is adding an other layer of difficulties.

2) It would be good to have a way to keep the explicit bottom friction stable
 (to prevent the model to blow-up) by limiting the magnitude of it, as you 
 proposed here - just that I did not like the names ...

3) Hydrostatic only: it turns out that it's very easy to treat implicitly the 
 bottom drag (and same for side drag), just involves to change the CG2D matrix
 the right way and the modified matrix remains symmetric (good !).

 The problem here is that this simple solution cannot be applied with implicit
 vertical viscosity. There should be a way to write it down with both implicit 
 viscosity and implicit bottom drag but it's no longer "very simple": need to 
 keep track of the LU decomposition of vertical viscosity operator and add it 
 inside the CG2D matrix. Still trying to figure out if some terms will cancel, 
 but it's not obvious that they should.
 
4) Otherwise, we could try to use an approximate method. 
 There are many options there but the simplest would be to define & use an 
 intermediate stage (u*): 
 a) solve vertical-viscosity (+ bottom drag) implicitly but using -Grad.(eta^n):
  u* = u^n + Delta_t { Gu + d/dz ( nu d/dz u* ) - Epsil_b Cd u* - g Grad.(eta^n) }
  where Epsil_b = 1 or 0 if at bottom level or not;
 b) and then solve (CG2D) for just the adjustment (increment) of eta, i.e.:
     eta^n+1 - eta^n
  including (or not) an implicit bottom drag contribution: - Cd ( u^n+1 - u* )  

 And just to clarity: currently, when (a) is done, the SSH gradient is missing,
 and since it's one of the larger term (if not the largest), it would be
 really wrong to include an implicit bottom drag term in our current implicit
 vertical viscosity solver.

 I have been thinking of this implementation mainly for the NH case (and had 
 this in mind in my previous answer) for which a truly (non-approximate) 
 implicit method is going to be complicated. 
 The advantage is that this method does not require to get the LU matrix inside 
 the 2-D (or 3-D solver) --> minor modifications. 
 But need to solve (in CG2D and CG3D) for the increment of eta / P_nh instead of
 the full field eta^n+1 / P_nh^n+1. This does not look bad but still significant
 changes (RHS normalisation ? CG convergence criteria ? ... ) and also to update
 how filter (FFT, Shapiro) are applied.

Cheers,
Jean-Michel

On Tue, Aug 30, 2016 at 11:22:41AM -0400, Jean-Michel Campin wrote:
> Hi Dan and Martin,
> 
> I have a problem regarding naming: 
> 1) I think that if one term is explicit, it shoud not be called implicit. 
>   This is a source of confusion that will even get worse when 
>   the real implicit treatment of bottom drag is implemented (not allowed with current
>   solve_for_pressure but plans are to change it to allow to try such implicit treatment).
> 2) The other thing is that I currently don't see a modified side-drag formulation in
>  this contrib dir: http://mitgcm.org/viewvc/MITgcm/MITgcm_contrib/dgoldberg/impl_drag/
>  and since we have a clear distinction between bottom drag and side drag, I think it 
>  would be better to keep it for the modified bottom drag term you propose, with a 
>  name that does not contain "SIDEDRAG" but only BOTTOMDRAG.
> 
> Cheers,
> Jean-Michel
> 
> On Tue, Aug 30, 2016 at 04:49:22PM +0200, Martin Losch wrote:
> > Hi Dan,
> > 
> > So far the drag is fully explicit (AFAIK).
> > I like your suggestion. I would actually make this ???more invasive??? (if it does not change results): remove all 
> > 
> > #ifndef IMPLICIT_BOTTOMSIDEDRAG
> >      &      * uFld(i,j)
> > #endif
> > 
> > and have an extra loop
> > 
> > #ifdef ALLOW_IMPLICIT_BOTTOMSIDEDRAG
> >      IF ( useImplictBottomSideDrag ) THEN
> >       DO j=1-OLy+1,sNy+OLy-1
> >        DO i=1-OLx,sNx+OLx-1
> >         uDragTerms(i,j) = uDragTerms(i,j)*uFld(i,j) / 
> >      &        (1. - deltaTmom*uDragTerms(i,j))
> >        ENDDO
> >       ENDDO
> >      ELSE
> > #else
> >      IF (.TRUE.)
> > #endif
> >       DO j=1-OLy+1,sNy+OLy-1
> >        DO i=1-OLx,sNx+OLx-1
> >         uDragTerms(i,j) = uDragTerms(i,j)*uFld(i,j)
> >        ENDDO
> >       ENDDO
> >     ENDIF
> > 
> > It would be more in-line with the general procedure (1) enable code with CPP flag, (2) turn it on with runtime flag. (1) is necessary to hide it from TAF/OpenAD if it create recomputations (division by active variable).
> > 
> > Martin
> > 
> > 
> > > On 18 Jul 2016, at 11:32, Daniel Goldberg <dngoldberg at gmail.com> wrote:
> > > 
> > > Hi All,
> > > 
> > > I am experimenting with some very small hFacMin values and experiencing instabilities as a result. 
> > > 
> > > Is it correct that there is no implicit formulation for bottom drag? If this is true, please see http://mitgcm.org/viewvc/MITgcm/MITgcm_contrib/dgoldberg/impl_drag/
> > > 
> > > The change is not invasive (and hidden with a CPP directive.) The drag formulations take the form
> > > gU -= F(u,v) * u
> > > and i replace by
> > > gU -=  F(u,v) / (1 + deltaTmom * F(u,v) ) * u
> > > which is not quite implicit but will help to stabilize when drag terms are high (e.g. when hFacW is small) and will hopefully not affect things too much when they are small.
> > > 
> > > If not I wonder if anyone has any thoughts/objections on this modification to allow for such a thing? In addition, it seems something similar can be done for side drag if biharmonic horizontal viscosity is not used?
> > > 
> > > Many thanks
> > > Dan
> > > 
> > > 
> > > 
> > > 
> > > -- 
> > > 
> > > Daniel Goldberg, PhD
> > > Lecturer in Glaciology
> > > School of Geosciences, University of Edinburgh
> > > Geography Building, Drummond Street, Edinburgh EH8 9XP
> > > 
> > > 
> > > em: Dan.Goldberg at ed.ac.uk
> > > web: http://ocean.mit.edu/~dgoldberg
> > > _______________________________________________
> > > MITgcm-devel mailing list
> > > MITgcm-devel at mitgcm.org
> > > http://mitgcm.org/mailman/listinfo/mitgcm-devel
> > 
> > 
> > _______________________________________________
> > MITgcm-devel mailing list
> > MITgcm-devel at mitgcm.org
> > http://mitgcm.org/mailman/listinfo/mitgcm-devel
> 
> _______________________________________________
> MITgcm-devel mailing list
> MITgcm-devel at mitgcm.org
> http://mitgcm.org/mailman/listinfo/mitgcm-devel