Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open boundary questions #293

Closed
kshedstrom opened this issue May 11, 2016 · 16 comments
Closed

Open boundary questions #293

kshedstrom opened this issue May 11, 2016 · 16 comments

Comments

@kshedstrom
Copy link
Collaborator

We are agreed that the open boundary points lie within the computational domain, aren't we? There is evidence that this was not the plan for the "Flather" points in the past, but I can clean it out as I find it.

If the open boundaries are inside the domain, I don't know why symmetric mode would be needed.

Are we changing the Flather name?

@adcroft
Copy link
Collaborator

adcroft commented May 13, 2016

I think the OB points can legitimately live on the edge of the [physical]
domain but only when the model is in symmetric mode because the
western/southern boundaries do not "belong" on any processor in
non-symmetric mode.

We've had a long standing question about why symmetric is required for OBCs
and I think it not related to the above point but more to do with the
timing of when things are halo-updated. Halo-updates in non-symmetric mode
would clobber any values imposed by a radiation condition on a
western/southern edge (referring to potentially interior tiles in a
parallel decomposition). In symmetric mode, the halo-updates do not touch
any edge values of the computational domain so we can specify values (via a
radiation routine) and then rely on the values being retained. It seems to
me we could find every halo-update and re-apply the radiation values but it
does seem like more work. I've bcc'd Mehmet Ilicak who once figured all
this out and explained it to us (convincingly then) even though I'm not
sure I am explaining it correctly now. Hopefully he'll chime in.

Flather is meant to refer to the barotropic radiation condition, yes? I
think the renaming idea was to cleanly separate the barotropic radiation
code from other [drafts] of radiation conditions which might have borrowed
the parameters/flags/code. My memory is poor so I might be confused about
what we discussed.

-A.

Dr Alistair Adcroft (Alistair.Adcroft@noaa.gov)
Princeton University Tel: (609) 987-5073
NOAA/GFDL, 201 Forrestal Road, Princeton, NJ 08540

On Wed, May 11, 2016 at 6:50 PM, Kate Hedstrom notifications@github.com
wrote:

We are agreed that the open boundary points lie within the computational
domain, aren't we? There is evidence that this was not the plan for the
"Flather" points in the past, but I can clean it out as I find it.

If the open boundaries are inside the domain, I don't know why symmetric
mode would be needed.

Are we changing the Flather name?


You are receiving this because you are subscribed to this thread.
Reply to this email directly or view it on GitHub
NOAA-GFDL#293

@kshedstrom
Copy link
Collaborator Author

OK, I've got something that runs (briefly) with "Flather" in CCS1. I had to
hack at things a bit and it changes the results of the circle_obc case. I
have to see what that case is all about...

I have no problem with staying in symmetric land. Whatever.

Yes, Flather wrote a paper about barotropic OBCs. In MOM6, the term is used
for any OBC - the 2D and 3D parts are all intertwined.

Kate

On Thu, May 12, 2016 at 4:35 PM, Alistair Adcroft notifications@github.com
wrote:

I think the OB points can legitimately live on the edge of the [physical]
domain but only when the model is in symmetric mode because the
western/southern boundaries do not "belong" on any processor in
non-symmetric mode.

We've had a long standing question about why symmetric is required for OBCs
and I think it not related to the above point but more to do with the
timing of when things are halo-updated. Halo-updates in non-symmetric mode
would clobber any values imposed by a radiation condition on a
western/southern edge (referring to potentially interior tiles in a
parallel decomposition). In symmetric mode, the halo-updates do not touch
any edge values of the computational domain so we can specify values (via a
radiation routine) and then rely on the values being retained. It seems to
me we could find every halo-update and re-apply the radiation values but it
does seem like more work. I've bcc'd Mehmet Ilicak who once figured all
this out and explained it to us (convincingly then) even though I'm not
sure I am explaining it correctly now. Hopefully he'll chime in.

Flather is meant to refer to the barotropic radiation condition, yes? I
think the renaming idea was to cleanly separate the barotropic radiation
code from other [drafts] of radiation conditions which might have borrowed
the parameters/flags/code. My memory is poor so I might be confused about
what we discussed.

-A.

Dr Alistair Adcroft (Alistair.Adcroft@noaa.gov)
Princeton University Tel: (609) 987-5073
NOAA/GFDL, 201 Forrestal Road, Princeton, NJ 08540

On Wed, May 11, 2016 at 6:50 PM, Kate Hedstrom notifications@github.com
wrote:

We are agreed that the open boundary points lie within the computational
domain, aren't we? There is evidence that this was not the plan for the
"Flather" points in the past, but I can clean it out as I find it.

If the open boundaries are inside the domain, I don't know why symmetric
mode would be needed.

Are we changing the Flather name?


You are receiving this because you are subscribed to this thread.
Reply to this email directly or view it on GitHub
NOAA-GFDL#293


You are receiving this because you authored the thread.
Reply to this email directly or view it on GitHub
NOAA-GFDL#293 (comment)

@kshedstrom
Copy link
Collaborator Author

The blow-up is happening in the barotropic mode, in the SE corner next to the land. Things go from ugly to unstable quite suddenly. Might be time for a simpler problem than CCS1.

@kshedstrom
Copy link
Collaborator Author

With the inside OBCs, circle_obcs blows up after ten hours in the two southern corners. This can be fixed by using OBC_SIMPLE (u, v = 0) for the eastern and western v boundaries, also northern and southern u boundaries.

@adcroft
Copy link
Collaborator

adcroft commented May 25, 2016

Hi,

Yes, Flather is a specific condition of Orlanski type radiation condition for specifically baratropic case since in the baratropic model the wave speed is sqrt(gh).
I would prefer to keep that name. Because other models such as ROMS is using the same nomenclature. It would be helpful for people if they come from another model environment.
Simply baroclinic open boundary conditions are Orlanski type radiation and Barotropic is Flather.

I think the update was the issue in non-symmetric mode for open BC. That’s why it did;t work with that.
The open BC has to be on the physical domain. If I remember correctly the wave will bounce back if you change that.

Best,
Mehmet

On May 13, 2016, at 2:34 AM, Alistair Adcroft - NOAA Affiliate alistair.adcroft@noaa.gov wrote:

I think the OB points can legitimately live on the edge of the [physical] domain but only when the model is in symmetric mode because the western/southern boundaries do not "belong" on any processor in non-symmetric mode.

We've had a long standing question about why symmetric is required for OBCs and I think it not related to the above point but more to do with the timing of when things are halo-updated. Halo-updates in non-symmetric mode would clobber any values imposed by a radiation condition on a western/southern edge (referring to potentially interior tiles in a parallel decomposition). In symmetric mode, the halo-updates do not touch any edge values of the computational domain so we can specify values (via a radiation routine) and then rely on the values being retained. It seems to me we could find every halo-update and re-apply the radiation values but it does seem like more work. I've bcc'd Mehmet Ilicak who once figured all this out and explained it to us (convincingly then) even though I'm not sure I am explaining it correctly now. Hopefully he'll chime in.

Flather is meant to refer to the barotropic radiation condition, yes? I think the renaming idea was to cleanly separate the barotropic radiation code from other [drafts] of radiation conditions which might have borrowed the parameters/flags/code. My memory is poor so I might be confused about what we discussed.

-A.

Dr Alistair Adcroft (Alistair.Adcroft@noaa.gov mailto:Alistair.Adcroft@noaa.gov)
Princeton University Tel: (609) 987-5073
NOAA/GFDL, 201 Forrestal Road, Princeton, NJ 08540

On Wed, May 11, 2016 at 6:50 PM, Kate Hedstrom <notifications@github.com mailto:notifications@github.com> wrote:
We are agreed that the open boundary points lie within the computational domain, aren't we? There is evidence that this was not the plan for the "Flather" points in the past, but I can clean it out as I find it.

If the open boundaries are inside the domain, I don't know why symmetric mode would be needed.

