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

Fix Diagnostics of Diagnostics #452

Merged
merged 25 commits into from
Nov 13, 2023
Merged

Fix Diagnostics of Diagnostics #452

merged 25 commits into from
Nov 13, 2023

Conversation

Witt-D
Copy link
Collaborator

@Witt-D Witt-D commented Oct 10, 2023

I went through the thermodynamic diagnostics and added the self_setup_theromdyancs line where appropriate.
Additionatlly I added Toms hack fix for why the potential energy is failing

@@ -27,10 +27,10 @@
tmax = dt
dumpfreq = 1
else:
nlayers = 10
columns = 150
tmax = 3600.
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Did you mean to add all of these changes?

diagnostic_fields = [CourantNumber(), Gradient("u"), Perturbation('theta'),
Gradient("theta_perturbation"), Perturbation('rho'),
RichardsonNumber("theta", parameters.g/Tsurf), Gradient("theta")]
diagnostic_fields = [CourantNumber(), Pressure(eqns), SteadyStateError('Pressure_Vt')]
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Are any of the changes to this file meant to be added to this PR?

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

yes im not sure why the examples have been added as i was unaware that these had even been pushed from my local area?

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think it makes sense to add something like Pressure(eqns) and SteadyStateError('Pressure_Vt') to a test. I don't think this is the right test as we don't expect this test to preserve a steady state.

I think it would make more sense to add SteadyStateError('relative_vorticity') [not sure if this is the right name so please check!] to the Williamson 2 test.

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Removed the changes to the skamarock test and implemented it in Williamson2

# By default set this new field to the current value
# This may be overwritten if picking up from a checkpoint
field2.assign(field1)
self.init_field_set = False
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Looks like we can get rid of field2 - is that right?

Copy link
Contributor

@tommbendall tommbendall left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Thanks for this Dan.

I've left a handful of comments, but actually I think there are some other things that we could do in this PR, that are harder to put in comments:

  • we should add a SteadyStateError of a diagnostic to another test (maybe Williamson 2)
  • we should add the PotentialEnergy diagnostic to an example to make sure that it's tested. Maybe the compressible DCMIP test?
  • while you've fixed the SteadyStateError of a diagnostic, we haven't fixed the Perturbation of a diagnostic. I think we should include a fix to that in this PR too.
  • a minor thing, but maybe you could update the name of the PR? I think it's more like "Fix diagnostics of diagnostics" -- you might have started by trying to making a minimum failing example but we don't want failing examples on main!

columns = 150
tmax = 3600.
dumpfreq = int(tmax / (2*dt))
nlayers = 5
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I don't know exactly what point you started this branch from, but I think that there should be almost no changes to this file


