[MITgcm-devel] implicit formulation for bottom (and side) drag?
Jean-Michel Campin
jmc at mit.edu
Fri Oct 28 12:01:56 EDT 2016
Hi Dan,
I've started to code few things there, and was going to start some basic
tests. If you feel like trying this very soon, I could check-in
these pieces of code. Otherwise I will wait until my basic tests pass
before putting this into CVS.
Cheers,
Jean-Michel
On Wed, Oct 19, 2016 at 12:03:12PM -0400, Jean-Michel Campin wrote:
> Hi Dan,
>
> I've asked around and seems that having an implicit treatment of
> bottom drag (with implicit vertical viscosity, i.e., solution 3) would be
> helpful also for other hydrostatic applications.
> So I am going to see practically what this involves and will
> let you know about any progress here.
>
> Cheers,
> Jean-Michel
>
> On Mon, Sep 26, 2016 at 04:28:44PM -0400, Jean-Michel Campin wrote:
> > 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
> >
> > _______________________________________________
> > 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
More information about the MITgcm-devel
mailing list