Are we changing the Flather name?


You are receiving this because you are subscribed to this thread.
Reply to this email directly or view it on GitHub NOAA-GFDL#293

@kshedstrom
Copy link
Collaborator Author

To summarize from the latest hangout and other thoughts on open boundaries:

  • The model boundary is some line through internal q points.
  • h points just outside the line are also part of the boundary conditions.
  • I'd like to use the land mask to separate wall parts of the boundary from
    open parts of the boundary. The open parts need valid h values out there,
    which isn't the case inside the land mask anyway.
  • These external h points will not be involved in the regular timestepping
  • perhaps a boundary mask can be used to shut that off. I got circle_obcs
    to run by using the OBC_SIMPLE with u,v=0 to shut off continuity out there
  • not exactly the general solution we want.
  • If the open boundary is more than one grid from the edge, points beyond
    the line of open boundary points can be marked as land.
  • I plan to have a setup which doesn't need messing with the user directory
    where the open boundaries are one grid in from the edge.

For the naming, I'm sure we'll have to go beyond simple Orlanski. The ROMS
world calls the barotropic "Chapman-Flather", the baroclinic
"radiation-nudging" (or Marchesiello). Sometimes Shchepetkin is mentioned
too. ROMS allows each field and each side to have various settings - T and
S are clamped east, radiation north, with nudging south, while dye sees a
wall all around. I don't know how general we need to be.

Kate

On Wed, May 25, 2016 at 4:22 AM, Alistair Adcroft notifications@github.com
wrote:

Hi,

Yes, Flather is a specific condition of Orlanski type radiation condition
for specifically baratropic case since in the baratropic model the wave
speed is sqrt(gh).
I would prefer to keep that name. Because other models such as ROMS is
using the same nomenclature. It would be helpful for people if they come
from another model environment.
Simply baroclinic open boundary conditions are Orlanski type radiation and
Barotropic is Flather.

I think the update was the issue in non-symmetric mode for open BC. That’s
why it did;t work with that.
The open BC has to be on the physical domain. If I remember correctly the
wave will bounce back if you change that.

Best,
Mehmet

On May 13, 2016, at 2:34 AM, Alistair Adcroft - NOAA Affiliate <
alistair.adcroft@noaa.gov> wrote:

I think the OB points can legitimately live on the edge of the
[physical] domain but only when the model is in symmetric mode because the
western/southern boundaries do not "belong" on any processor in
non-symmetric mode.

We've had a long standing question about why symmetric is required for
OBCs and I think it not related to the above point but more to do with the
timing of when things are halo-updated. Halo-updates in non-symmetric mode
would clobber any values imposed by a radiation condition on a
western/southern edge (referring to potentially interior tiles in a
parallel decomposition). In symmetric mode, the halo-updates do not touch
any edge values of the computational domain so we can specify values (via a
radiation routine) and then rely on the values being retained. It seems to
me we could find every halo-update and re-apply the radiation values but it
does seem like more work. I've bcc'd Mehmet Ilicak who once figured all
this out and explained it to us (convincingly then) even though I'm not
sure I am explaining it correctly now. Hopefully he'll chime in.

Flather is meant to refer to the barotropic radiation condition, yes? I
think the renaming idea was to cleanly separate the barotropic radiation
code from other [drafts] of radiation conditions which might have borrowed
the parameters/flags/code. My memory is poor so I might be confused about
what we discussed.

-A.

Dr Alistair Adcroft (Alistair.Adcroft@noaa.gov <mailto:
Alistair.Adcroft@noaa.gov>)
Princeton University Tel: (609) 987-5073
NOAA/GFDL, 201 Forrestal Road, Princeton, NJ 08540

On Wed, May 11, 2016 at 6:50 PM, Kate Hedstrom <notifications@github.com
mailto:notifications@github.com> wrote:
We are agreed that the open boundary points lie within the computational
domain, aren't we? There is evidence that this was not the plan for the
"Flather" points in the past, but I can clean it out as I find it.

If the open boundaries are inside the domain, I don't know why symmetric
mode would be needed.

Are we changing the Flather name?


You are receiving this because you are subscribed to this thread.
Reply to this email directly or view it on GitHub <
https://github.com/NOAA-GFDL/MOM6/issues/293>


You are receiving this because you authored the thread.
Reply to this email directly or view it on GitHub
NOAA-GFDL#293 (comment)

@MJHarrison-GFDL
Copy link
Contributor

Hi Kate,

Sorry , but I'm catching up here

In symmetric mode, q points

On Wed, May 25, 2016 at 3:30 PM, Kate Hedstrom notifications@github.com
wrote:

To summarize from the latest hangout and other thoughts on open boundaries:

  • The model boundary is some line through internal q points.
  • h points just outside the line are also part of the boundary conditions.
  • I'd like to use the land mask to separate wall parts of the boundary from
    open parts of the boundary. The open parts need valid h values out there,
    which isn't the case inside the land mask anyway.
  • These external h points will not be involved in the regular timestepping
  • perhaps a boundary mask can be used to shut that off. I got circle_obcs
    to run by using the OBC_SIMPLE with u,v=0 to shut off continuity out there
  • not exactly the general solution we want.
  • If the open boundary is more than one grid from the edge, points beyond
    the line of open boundary points can be marked as land.
  • I plan to have a setup which doesn't need messing with the user directory
    where the open boundaries are one grid in from the edge.

For the naming, I'm sure we'll have to go beyond simple Orlanski. The ROMS
world calls the barotropic "Chapman-Flather", the baroclinic
"radiation-nudging" (or Marchesiello). Sometimes Shchepetkin is mentioned
too. ROMS allows each field and each side to have various settings - T and
S are clamped east, radiation north, with nudging south, while dye sees a
wall all around. I don't know how general we need to be.

Kate

On Wed, May 25, 2016 at 4:22 AM, Alistair Adcroft <
notifications@github.com>
wrote:

Hi,

Yes, Flather is a specific condition of Orlanski type radiation condition
for specifically baratropic case since in the baratropic model the wave
speed is sqrt(gh).
I would prefer to keep that name. Because other models such as ROMS is
using the same nomenclature. It would be helpful for people if they come
from another model environment.
Simply baroclinic open boundary conditions are Orlanski type radiation
and
Barotropic is Flather.

I think the update was the issue in non-symmetric mode for open BC.
That’s
why it did;t work with that.
The open BC has to be on the physical domain. If I remember correctly the
wave will bounce back if you change that.

Best,
Mehmet

On May 13, 2016, at 2:34 AM, Alistair Adcroft - NOAA Affiliate <
alistair.adcroft@noaa.gov> wrote:

I think the OB points can legitimately live on the edge of the
[physical] domain but only when the model is in symmetric mode because
the
western/southern boundaries do not "belong" on any processor in
non-symmetric mode.

We've had a long standing question about why symmetric is required for
OBCs and I think it not related to the above point but more to do with
the
timing of when things are halo-updated. Halo-updates in non-symmetric
mode
would clobber any values imposed by a radiation condition on a
western/southern edge (referring to potentially interior tiles in a
parallel decomposition). In symmetric mode, the halo-updates do not touch
any edge values of the computational domain so we can specify values
(via a
radiation routine) and then rely on the values being retained. It seems
to
me we could find every halo-update and re-apply the radiation values but
it
does seem like more work. I've bcc'd Mehmet Ilicak who once figured all
this out and explained it to us (convincingly then) even though I'm not
sure I am explaining it correctly now. Hopefully he'll chime in.

Flather is meant to refer to the barotropic radiation condition, yes? I
think the renaming idea was to cleanly separate the barotropic radiation
code from other [drafts] of radiation conditions which might have
borrowed
the parameters/flags/code. My memory is poor so I might be confused about
what we discussed.

-A.