# Dumping point data using legacy PointDataOutput is not supported in parallel
if COMM_WORLD.size == 1:
output = OutputParameters(
dirname=dirname,
dumpfreq=dumpfreq,
pddumpfreq=dumpfreq,
dumplist=['u'],
point_data=[('theta_perturbation', points)],
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

please keep the point data outputting in here for now!

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

So i had no intention of pushing the change to this example and will move to testing through the examples suggested above. For the reference when including point data the examples wouldn't run so will most likely revert a lot of the changes to skamarock and move to testing in Williams 2 and dcmip etc

diagnostic_fields = [CourantNumber(), Gradient("u"), Perturbation('theta'),
Gradient("theta_perturbation"), Perturbation('rho'),
RichardsonNumber("theta", parameters.g/Tsurf), Gradient("theta")]
diagnostic_fields = [CourantNumber(), Pressure(eqns), SteadyStateError('Pressure_Vt')]
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think it makes sense to add something like Pressure(eqns) and SteadyStateError('Pressure_Vt') to a test. I don't think this is the right test as we don't expect this test to preserve a steady state.

I think it would make more sense to add SteadyStateError('relative_vorticity') [not sure if this is the right name so please check!] to the Williamson 2 test.

@@ -1200,7 +1221,10 @@ def setup(self, domain, state_fields):
state_fields (:class:`StateFields`): the model's field container.
"""
x = SpatialCoordinate(domain.mesh)
self.expr = self.rho_averaged * (1 + self.r_t) * self.parameters.g * dot(x, domain.k)
self._setup_thermodynamics(domain, state_fields)
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Thanks for fixing this!


self.init_field_set = False

#HACK Defo not the best way to do this perhaps we want a flag for checkpoint or prognostic
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

As discussed offline, I think we can avoid a hack like this by adding things to the else part of this statement:

field1 = state_fields(self.field_name1)
field2.assign(field1)

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

changed the fix to this by placing a check in the checkpointing routine to avoid it resetting as discuessed

@Witt-D Witt-D changed the title Diagnostic minimal fail Fix Diagnostics of Disagnosics Oct 11, 2023
Copy link
Contributor

@tommbendall tommbendall left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Thanks Dan! Good that the tests now pass.

A few more minor things to sort.

Do we also want to fix the Perturbation diagnostic as part of this PR?

PS: you might want to fix the typo in the name of the pull request!

examples/compressible/skamarock_klemp_nonhydrostatic.py Outdated Show resolved Hide resolved
examples/compressible/skamarock_klemp_nonhydrostatic.py Outdated Show resolved Hide resolved
examples/compressible/skamarock_klemp_nonhydrostatic.py Outdated Show resolved Hide resolved
examples/compressible/skamarock_klemp_nonhydrostatic.py Outdated Show resolved Hide resolved
examples/compressible/skamarock_klemp_nonhydrostatic.py Outdated Show resolved Hide resolved
@@ -1033,12 +1033,30 @@ def setup(self, domain, state_fields):
field1 = state_fields(self.field_name1)
field2 = state_fields(self.field_name2, space=field1.function_space(),
pick_up=True, dump=False)
# By default set this new field to the current value
# This may be overwritten if picking up from a checkpoint
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think we still want these deleted comments in here...
and this raises the point that we should make sure to comment our changes clearly so that we can still understand this change in 6 months time...

I suggest adding something like (you can use your own words!):

# The initial value for this field may not have been set yet, so we should postpone setting
# the initial field until the first time that the compute method is called


super().setup(domain, state_fields)

def compute(self):
if not self.init_field_set:
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Add comment here along the lines of:

# The first time the compute method is called, we should ensure the initial field is set
# However we don't want to do this if picking up from a checkpoint

@@ -630,6 +630,10 @@ def pick_up_from_checkpoint(self, state_fields):
self.output.checkpoint_pickup_filename = None
else:
raise ValueError("Must set checkpoint True if picking up")
# Prevent any steady-state diagnostics overwriting their original fields
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Very picky, but it would be good to have a blank line above this comment to help illustrate that the next segment of code is doing something separate from the above bits of code

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

So ive addressed the points regarding the example and some of the comments but was a bit uncertain on my explaining comments so they will be worth checking.

I think including perturbation in this will be a good idea too so will start having a dig around in that now

@Witt-D Witt-D changed the title Fix Diagnostics of Disagnosics Fix Diagnostics of Diagnostics Oct 12, 2023
Copy link
Contributor

@tommbendall tommbendall left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Some typos to fix

gusto/diagnostics.py Outdated Show resolved Hide resolved
gusto/diagnostics.py Outdated Show resolved Hide resolved
gusto/diagnostics.py Outdated Show resolved Hide resolved
@@ -216,6 +216,7 @@ def set_reference_profiles(self, reference_profiles):
set the reference field.
"""
for field_name, profile in reference_profiles:
breakpoint()
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

this is why you're breaking the tests

Copy link
Contributor

@tommbendall tommbendall left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This is very very nearly ready, you just need to fix the typos in a comment

gusto/timeloop.py Outdated Show resolved Hide resolved
Copy link
Contributor

@tommbendall tommbendall left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Thanks for doing this Dan! Good to go now

@tommbendall tommbendall merged commit 6a1678b into main Nov 13, 2023
4 checks passed
@tommbendall tommbendall deleted the diagnostic_minimal_fail branch November 13, 2023 16:51
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

Successfully merging this pull request may close these issues.

3 participants