miniufo / xinvert Goto Github PK
View Code? Open in Web Editor NEWInvert geophysical fluid dynamic problems (elliptic partial differential equations) using SOR iteration method.
Home Page: https://xinvert.readthedocs.io/
License: MIT License
Invert geophysical fluid dynamic problems (elliptic partial differential equations) using SOR iteration method.
Home Page: https://xinvert.readthedocs.io/
License: MIT License
This example is quite similar to the Gill-Matsuno model in the sense that the mass source/sink drives a horizontal circulation.
xinvert
uses many constants within the package. However, these constants need to be modified in some cases. For example, a numerical model may use 6371,000 as the radius of the earth, while here 6371,200 is used as default. We need to add a constants
module, so that these constants can be modified by users across the whole package, based on their needs.
Fofonoff (1954) proposed a nonlinear steady-state solution of an ocean, with an initial energy but no forcing and dissipation (conservative flow). This model fits the solver well.
To cite the packages windspharm and Dedalus, please use the citations to the software papers (rather than just linking to the package github repos):
This is part of the JOSS review process at openjournals/joss-reviews#5510.
Hello, I'm very new to python and trying to implement the QG omega test case using the sample data: atmos3D.nc provided, but I keep running into the following error upon executing the testOmegaEq.py script:
Traceback (most recent call last):
File "", line 1, in
File "", line 18, in
File "C:\Users\user-pc\miniconda3\Lib\site-packages\xinvert\tests\GeoApps\GridUtils.py", line 70, in add_latlon_metrics
dlonG = grid.diff(ds[lon ], 'X', boundary_discontinuity=360)
File "C:\Users\user-pc\miniconda3\lib\site-packages\xgcm\grid.py", line 2093, in diff
return self._1d_grid_ufunc_dispatch("diff", da, axis, **kwargs)
File "C:\Users\user-pc\miniconda3\lib\site-packages\xgcm\grid.py", line 1836, in _1d_grid_ufunc_dispatch
array = grid_ufunc(
File "C:\Users\user-pc\miniconda3\lib\site-packages\xgcm\grid_ufunc.py", line 460, in call
return apply_as_grid_ufunc(
File "C:\Users\user-pc\miniconda3\lib\site-packages\xgcm\grid_ufunc.py", line 770, in apply_as_grid_ufunc
results = _apply(
File "C:\Users\user-pc\miniconda3\lib\site-packages\xgcm\grid_ufunc.py", line 837, in _apply
results = xr.apply_ufunc(
TypeError: apply_ufunc() got an unexpected keyword argument 'boundary_discontinuity'
Can you please assist in figuring out what is causing this?
Hi,
I am trying to work on the 11_Omega_equation.ipynb notebook. I am currently experiencing "Exception: coordinate lev is non-uniform". The dataset I am using has the following pressure levels array([100000., 97500., 95000., 92500., 90000., 87500., 85000., 82500., 80000., 77500., 75000., 72500., 70000., 65000., 60000., 55000.,50000., 45000., 40000., 35000., 30000., 25000., 20000.])
Do I need to make those levels uniform or is there any other way I could work on this issue without changing the current levels?
Recently, I find that the steady state for 2D shallow-water model can be cast as a 2nd-order PDE that fits the solver. So just do it.
As I want to fit more problems into the architecture of xinvert
, the current design of this package makes it a bit hard for some problems. Specifically, it is much better to refactor the internal design of this package, and also change the user APIs accordingly. Specifically, the problems are:
mParams
that prescribed in the model or equation. For example, Coriolis paramter iParams
that controls the iteration and has nothing to do with the model. For example, maximum loops or tolerance for stop looping.xgcm
. However, xgcm
is still in heavy development and is targeted at Arakawa C-grid. Since I don't want this package to be fully relied on xgcm
, I need to add this finite difference capability.animate_iteration
function fo all the models.The Q-vector form of the Omega equation doesn't have an adiabatic forcing term in your calculation. Is there any easy way to incorporate that term while calculating the Omega value from the Qvector form of the Omega equation?
This is an example I recently ran into, which adopted a variational calculas to get the steady-state equation between a 2D flow with total kinetic energy
Once the topography is specified, the steady state of the flow with
Hello,
I'm working with the dynamic pressure perturbation equation (Laplacian of p = -roh*((dudx)**2+(dvdy)**2+(dwdz)**2)). At small scales, what boundary conditions, mxLoop, and tolerance would you reccomend?
A summary table of all the current problems would be good for the users to see the functionality of xinvert
. Also, such a table makes an example for the users and they can also add their own problems.
Hi, I've run xinvert for the Gill-Matsuno problem over the maritime continent, as shown in your example. The results I've got regarding h1, u1 and v1 are constrained just to the heating source region (see attached figure), but not to the entire Indo-Pacific region as in your figures. I think the problem can be in the plotting of the wind and streamfunction, which I include below. I appreciate any help. Many thanks in advance.
##########################
import os
import matplotlib.pyplot as plt
from netCDF4 import Dataset as netcdf_dataset
import numpy as np
from cartopy import config
import cartopy.crs as ccrs
import cartopy.feature as cfeature
from cartopy.mpl.ticker import (LongitudeFormatter, LatitudeFormatter,
LatitudeLocator)
plt.figure(figsize=(40,20))
heat = Q1
lats = lat
lons = lon
skip = 1
fontsize = 16
ax = plt.axes(projection=ccrs.PlateCarree())
ax.coastlines()
ax.add_feature(cfeature.BORDERS)
ax.gridlines(crs=ccrs.PlateCarree(), draw_labels=True,
xlocs=np.arange(-140,100, 80), ylocs=np.arange(-40, 40, 40),
x_inline=False, y_inline=False, linewidth=0.33, color='k',alpha=0.5)
cmap = plt.colormaps['jet']
ax.contourf(lons, lats, Q1, 10, transform=ccrs.PlateCarree(), cmap=cmap , levels=np.linspace(-0.05, 0.05, 10))
ax.contour(lons, lats, h1, transform=ccrs.PlateCarree(), cmap='jet')
windspeed = (u1 ** 2 + v1 ** 2) ** 0.5
windspeed.rename('windspeed')
ax.contour(lons, lats, windspeed, 20, cmap='jet')
plt.quiver(lons, lats, u1, v1, pivot='middle', width=0.0015, headwidth=1., headlength=1., minlength=1)
#Limits for Maritime Continent
ax.set_ylim([-40, 40])
ax.set_xlim([40, 200])
plt.show()
Hi,
Is it possible to pass different values of N2 per pressure level, and f0 per latitude, in the form of coordinate-aware arrays? In particular, N2 is quite different for the upper and lower levels.
Thanks!
Vertical mode decomposition (VMD) is a common way to analyze the vertical structure of many phenomena. Generally, given a vertical profiles of
Hi, I am trying to reproduce the vertical velocity using QG Omega inversion. If I need to invert the equation for the western hemisphere only, how do I need to change the boundary conditions?
@miniufo - well done on your submission to JOSS.
There is a requirement on the reviewer checklist relating to the software paper which reads as follows:
State of the field: Do the authors describe how this software compares to other commonly-used packages?
(See openjournals/joss-reviews#5510 for details.)
You make a comparison against windspharm in your README file...
Inversion on the spherical earth, like some meteorological problems, could utilize the spherical harmonics like windspharm, which would be more efficient using FFT than SOR used here. However, in the case of ocean, SOR method is definitely a better choice in the presence of irregular land/sea mask.
... but not in the software paper itself. I think it would be worth mentioning windspharm (and any other relevant commonly-used packages) in your software paper too.
@miniufo - well done on your submission to JOSS.
As part of my review of your submission (see openjournals/joss-reviews#5510), I have some questions/suggestions regarding software installation:
Hi, this package looks really neat!
I was wondering if you had a reference for how you are calculating the finite differences with spherical coordinates? From the looks of it, you’re simply multiplying the Y-axis by a scalar factor of degrees to meters and the X-axis by the same plus and a cosine of the Y-coords. I’m not familiar with this but this seems really, really useful in general with finite differences for lat-lon.
Do you have a reference or something I could look at because I’m not familiar with this but I would really like to use it.
Thanks!
When trying to re-run the notebooks, three of them fail:
11_Omega_equation.ipynb
the following fails because the netcdf file is not availableds = xr.open_dataset('I:/Omega/OFES30_20011206_Qian/data.nc',
chunks={'lev':6}).astype('f4')
04_Eliassen_model.ipynb
the following fails:ds = xr.open_dataset('../../../Data/Zonalmean.nc')
The latter is simply a typo in the filename; it needs to be changed to ./../../Data/ZonalMean.nc
test_01_Stommel_Arons_model.ipynb
the following fails:eps = xr.where(basin, eps, np.nan)
with the error
ValueError: cannot align objects with join='exact' where index/labels/sizes are not equal along these coordinates (dimensions): 'lon' ('lon',)
This is part of the review process at openjournals/joss-reviews#5510.
@miniufo - well done on your submission to JOSS.
As it stands, your package doesn't meet the following requirement on the reviewer checklist:
Automated tests: Are there automated tests or manual steps described so that the functionality of the software can be verified?
(See openjournals/joss-reviews#5510 for details.)
You've got some test code available at tests/
and associated data files at Data/
, but there is no guidance provided on how to run those tests manually or better still there is no automated testing.
In addition, it looks like your test code culminates in the plotting of figures. While that can be useful for visual verification, common practice would also include a suite of tests that have the following components:
The comparison between the actual result and expected result would take the form of an assertion, so that the test can be said to pass or fail. You can then use a test runner like pytest in conjunction with GitHub actions to setup automated testing.
You may already be familiar with these testing concepts, but if not let me know and I can point you to some resources that explain how to setup automated testing.
Hi,
I'm trying to use the method in a domain which has 19 levels in depth (totally 1000m), and I found it makes a poor inversion result. However, it comes a much better result after I set my vertical levels as the data used in oceanic example (200 levels for 1000m), which convinced me that 'dz' counts a lot in calculation. But with the improvement comes a fat data array ( [360,200,1700,650] ) that is so hard to run forward. So I’d like to know, from your perspective, what value of 'dz' (or levels in depth) can allows me to reduce the size of data while maintaining a good accuracy?
Thanks a lot!
It would be helpful to have the installation instructions not only in the README, but also in the documentation.
This is part of the JOSS review process at openjournals/joss-reviews#5510.
Although the current implementation is based on xarray
and dask
. The solution is first loaded completely into memory. During the multi timesteps calculation, xinvert
does not accelerate on a multi-core PC. This is a little strange to me and need some times to make it parallelable.
Is there any DOI to cite this code?
@miniufo - well done on your submission to JOSS.
As it stands, your package doesn't meet the following requirement on the reviewer checklist:
Community guidelines: Are there clear guidelines for third parties wishing to 1) Contribute to the software 2) Report issues or problems with the software 3) Seek support
(See openjournals/joss-reviews#5510 for details.)
I'm opening a new issue here to continue the discussion on the issue of the xcontour dependence that we started here.
The dependence of xcontour is a bit complicated. It is another package I've written here. So the example of reference state, which is newly added, depends on xcontour to prepare the forcing (actually preparing the forcing data is not the core part of xinvert). This is not necessary for other examples. So far, I don't have any plan to get xcontour released as a conda package, so it is not possible to add this to conda dependency. Do you have any suggestion? I really like this example to remain in the problem list.
Is there a way to precompute M(Q) and C(Q) with xcontour
and in the https://xinvert.readthedocs.io/en/latest/notebooks/05_reference_SWM.html example, and just load these states from netcdf if xcontour
is not available? You could then adapt the notebook as follows:
try:
import xcontour
# do xcontour computation in this notebook, as you do now
except ImportError:
ds = xr.open_dataset(your_precomputed_dataset)
This is part of the review process at openjournals/joss-reviews#5510.
Currently, the boundary conditions are specified as:
sf = invert_Poisson(vor, dims=['lat', 'lon'], BCs=['fixed', 'periodic'])
which means the north-south BCs are all fixed and east-west BCs are periodic. This is not flexible enough if one wants the north BC as the 'fixed' one while the south BC as the 'extend' or 'reflect' one.
One possible solution is to use dict to specify BCs as:
sf = invert_Poisson(vor, dims=['lat', 'lon'], BCs={'lat':('fixed'. 'extend'), 'lon':('periodic', 'periodic')})
Also, I am thinking if the underlying algorithm should be refactored as two steps: 1) pad the domain with proper BCs and 2) iterate for the solution. But this may cause some problems if BCs are at 1) poles and 2) vortex center where coordinates are undefined.
It would be convenient to submit a JOSS paper for later citation.
A declarative, efficient, and flexible JavaScript library for building user interfaces.
🖖 Vue.js is a progressive, incrementally-adoptable JavaScript framework for building UI on the web.
TypeScript is a superset of JavaScript that compiles to clean JavaScript output.
An Open Source Machine Learning Framework for Everyone
The Web framework for perfectionists with deadlines.
A PHP framework for web artisans
Bring data to life with SVG, Canvas and HTML. 📊📈🎉
JavaScript (JS) is a lightweight interpreted programming language with first-class functions.
Some thing interesting about web. New door for the world.
A server is a program made to process requests and deliver data to clients.
Machine learning is a way of modeling and interpreting data that allows a piece of software to respond intelligently.
Some thing interesting about visualization, use data art
Some thing interesting about game, make everyone happy.
We are working to build community through open source technology. NB: members must have two-factor auth.
Open source projects and samples from Microsoft.
Google ❤️ Open Source for everyone.
Alibaba Open Source for everyone
Data-Driven Documents codes.
China tencent open source team.