Dr Alistair Adcroft (Alistair.Adcroft@noaa.gov <mailto:
Alistair.Adcroft@noaa.gov>)
Princeton University Tel: (609) 987-5073
NOAA/GFDL, 201 Forrestal Road, Princeton, NJ 08540

On Wed, May 11, 2016 at 6:50 PM, Kate Hedstrom <
notifications@github.com
mailto:notifications@github.com> wrote:
We are agreed that the open boundary points lie within the
computational
domain, aren't we? There is evidence that this was not the plan for the
"Flather" points in the past, but I can clean it out as I find it.

If the open boundaries are inside the domain, I don't know why
symmetric
mode would be needed.

Are we changing the Flather name?


You are receiving this because you are subscribed to this thread.
Reply to this email directly or view it on GitHub <
https://github.com/NOAA-GFDL/MOM6/issues/293>


You are receiving this because you authored the thread.
Reply to this email directly or view it on GitHub
NOAA-GFDL#293 (comment)


You are receiving this because you are subscribed to this thread.
Reply to this email directly or view it on GitHub
NOAA-GFDL#293 (comment)

@MJHarrison-GFDL
Copy link
Contributor

Oops, sorry. Still getting used to my Chromebook.

I was going to discuss the question of whether the OBC locations need to be
2:N-1 instead of 1:N

Upon further thought ... I will try to get some clarification from Bob or
Alistair tomorrow.

Regarding the naming convention, I think we should stick with generic names
which reflect the wave-type , e.g. barotropic, baroclinic, normal, oblique.

Matt

On Wed, May 25, 2016 at 11:36 PM, Matthew Harrison - NOAA Federal <
matthew.harrison@noaa.gov> wrote:

Hi Kate,

Sorry , but I'm catching up here

In symmetric mode, q points

On Wed, May 25, 2016 at 3:30 PM, Kate Hedstrom notifications@github.com
wrote:

To summarize from the latest hangout and other thoughts on open
boundaries:

  • The model boundary is some line through internal q points.
  • h points just outside the line are also part of the boundary conditions.
  • I'd like to use the land mask to separate wall parts of the boundary
    from
    open parts of the boundary. The open parts need valid h values out there,
    which isn't the case inside the land mask anyway.
  • These external h points will not be involved in the regular timestepping
  • perhaps a boundary mask can be used to shut that off. I got circle_obcs
    to run by using the OBC_SIMPLE with u,v=0 to shut off continuity out there
  • not exactly the general solution we want.
  • If the open boundary is more than one grid from the edge, points beyond
    the line of open boundary points can be marked as land.
  • I plan to have a setup which doesn't need messing with the user
    directory
    where the open boundaries are one grid in from the edge.

For the naming, I'm sure we'll have to go beyond simple Orlanski. The ROMS
world calls the barotropic "Chapman-Flather", the baroclinic
"radiation-nudging" (or Marchesiello). Sometimes Shchepetkin is mentioned
too. ROMS allows each field and each side to have various settings - T and
S are clamped east, radiation north, with nudging south, while dye sees a
wall all around. I don't know how general we need to be.

Kate

On Wed, May 25, 2016 at 4:22 AM, Alistair Adcroft <
notifications@github.com>
wrote:

Hi,

Yes, Flather is a specific condition of Orlanski type radiation
condition
for specifically baratropic case since in the baratropic model the wave
speed is sqrt(gh).
I would prefer to keep that name. Because other models such as ROMS is
using the same nomenclature. It would be helpful for people if they come
from another model environment.
Simply baroclinic open boundary conditions are Orlanski type radiation
and
Barotropic is Flather.

I think the update was the issue in non-symmetric mode for open BC.
That’s
why it did;t work with that.
The open BC has to be on the physical domain. If I remember correctly
the
wave will bounce back if you change that.

Best,
Mehmet

On May 13, 2016, at 2:34 AM, Alistair Adcroft - NOAA Affiliate <
alistair.adcroft@noaa.gov> wrote:

I think the OB points can legitimately live on the edge of the
[physical] domain but only when the model is in symmetric mode because
the
western/southern boundaries do not "belong" on any processor in
non-symmetric mode.

We've had a long standing question about why symmetric is required for
OBCs and I think it not related to the above point but more to do with
the
timing of when things are halo-updated. Halo-updates in non-symmetric
mode
would clobber any values imposed by a radiation condition on a
western/southern edge (referring to potentially interior tiles in a
parallel decomposition). In symmetric mode, the halo-updates do not
touch
any edge values of the computational domain so we can specify values
(via a
radiation routine) and then rely on the values being retained. It seems
to
me we could find every halo-update and re-apply the radiation values
but it
does seem like more work. I've bcc'd Mehmet Ilicak who once figured all
this out and explained it to us (convincingly then) even though I'm not
sure I am explaining it correctly now. Hopefully he'll chime in.

Flather is meant to refer to the barotropic radiation condition, yes?
I
think the renaming idea was to cleanly separate the barotropic radiation
code from other [drafts] of radiation conditions which might have
borrowed
the parameters/flags/code. My memory is poor so I might be confused
about
what we discussed.

-A.

Dr Alistair Adcroft (Alistair.Adcroft@noaa.gov <mailto:
Alistair.Adcroft@noaa.gov>)
Princeton University Tel: (609) 987-5073
NOAA/GFDL, 201 Forrestal Road, Princeton, NJ 08540

On Wed, May 11, 2016 at 6:50 PM, Kate Hedstrom <
notifications@github.com
mailto:notifications@github.com> wrote:
We are agreed that the open boundary points lie within the
computational
domain, aren't we? There is evidence that this was not the plan for the
"Flather" points in the past, but I can clean it out as I find it.

If the open boundaries are inside the domain, I don't know why
symmetric
mode would be needed.

Are we changing the Flather name?


You are receiving this because you are subscribed to this thread.
Reply to this email directly or view it on GitHub <
https://github.com/NOAA-GFDL/MOM6/issues/293>


You are receiving this because you authored the thread.
Reply to this email directly or view it on GitHub
NOAA-GFDL#293 (comment)


You are receiving this because you are subscribed to this thread.
Reply to this email directly or view it on GitHub
NOAA-GFDL#293 (comment)

@StephenGriffies
Copy link
Contributor

Hi

I agree with the names; they should be generic. I also think they should
be prefixed with something like "obc_" to better allow for grepping for
anything related to OBC code.

Steve

On Thu, May 26, 2016 at 1:43 PM, Matthew Harrison notifications@github.com
wrote:

Oops, sorry. Still getting used to my Chromebook.

I was going to discuss the question of whether the OBC locations need to be
2:N-1 instead of 1:N

Upon further thought ... I will try to get some clarification from Bob or
Alistair tomorrow.

Regarding the naming convention, I think we should stick with generic names
which reflect the wave-type , e.g. barotropic, baroclinic, normal, oblique.

Matt

On Wed, May 25, 2016 at 11:36 PM, Matthew Harrison - NOAA Federal <
matthew.harrison@noaa.gov> wrote:

Hi Kate,

Sorry , but I'm catching up here

In symmetric mode, q points

On Wed, May 25, 2016 at 3:30 PM, Kate Hedstrom <notifications@github.com

wrote:

To summarize from the latest hangout and other thoughts on open
boundaries:

  • The model boundary is some line through internal q points.
  • h points just outside the line are also part of the boundary
    conditions.
  • I'd like to use the land mask to separate wall parts of the boundary
    from
    open parts of the boundary. The open parts need valid h values out
    there,
    which isn't the case inside the land mask anyway.
  • These external h points will not be involved in the regular
    timestepping
  • perhaps a boundary mask can be used to shut that off. I got
    circle_obcs
    to run by using the OBC_SIMPLE with u,v=0 to shut off continuity out
    there
  • not exactly the general solution we want.
  • If the open boundary is more than one grid from the edge, points
    beyond
    the line of open boundary points can be marked as land.
  • I plan to have a setup which doesn't need messing with the user
    directory
    where the open boundaries are one grid in from the edge.

For the naming, I'm sure we'll have to go beyond simple Orlanski. The
ROMS
world calls the barotropic "Chapman-Flather", the baroclinic
"radiation-nudging" (or Marchesiello). Sometimes Shchepetkin is
mentioned
too. ROMS allows each field and each side to have various settings - T
and
S are clamped east, radiation north, with nudging south, while dye sees
a
wall all around. I don't know how general we need to be.

Kate

On Wed, May 25, 2016 at 4:22 AM, Alistair Adcroft <
notifications@github.com>
wrote:

Hi,

Yes, Flather is a specific condition of Orlanski type radiation
condition
for specifically baratropic case since in the baratropic model the
wave
speed is sqrt(gh).
I would prefer to keep that name. Because other models such as ROMS is
using the same nomenclature. It would be helpful for people if they
come
from another model environment.
Simply baroclinic open boundary conditions are Orlanski type radiation
and
Barotropic is Flather.

I think the update was the issue in non-symmetric mode for open BC.
That’s
why it did;t work with that.
The open BC has to be on the physical domain. If I remember correctly
the
wave will bounce back if you change that.

Best,
Mehmet

On May 13, 2016, at 2:34 AM, Alistair Adcroft - NOAA Affiliate <
alistair.adcroft@noaa.gov> wrote:

I think the OB points can legitimately live on the edge of the
[physical] domain but only when the model is in symmetric mode because
the
western/southern boundaries do not "belong" on any processor in
non-symmetric mode.

We've had a long standing question about why symmetric is required
for
OBCs and I think it not related to the above point but more to do with
the
timing of when things are halo-updated. Halo-updates in non-symmetric
mode
would clobber any values imposed by a radiation condition on a
western/southern edge (referring to potentially interior tiles in a
parallel decomposition). In symmetric mode, the halo-updates do not
touch
any edge values of the computational domain so we can specify values
(via a
radiation routine) and then rely on the values being retained. It
seems
to
me we could find every halo-update and re-apply the radiation values
but it
does seem like more work. I've bcc'd Mehmet Ilicak who once figured
all
this out and explained it to us (convincingly then) even though I'm
not
sure I am explaining it correctly now. Hopefully he'll chime in.

Flather is meant to refer to the barotropic radiation condition,
yes?
I
think the renaming idea was to cleanly separate the barotropic
radiation
code from other [drafts] of radiation conditions which might have
borrowed
the parameters/flags/code. My memory is poor so I might be confused
about
what we discussed.

-A.

Dr Alistair Adcroft (Alistair.Adcroft@noaa.gov <mailto:
Alistair.Adcroft@noaa.gov>)
Princeton University Tel: (609) 987-5073
NOAA/GFDL, 201 Forrestal Road, Princeton, NJ 08540

On Wed, May 11, 2016 at 6:50 PM, Kate Hedstrom <
notifications@github.com
mailto:notifications@github.com> wrote:
We are agreed that the open boundary points lie within the
computational
domain, aren't we? There is evidence that this was not the plan for
the
"Flather" points in the past, but I can clean it out as I find it.

If the open boundaries are inside the domain, I don't know why
symmetric
mode would be needed.

Are we changing the Flather name?


You are receiving this because you are subscribed to this thread.
Reply to this email directly or view it on GitHub <
https://github.com/NOAA-GFDL/MOM6/issues/293>


You are receiving this because you authored the thread.
Reply to this email directly or view it on GitHub
NOAA-GFDL#293 (comment)


You are receiving this because you are subscribed to this thread.
Reply to this email directly or view it on GitHub
NOAA-GFDL#293 (comment)


You are receiving this because you are subscribed to this thread.
Reply to this email directly or view it on GitHub
NOAA-GFDL#293 (comment)

Dr. Stephen M. Griffies
NOAA Geophysical Fluid Dynamics Lab
201 Forrestal Road
Princeton, NJ 08542
USA

@adcroft
Copy link
Collaborator

adcroft commented Jun 4, 2016

I finally had time to look through the OBC code myself. I've not discussed this response with @Hallberg-NOAA yet, who we know has the better understanding of the OBCs as they exist, but here are my comments anyway - they might need to be corrected by Bob.

@kshedstrom wrote:

  • The model boundary is some line through internal q points.

To be precise, the model boundary must be through q-points in the global computational domain of q, which does include the edges of the physical domain. Here's the rationale (I know you know this but this is for the wider audience):

Speaking in global terms only (no decomposition) we consider the physical domain to be defined by the volume of all the h cells and the edges of those h-volumes are part of the physical domain in that they are boundaries. This then includes (as part of the physical domain) the q-points which are on the boundary of the physical domain. In global indices, h would be declared h(isg:ieg,jsg:jeg). See comment on global indices for fuller explanation of isg, ieg, etc. but isg=1, ieg=25 for circle_obcs. In symmetric memory mode, the q-points would be declared q(isg-1:ieg,jsg-1:jeg), u-points would be u(isg-1:ieg,jsg:jeg), and v-points would be v(isg:ieg,jsg-1:jeg). For a solid-walled box, no normal flow would give u(isg-1,:)=0, u(ieg,:)=0, v(:,jsg-1)=0 and v(:,jeg)=0. The physical domain would only be reduced by masking some of those volumes as land-points. Any land-masked cells in the interior would have similar boundary conditions. The "zero" part makes the model immune to values outside of the physical domain (land values). An open boundary condition (OBC) on u or v has a non-zero value on the physical boundary and this brings us to your next comment about values outside of the boundary position.

  • h points just outside the line are also part of the boundary conditions.

Yes, they might be part of the boundary conditions but those points are outside of the physical domain (because they are beyond the boundary). The model expects the boundary condition code to either override the fluxes that would propagate those exterior values in, or to set those exterior values so that their inward propagation is sensible. So if the boundary condition is on the edge of the computational domain (not mandatory) we would be using a halo location for this data. In short, the data immediately outside of the boundary is expected to be set by the boundary condition.

  • I'd like to use the land mask to separate wall parts of the boundary from open parts of the boundary. The open parts need valid h values out there, which isn't the case inside the land mask anyway.

To have patches of open-boundary we need to set the elements of the arrays OBC%OBC_mask_u(:,:) and OBC%OBC_mask_v(:,:) appropriately. Currently, the flag APPLY_OBC_U_FLATHER_EAST sets OBC%OBC_mask_u(isg,:) = .true. (i.e. the entire edge) at lines 1822-1838 or MOM_state_initialization.F90. We need a more flexible way to describe where OBCs will be applied. The comment at L1820 says this (bless the comment writer).

The initial h values immediately outside the open boundaries are set to be equal to the first interior value at the bottom of that same subroutine. i.e. The initial normal gradient of h is set to zero by setting the value immediately outside the physical domain.

The bathymetry and masks are similarly set up to be consistent with flather. Halos of G%bathyT are first filled in from the neighboring PEs at line 1382 of MOM_grid_initialize.F90 and then lines 1388-1394 set G%bathyT(isg-1,:)=G%bathyT(isg,:), i.e. d/dx D = 0 alond the entire edge of the domain. Unfortunately again, there is a comment saying this should be generalized (bless the comment writer).

For circle_obcs, a square box, the above works and yields a T-mask shown below. Note that the highlighted values in the halo region (halo width=4) have been "opened up" so the model does not treat the boundary as solid.
image

To generalize the above bits of code I've mentioned, we need to generalize the way OBC%OBC_mask_u(:,:) is set and to only set 'G%bathyT(i-1,j)=G%bathyT(i,j) only where OBC%OBC_mask_u(i-1,j) = .true.. The problem is that OBC%OBC_mask_u(:,:) is currently set up AFTER the bathymetry/masks have been set up. So the first thing we need to do is to break out the setup of OBC positions and call their set up at the point where they should subsequently update the bathymetry/masks.

  • These external h points will not be involved in the regular time stepping
  • perhaps a boundary mask can be used to shut that off. I got circle_obcs to run by using the OBC_SIMPLE with u,v=0 to shut off continuity out there
  • not exactly the general solution we want.

These points are time-stepped but the results are ignored or overridden by the OBC calls. For example, where the zonal mass flux has just been calculated, ignoring the presence of OBC, lines 442-444 of MOM_continuity_PPM.F90 then replace the value of the mass flux with the prescribed value.

  • I plan to have a setup which doesn't need messing with the user directory
    where the open boundaries are one grid in from the edge.

Yes. We need to design a good way to describe to the model where open boundaries will be located, possibly by just some input parameters rather than binary files. Here's a suggestion we've recently thought about:
unnamed

The above sketched open boundary locations can be described as a set of numbered segments which can be described to the model by a few short lists of indices. With these indices the model could figure out the maps of OBC_mask_u and OBC_mask_v as needed. The only restriction is that each segment must be either a line of constant i or constant j. Hence the three segments 5, 6, 7 are separate segments even though physically joined. This segment labelling also provides a rational way to make the input data more concise - just the T(x,z,t) for each segment need be provided rather than provide the full 4D data for the entire domain. @MJHarrison-GFDL is looking into what features of mpp_io could be used to permit this.

@StephenGriffies
Copy link
Contributor

Can we please please please bring these comments into a form that is more
permanent than email? This sort of material should go into an online
manual for the OBC code.

On Sat, Jun 4, 2016 at 4:46 PM, Alistair Adcroft notifications@github.com
wrote:

I finally had time to look through the OBC code myself. I've not discussed
this response with @Hallberg-NOAA https://github.com/Hallberg-NOAA yet,
who we know has the better understanding of the OBCs as they exist, but
here are my comments anyway - they might need to be corrected by Bob.

@kshedstrom https://github.com/kshedstrom wrote:

  • The model boundary is some line through internal q points.

To be precise, the model boundary must be through q-points in the global
computational domain of q, which does include the edges of the physical
domain. Here's the rationale (I know you know this but this is for the
wider audience):

Speaking in global terms only (no decomposition) we consider the
physical domain to be defined by the volume of all the h cells and the
edges of those h-volumes are part of the physical domain in that they are
boundaries. This then includes (as part of the physical domain) the
q-points which are on the boundary of the physical domain. In global
indices, h would be declared h(isg:ieg,jsg:jeg). See comment on global
indices
NOAA-GFDL/MOM6-examples#42 (comment)
for fuller explanation of isg, ieg, etc. but isg=1, ieg=25 for circle_obcs.
In symmetric memory mode, the q-points would be declared
q(isg-1:ieg,jsg-1:jeg), u-points would be u(isg-1:ieg,jsg:jeg), and
v-points would be v(isg:ieg,jsg-1:jeg). For a solid-walled box, no normal
flow would give u(isg-1,:)=0, u(ieg,:)=0, v(:,jsg-1)=0 and v(:,jeg)=0.
The ph ysical domain would only be reduced by masking some of those volumes
as land-points. Any land-masked cells in the interior would have similar
boundary conditions. The "zero" part makes the model immune to values
outside of the physical domain (land values). An open boundary condition
(OBC) on u or v has a non-zero value on the physical boundary and this
brings us to your next comment about values outside of the boundary
position.

  • h points just outside the line are also part of the boundary
    conditions.

Yes, they might be part of the boundary conditions but those points are
outside of the physical domain (because they are beyond the boundary). The
model expects the boundary condition code to either override the fluxes
that would propagate those exterior values in, or to set those exterior
values so that their inward propagation is sensible. So if the boundary
condition is on the edge of the computational domain (not mandatory) we
would be using a halo location for this data. In short, the data
immediately outside of the boundary is expected to be set by the boundary
condition.

  • I'd like to use the land mask to separate wall parts of the boundary
    from open parts of the boundary. The open parts need valid h values out
    there, which isn't the case inside the land mask anyway.

To have patches of open-boundary we need to set the elements of the arrays
OBC%OBC_mask_u(:,:) and OBC%OBC_mask_v(:,:) appropriately. Currently, the
flag APPLY_OBC_U_FLATHER_EAST sets OBC%OBC_mask_u(isg,:) = .true. (i.e.
the entire edge
) at lines 1822-1838 or MOM_state_initialization.F90
https://github.com/NOAA-GFDL/MOM6/blob/4571bffbb031e5bab93da05796b30cb619c3424e/src/initialization/MOM_state_initialization.F90#L1822-L1838.
We need a more flexible way to describe where OBCs will be applied. The
comment at L1820
https://github.com/NOAA-GFDL/MOM6/blob/4571bffbb031e5bab93da05796b30cb619c3424e/src/initialization/MOM_state_initialization.F90#L1820
says this (bless the comment writer).

The initial h values immediately outside the open boundaries are set to
be equal to the first interior value at the bottom of that same subroutine.
i.e. The initial normal gradient of h is set to zero by setting the value
immediately outside the physical domain.

The bathymetry and masks are similarly set up to be consistent with
flather. Halos of G%bathyT are first filled in from the neighboring PEs
at line 1382 of MOM_grid_initialize.F90
https://github.com/NOAA-GFDL/MOM6/blob/4571bffbb031e5bab93da05796b30cb619c3424e/src/initialization/MOM_grid_initialize.F90#L1382
and then lines 1388-1394
https://github.com/NOAA-GFDL/MOM6/blob/4571bffbb031e5bab93da05796b30cb619c3424e/src/initialization/MOM_grid_initialize.F90#L1388-L1394
set G%bathyT(isg-1,:)=G%bathyT(isg,:), i.e. d/dx D = 0 alond the entire
edge of the domain. Unfortunately again, there is a comment saying this
should be generalized (bless the comment writer).

For circle_obcs, a square box, the above works and yields a T-mask shown
below. Note that the highlighted values in the halo region (halo width=4)
have been "opened up" so the model does not treat the boundary as solid.
[image: image]
https://cloud.githubusercontent.com/assets/5859571/15801446/eec7eb34-2a62-11e6-87b3-52e3cdd01ee0.png

To generalize the above bits of code I've mentioned, we need to generalize
the way OBC%OBC_mask_u(:,:) is set and to only set '
G%bathyT(i-1,j)=G%bathyT(i,j) only where OBC%OBC_mask_u(i-1,j) = .true..
The problem is that OBC%OBC_mask_u(:,:) is currently set up AFTER the
bathymetry/masks have been set up. So the first thing we need to do is to
break out the setup of OBC positions and call their set up at the point
where they should subsequently update the bathymetry/masks.

  • These external h points will not be involved in the regular time
    stepping
  • perhaps a boundary mask can be used to shut that off. I got
    circle_obcs to run by using the OBC_SIMPLE with u,v=0 to shut off
    continuity out there
  • not exactly the general solution we want.

These points are time-stepped but the results are ignored or overridden by
the OBC calls. For example, where the zonal mass flux has just been
calculated, ignoring the presence of OBC, lines 442-444 of
MOM_continuity_PPM.F90
https://github.com/NOAA-GFDL/MOM6/blob/4571bffbb031e5bab93da05796b30cb619c3424e/src/core/MOM_continuity_PPM.F90#L442-L444
then replace the value of the mass flux with the prescribed value.

  • I plan to have a setup which doesn't need messing with the user
    directory where the open boundaries are one grid in from the edge.

Yes. We need to design a good way to describe to the model where open
boundaries will be located, possibly by just some input parameters rather
than binary files. Here's a suggestion we've recently thought about:
[image: unnamed]
https://cloud.githubusercontent.com/assets/5859571/15801848/bfd69106-2a6f-11e6-9d52-c49133165a78.png

The above sketched open boundary locations can be described as a set of
numbered segments which can be described to the model by a few short lists
of indices. With these indices the model could figure out the maps of
OBC_mask_u and OBC_mask_v as needed. The only restriction is that each
segment must be either a line of constant i or constant j. Hence the
three segments 5, 6, 7 are separate segments even though physically joined.
This segment labelling also provides a rational way to make the input data
more concise - just the T(x,z,t) for each segment need be provided rather
than provide the full 4D data for the entire domain. @MJHarrison-GFDL
https://github.com/MJHarrison-GFDL is looking into what features of
mpp_io could be used to permit this.


You are receiving this because you commented.
Reply to this email directly, view it on GitHub
NOAA-GFDL#293 (comment), or mute
the thread
https://github.com/notifications/unsubscribe/ACBam87t0j3BeXj8oUImIYslehwA_B7nks5qIeQdgaJpZM4IckNe
.

Dr. Stephen M. Griffies
NOAA Geophysical Fluid Dynamics Lab
201 Forrestal Road
Princeton, NJ 08542
USA

@MJHarrison-GFDL
Copy link
Contributor

Thanks @adcroft,

My two cents:

In addition to the end points, we need to specify the orientation of the
OBC, i.e. whether the OBC_mask points are at +/-1

The discarded part of the domain should be flooded via ICE9, setting the h
values to zero.

There are two options for filling in the OBC data: (A) supply gridded data
which encompasses the domain of the OBC segments, and passing through
space-time interpolation routine at runtime (B) supply preprocessed data
(tracers and transports) for each segment. (A) of course eliminates the
need for a preprocessing step, but suffers from potentially HUGE volumes of
data . (B) simplifies the code - requiring time interpolation only, but
introduces additional preprocessing steps (which we have been trying to
avoid).

Matt

On Sat, Jun 4, 2016 at 4:46 PM, Alistair Adcroft notifications@github.com
wrote:

I finally had time to look through the OBC code myself. I've not discussed
this response with @Hallberg-NOAA https://github.com/Hallberg-NOAA yet,
who we know has the better understanding of the OBCs as they exist, but
here are my comments anyway - they might need to be corrected by Bob.

@kshedstrom https://github.com/kshedstrom wrote:

  • The model boundary is some line through internal q points.

To be precise, the model boundary must be through q-points in the global
computational domain of q, which does include the edges of the physical
domain. Here's the rationale (I know you know this but this is for the
wider audience):

Speaking in global terms only (no decomposition) we consider the
physical domain to be defined by the volume of all the h cells and the
edges of those h-volumes are part of the physical domain in that they are
boundaries. This then includes (as part of the physical domain) the
q-points which are on the boundary of the physical domain. In global
indices, h would be declared h(isg:ieg,jsg:jeg). See comment on global
indices
NOAA-GFDL/MOM6-examples#42 (comment)
for fuller explanation of isg, ieg, etc. but isg=1, ieg=25 for circle_obcs.
In symmetric memory mode, the q-points would be declared
q(isg-1:ieg,jsg-1:jeg), u-points would be u(isg-1:ieg,jsg:jeg), and
v-points would be v(isg:ieg,jsg-1:jeg). For a solid-walled box, no normal
flow would give u(isg-1,:)=0, u(ieg,:)=0, v(:,jsg-1)=0 and v(:,jeg)=0.
The ph ysical domain would only be reduced by masking some of those volumes
as land-points. Any land-masked cells in the interior would have similar
boundary conditions. The "zero" part makes the model immune to values
outside of the physical domain (land values). An open boundary condition
(OBC) on u or v has a non-zero value on the physical boundary and this
brings us to your next comment about values outside of the boundary
position.

  • h points just outside the line are also part of the boundary
    conditions.

Yes, they might be part of the boundary conditions but those points are
outside of the physical domain (because they are beyond the boundary). The
model expects the boundary condition code to either override the fluxes
that would propagate those exterior values in, or to set those exterior
values so that their inward propagation is sensible. So if the boundary
condition is on the edge of the computational domain (not mandatory) we
would be using a halo location for this data. In short, the data
immediately outside of the boundary is expected to be set by the boundary
condition.

  • I'd like to use the land mask to separate wall parts of the boundary
    from open parts of the boundary. The open parts need valid h values out
    there, which isn't the case inside the land mask anyway.

To have patches of open-boundary we need to set the elements of the arrays
OBC%OBC_mask_u(:,:) and OBC%OBC_mask_v(:,:) appropriately. Currently, the
flag APPLY_OBC_U_FLATHER_EAST sets OBC%OBC_mask_u(isg,:) = .true. (i.e.
the entire edge
) at lines 1822-1838 or MOM_state_initialization.F90
https://github.com/NOAA-GFDL/MOM6/blob/4571bffbb031e5bab93da05796b30cb619c3424e/src/initialization/MOM_state_initialization.F90#L1822-L1838.
We need a more flexible way to describe where OBCs will be applied. The
comment at L1820
https://github.com/NOAA-GFDL/MOM6/blob/4571bffbb031e5bab93da05796b30cb619c3424e/src/initialization/MOM_state_initialization.F90#L1820
says this (bless the comment writer).

The initial h values immediately outside the open boundaries are set to
be equal to the first interior value at the bottom of that same subroutine.
i.e. The initial normal gradient of h is set to zero by setting the value
immediately outside the physical domain.

The bathymetry and masks are similarly set up to be consistent with
flather. Halos of G%bathyT are first filled in from the neighboring PEs
at line 1382 of MOM_grid_initialize.F90
https://github.com/NOAA-GFDL/MOM6/blob/4571bffbb031e5bab93da05796b30cb619c3424e/src/initialization/MOM_grid_initialize.F90#L1382
and then lines 1388-1394
https://github.com/NOAA-GFDL/MOM6/blob/4571bffbb031e5bab93da05796b30cb619c3424e/src/initialization/MOM_grid_initialize.F90#L1388-L1394
set G%bathyT(isg-1,:)=G%bathyT(isg,:), i.e. d/dx D = 0 alond the entire
edge of the domain. Unfortunately again, there is a comment saying this
should be generalized (bless the comment writer).

For circle_obcs, a square box, the above works and yields a T-mask shown
below. Note that the highlighted values in the halo region (halo width=4)
have been "opened up" so the model does not treat the boundary as solid.
[image: image]
https://cloud.githubusercontent.com/assets/5859571/15801446/eec7eb34-2a62-11e6-87b3-52e3cdd01ee0.png

To generalize the above bits of code I've mentioned, we need to generalize
the way OBC%OBC_mask_u(:,:) is set and to only set '
G%bathyT(i-1,j)=G%bathyT(i,j) only where OBC%OBC_mask_u(i-1,j) = .true..
The problem is that OBC%OBC_mask_u(:,:) is currently set up AFTER the
bathymetry/masks have been set up. So the first thing we need to do is to
break out the setup of OBC positions and call their set up at the point
where they should subsequently update the bathymetry/masks.

  • These external h points will not be involved in the regular time
    stepping
  • perhaps a boundary mask can be used to shut that off. I got
    circle_obcs to run by using the OBC_SIMPLE with u,v=0 to shut off
    continuity out there
  • not exactly the general solution we want.

These points are time-stepped but the results are ignored or overridden by
the OBC calls. For example, where the zonal mass flux has just been
calculated, ignoring the presence of OBC, lines 442-444 of
MOM_continuity_PPM.F90
https://github.com/NOAA-GFDL/MOM6/blob/4571bffbb031e5bab93da05796b30cb619c3424e/src/core/MOM_continuity_PPM.F90#L442-L444
then replace the value of the mass flux with the prescribed value.

  • I plan to have a setup which doesn't need messing with the user
    directory where the open boundaries are one grid in from the edge.

Yes. We need to design a good way to describe to the model where open
boundaries will be located, possibly by just some input parameters rather
than binary files. Here's a suggestion we've recently thought about:
[image: unnamed]
https://cloud.githubusercontent.com/assets/5859571/15801848/bfd69106-2a6f-11e6-9d52-c49133165a78.png

The above sketched open boundary locations can be described as a set of
numbered segments which can be described to the model by a few short lists
of indices. With these indices the model could figure out the maps of
OBC_mask_u and OBC_mask_v as needed. The only restriction is that each
segment must be either a line of constant i or constant j. Hence the
three segments 5, 6, 7 are separate segments even though physically joined.
This segment labelling also provides a rational way to make the input data
more concise - just the T(x,z,t) for each segment need be provided rather
than provide the full 4D data for the entire domain. @MJHarrison-GFDL
https://github.com/MJHarrison-GFDL is looking into what features of
mpp_io could be used to permit this.


You are receiving this because you were mentioned.
Reply to this email directly, view it on GitHub
NOAA-GFDL#293 (comment), or mute
the thread
https://github.com/notifications/unsubscribe/AGT88MR295ocKOhzDAAZzm_1o1QIf7bSks5qIeQegaJpZM4IckNe
.

@kshedstrom
Copy link
Collaborator Author

Here's an example of something we've run. There are open boundaries to the
south and west with values coming from a parent grid at hourly resolution -
this resolves the tides and avoids the need for driving this domain with
additional tidal information. However, we are also feeding it 1667 "simple"
boundary values with daily fields coming from a hydrology model, each being
one grid cell in length. As @MJHarrison-NOAA notes, each land-sea edge
needs to know its orientation. This is so far in the preprocessing camp.
ROMS by default expects tracer values for each depth for each river, but I
introduced an option whereby each river has the same tracer values (they
can still be a function of time).

[image: Inline image 1]

For @adcroft, your band of highlighted T points is part of the ROMS
(symmetric) version of the "whole grid". That's why we can count on knowing
the land mask values out there. For CCS1, the western edge can be
considered to be at isg-1 because there are no masked points along it.

Our final solution is going to require things that are currently part of
Flather as well as things that are in the other boundary code. In ROMS, the
biharmonic operators need an extra boundary condition on the intermediate
Laplacian of the field, setting gradients to zero. I don't know off-hand
what is done for other terms in the equations requiring a stencil larger
than +/- 1.

Kate

On Sun, Jun 5, 2016 at 8:28 AM, Matthew Harrison notifications@github.com
wrote:

Thanks @adcroft,

My two cents:

In addition to the end points, we need to specify the orientation of the
OBC, i.e. whether the OBC_mask points are at +/-1

The discarded part of the domain should be flooded via ICE9, setting the h
values to zero.

There are two options for filling in the OBC data: (A) supply gridded data
which encompasses the domain of the OBC segments, and passing through
space-time interpolation routine at runtime (B) supply preprocessed data
(tracers and transports) for each segment. (A) of course eliminates the
need for a preprocessing step, but suffers from potentially HUGE volumes of
data . (B) simplifies the code - requiring time interpolation only, but
introduces additional preprocessing steps (which we have been trying to
avoid).

Matt

On Sat, Jun 4, 2016 at 4:46 PM, Alistair Adcroft <notifications@github.com

wrote:

I finally had time to look through the OBC code myself. I've not
discussed
this response with @Hallberg-NOAA https://github.com/Hallberg-NOAA
yet,
who we know has the better understanding of the OBCs as they exist, but
here are my comments anyway - they might need to be corrected by Bob.

@kshedstrom https://github.com/kshedstrom wrote:

  • The model boundary is some line through internal q points.

To be precise, the model boundary must be through q-points in the global
computational domain of q, which does include the edges of the physical
domain. Here's the rationale (I know you know this but this is for the
wider audience):

Speaking in global terms only (no decomposition) we consider the
physical domain to be defined by the volume of all the h cells and the
edges of those h-volumes are part of the physical domain in that they are
boundaries. This then includes (as part of the physical domain) the
q-points which are on the boundary of the physical domain. In global
indices, h would be declared h(isg:ieg,jsg:jeg). See comment on global
indices
<
NOAA-GFDL/MOM6-examples#42 (comment)

for fuller explanation of isg, ieg, etc. but isg=1, ieg=25 for
circle_obcs.
In symmetric memory mode, the q-points would be declared
q(isg-1:ieg,jsg-1:jeg), u-points would be u(isg-1:ieg,jsg:jeg), and
v-points would be v(isg:ieg,jsg-1:jeg). For a solid-walled box, no normal
flow would give u(isg-1,:)=0, u(ieg,:)=0, v(:,jsg-1)=0 and v(:,jeg)=0.
The ph ysical domain would only be reduced by masking some of those
volumes
as land-points. Any land-masked cells in the interior would have similar
boundary conditions. The "zero" part makes the model immune to values
outside of the physical domain (land values). An open boundary condition
(OBC) on u or v has a non-zero value on the physical boundary and this
brings us to your next comment about values outside of the boundary
position.

  • h points just outside the line are also part of the boundary
    conditions.

Yes, they might be part of the boundary conditions but those points are
outside of the physical domain (because they are beyond the boundary).
The
model expects the boundary condition code to either override the fluxes
that would propagate those exterior values in, or to set those exterior
values so that their inward propagation is sensible. So if the boundary
condition is on the edge of the computational domain (not mandatory) we
would be using a halo location for this data. In short, the data
immediately outside of the boundary is expected to be set by the boundary
condition.

  • I'd like to use the land mask to separate wall parts of the boundary
    from open parts of the boundary. The open parts need valid h values out
    there, which isn't the case inside the land mask anyway.

To have patches of open-boundary we need to set the elements of the
arrays
OBC%OBC_mask_u(:,:) and OBC%OBC_mask_v(:,:) appropriately. Currently, the
flag APPLY_OBC_U_FLATHER_EAST sets OBC%OBC_mask_u(isg,:) = .true. (i.e.
the entire edge
) at lines 1822-1838 or MOM_state_initialization.F90
<
https://github.com/NOAA-GFDL/MOM6/blob/4571bffbb031e5bab93da05796b30cb619c3424e/src/initialization/MOM_state_initialization.F90#L1822-L1838
.
We need a more flexible way to describe where OBCs will be applied. The
comment at L1820
<
https://github.com/NOAA-GFDL/MOM6/blob/4571bffbb031e5bab93da05796b30cb619c3424e/src/initialization/MOM_state_initialization.F90#L1820

says this (bless the comment writer).

The initial h values immediately outside the open boundaries are set to
be equal to the first interior value at the bottom of that same
subroutine.
i.e. The initial normal gradient of h is set to zero by setting the value
immediately outside the physical domain.

The bathymetry and masks are similarly set up to be consistent with
flather. Halos of G%bathyT are first filled in from the neighboring PEs
at line 1382 of MOM_grid_initialize.F90
<
https://github.com/NOAA-GFDL/MOM6/blob/4571bffbb031e5bab93da05796b30cb619c3424e/src/initialization/MOM_grid_initialize.F90#L1382

and then lines 1388-1394
<
https://github.com/NOAA-GFDL/MOM6/blob/4571bffbb031e5bab93da05796b30cb619c3424e/src/initialization/MOM_grid_initialize.F90#L1388-L1394

set G%bathyT(isg-1,:)=G%bathyT(isg,:), i.e. d/dx D = 0 alond the entire
edge of the domain. Unfortunately again, there is a comment saying this
should be generalized (bless the comment writer).

For circle_obcs, a square box, the above works and yields a T-mask shown
below. Note that the highlighted values in the halo region (halo width=4)
have been "opened up" so the model does not treat the boundary as solid.
[image: image]
<
https://cloud.githubusercontent.com/assets/5859571/15801446/eec7eb34-2a62-11e6-87b3-52e3cdd01ee0.png

To generalize the above bits of code I've mentioned, we need to
generalize
the way OBC%OBC_mask_u(:,:) is set and to only set '
G%bathyT(i-1,j)=G%bathyT(i,j) only where OBC%OBC_mask_u(i-1,j) = .true..
The problem is that OBC%OBC_mask_u(:,:) is currently set up AFTER the
bathymetry/masks have been set up. So the first thing we need to do is to
break out the setup of OBC positions and call their set up at the point
where they should subsequently update the bathymetry/masks.

  • These external h points will not be involved in the regular time
    stepping
  • perhaps a boundary mask can be used to shut that off. I got
    circle_obcs to run by using the OBC_SIMPLE with u,v=0 to shut off
    continuity out there
  • not exactly the general solution we want.

These points are time-stepped but the results are ignored or overridden
by
the OBC calls. For example, where the zonal mass flux has just been
calculated, ignoring the presence of OBC, lines 442-444 of
MOM_continuity_PPM.F90
<
https://github.com/NOAA-GFDL/MOM6/blob/4571bffbb031e5bab93da05796b30cb619c3424e/src/core/MOM_continuity_PPM.F90#L442-L444

then replace the value of the mass flux with the prescribed value.

  • I plan to have a setup which doesn't need messing with the user
    directory where the open boundaries are one grid in from the edge.

Yes. We need to design a good way to describe to the model where open
boundaries will be located, possibly by just some input parameters rather
than binary files. Here's a suggestion we've recently thought about:
[image: unnamed]
<
https://cloud.githubusercontent.com/assets/5859571/15801848/bfd69106-2a6f-11e6-9d52-c49133165a78.png

The above sketched open boundary locations can be described as a set of
numbered segments which can be described to the model by a few short
lists
of indices. With these indices the model could figure out the maps of
OBC_mask_u and OBC_mask_v as needed. The only restriction is that each
segment must be either a line of constant i or constant j. Hence the
three segments 5, 6, 7 are separate segments even though physically
joined.
This segment labelling also provides a rational way to make the input
data
more concise - just the T(x,z,t) for each segment need be provided rather
than provide the full 4D data for the entire domain. @MJHarrison-GFDL
https://github.com/MJHarrison-GFDL is looking into what features of
mpp_io could be used to permit this.


You are receiving this because you were mentioned.
Reply to this email directly, view it on GitHub
NOAA-GFDL#293 (comment),
or mute
the thread
<
https://github.com/notifications/unsubscribe/AGT88MR295ocKOhzDAAZzm_1o1QIf7bSks5qIeQegaJpZM4IckNe

.


You are receiving this because you were mentioned.
Reply to this email directly, view it on GitHub
NOAA-GFDL#293 (comment), or mute
the thread
https://github.com/notifications/unsubscribe/AAbIHgrXD2-_cQFJW3ehIJhr7Myjgtbhks5qIvktgaJpZM4IckNe
.

@kshedstrom
Copy link
Collaborator Author

@StephenGriffies
Here's a stab at something: https://github.com/NOAA-GFDL/MOM6-examples/wiki/Open-Boundary-Conditions

@StephenGriffies
Copy link
Contributor

Great!

An image with halo issues examined will be useful too. But great to have a
place to start putting these points, with eventual cleanup when mature.

Thanks,
Steve

On Mon, Jun 6, 2016 at 8:27 PM, Kate Hedstrom notifications@github.com
wrote:

@StephenGriffies https://github.com/StephenGriffies
Here's a stab at something:
https://github.com/NOAA-GFDL/MOM6-examples/wiki/Open-Boundary-Conditions


You are receiving this because you were mentioned.
Reply to this email directly, view it on GitHub
NOAA-GFDL#293 (comment), or mute
the thread
https://github.com/notifications/unsubscribe/ACBam6JErO2OWy5ATQz6Iw0qbXMHXn63ks5qJLrngaJpZM4IckNe
.

Dr. Stephen M. Griffies
NOAA Geophysical Fluid Dynamics Lab
201 Forrestal Road
Princeton, NJ 08542
USA

@kshedstrom
Copy link
Collaborator Author

open_shut
Is this the best way to get images onto github?

Kate

marshallward pushed a commit to MJHarrison-GFDL/MOM6 that referenced this issue Jan 9, 2023
* Add GL90 parameterization in stacked shallow water

This adds a new vertical viscosity parameterization as in Greatbatch and Lamb (1990),
Ferreira & Marshall (2006) and Zhao & Vallis (2008), hereafter referred to as the GL90
vertical viscosity parameterization. This vertical viscosity scheme redistributes momentum
in the vertical, and is the equivalent of the Gent & McWilliams (1990) parameterization,
but in a TWA (thickness-weighted averaged) set of equations. The vertical viscosity coefficient
nu is computed from kappa_GM via thermal wind balance, and the following relation:

nu = kappa_GM * f^2 / N^2.

The vertical viscosity del_z ( nu del_z u) is applied to the momentum equation with stress-free
boundary conditions at the top and bottom.

In the current implementation, kappa_GM is assumed either (a) constant or as (b) having an
EBT structure. A third possible formulation of nu is depth-independent:

nu = f^2 * alpha

The latter formulation would be equivalent to a kappa_GM that varies as N^2 with depth.
Currently, the GL90 parameterization is only implemented in stacked shallow water (SSW) mode,
in which case we have 1/N^2 = h/g'.

More specifically, this commit adds a new subroutine that computes the couping coefficient
associated with GL90 via

a_cpl_gl90 = nu / h = kappa_GM * f^2 / g'

or

a_cpl_gl90 = nu / h = f^2 * alpha / h.

Further, a_cpl_gl90 is multiplied by a function (botfn), which is 0 within the GL90 bottom boundary
layer, whose depth is set by Hbbl_gl90, and 1 otherwise. This modification is necessary to avlid fluxing
momentum into vanished layers that ride over steep topography.

Finally, a_cpl_gl90 is added to a_cpl, where the latter is the coupling coefficient associated with the
remaining vertical stresses, used in the vertical viscosity solver.

More information can be found in Loose et al. (https://www.essoar.org/doi/abs/10.1002/essoar.10512867.1),
Appendix B.

New diagnostics:
* au_gl90_visc: zonal viscous coupling coefficient associated with GL90, is contained in au_visc
* av_gl90_visc: meridional viscous coupling coefficient associated with GL90, is contained in av_visc
* Kv_gl90_u: GL90 vertical viscosity at u-points, is contained in Kv_u
* Kv_gl90_v: GL90 vertical viscosity at v-points, is contained in Kv_v
* du_dt_visc_gl90: zonal acceleration due to GL90 vertical viscosity, included in du_dt_visc
* dv_dt_visc_gl90: meridional acceleration due to GL90 vertical viscosity, included in dv_dt_visc
* GLwork: Kinetic Energy Source from GL90 Vertical Viscosity

The energetics of the GL90 parameterization (named "GLwork") are intentionally computed in
MOM_vert_friction, rather than in MOM_diagnostics, where the reamining kinetic energy budget
terms are computed. We have to do the computation in MOM_vert_friction to ensure sign-
definiteness when GLwork is summed in the vertical. Indeed, MOM_diagnostics does not have
access to the velocities and thicknesses used in the vertical solver, but rather uses a time-mean
barotropic transport [uv]h to compute the energy budget diagnostics. A detailed discussion and
exploration of this issue can be found in ocean-eddy-cpt#25.

As a result of not computing the energetics in MOM_diagnostics, GLwork is not exactly contained
in KE_visc. KE_visc represents the energetics of all vertical viscosity contributions, including
the GL90 vertical viscosity. We could implement a term "KE_visc_gl90" that can be 1-to-1 compared
to KE_visc; that is, KE_visc - KE_visc_gl90 would represent exactly the energetics of all viscosity
contributions EXCEPT the GL90 viscosity. If we implemented KE_visc_gl90, this term would in practice
be very similar as GLwork, but sign-definiteness is not ensured, see above.
iangrooms pushed a commit to iangrooms/MOM6 that referenced this issue Aug 20, 2024
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

No branches or pull requests

4 participants