ClimateGlobalChange / tempestremap

TempestRemap: Remapping software for climate applications
Other
39 stars 28 forks source link

FV->FV bilin and related algorithms dont work with MPAS oEC60to30v3 source grid #107

Open mt5555 opened 1 year ago

mt5555 commented 1 year ago

error with intbilin. But bilin and intbilinbg both work.

Parameters:
  --in_mesh <string> ["/Users/mt/scratch1/mapping/grids/TEMPEST_ne30pg2.g"]
  --out_mesh <string> ["/Users/mt/scratch1/mapping/grids/ocean.oEC60to30v3.scrip.181106.nc"]
  --ov_mesh <string> ["/Users/mt/scratch1/mapping/maps/overlap_oEC60to30v3_ne30pg2.g"]
  --in_type <string> ["fv"] [fv|cgll|dgll]
  --out_type <string> ["fv"] [fv|cgll|dgll]
  --out_map <string> ["/Users/mt/scratch1/mapping/maps/map_ne30pg2_to_oEC60to30v3_intbilin.nc"]
  --in_meta <string> [""]
  --out_meta <string> [""]
  --in_concave <bool> [false]
  --out_concave <bool> [false]
  --in_np <integer> [1]
  --out_np <integer> [1]
  --method <string> ["intbilin"]
  --mono <bool> [false]
  --nobubble <bool> [false]
  --correct_areas <bool> [false]
  --nocorrectareas <bool> [false]
  --noconserve <bool> [true]
  --nocheck <bool> [false]
  --sparse_constraints <bool> [false]
  --volumetric <bool> [false]
  --mono2 <bool> [false]
  --mono3 <bool> [false]
  --in_data <string> [""]
  --out_data <string> [""]
  --var <string> [""]
  --ncol_name <string> ["ncol"]
  --out_double <bool> [true]
  --preserve <string> [""]
  --preserveall <bool> [false]
  --fillvalue <double> [0.000000]
  --out_format <string> ["Netcdf4"] [Classic|Offset64Bits|Netcdf4|Netcdf4Classic]
------------------------------------------------------------
Initializing dimensions of map
..Input mesh
..Output mesh
Loading input mesh
..Mesh size: Nodes [21602] Elements [21600]
Loading output mesh
..SCRIP Format File detected
..1167285 duplicate nodes detected
..Mesh size: Nodes [478835] Elements [235160]
Loading overlap mesh
..Mesh size: Nodes [632018] Elements [386577]
..Calculating input mesh Face areas
....Input Mesh Geometric Area: 1.256637061435917e+01 (1.000000000000000e+00)
..Calculating output mesh Face areas
....Output Mesh Geometric Area: 8.873230289717204e+00 (7.061092309006118e-01)
..Overlap mesh reverse correspondence found (reversing)
Calculating overlap mesh Face areas
..Overlap Mesh Area: 8.873230289716783e+00 (7.061092309005783e-01)
Correcting source/target areas to overlap mesh areas
WARNING: Significant mismatch between overlap mesh area and input mesh area.
  Disabling checks for consistency.
Mesh size: Edges [43200]
Calculating offline map (intbilin)
..Mesh size: Edges [43200]
..Element 0/21600
..EXCEPTION (src/LinearRemapFV.cpp, Line 2523) Non-monotone weight
mt5555 commented 1 year ago

for ocn->atm maps, bilin, intbilin and intbilinbg all fail:

Parameters:
  --in_mesh <string> ["/Users/mt/scratch1/mapping/grids/ocean.oEC60to30v3.scrip.181106.nc"]
  --out_mesh <string> ["/Users/mt/scratch1/mapping/grids/TEMPEST_ne30pg2.scrip.nc"]
  --ov_mesh <string> ["/Users/mt/scratch1/mapping/maps/overlap_oEC60to30v3_ne30pg2.g"]
  --in_type <string> ["fv"] [fv|cgll|dgll]
  --out_type <string> ["fv"] [fv|cgll|dgll]
  --out_map <string> ["/Users/mt/scratch1/mapping/maps/map_oEC60to30v3_to_ne30pg2_bilin.nc"]
  --in_meta <string> [""]
  --out_meta <string> [""]
  --in_concave <bool> [false]
  --out_concave <bool> [false]
  --in_np <integer> [1]
  --out_np <integer> [1]
  --method <string> ["bilin"]
  --mono <bool> [false]
  --nobubble <bool> [false]
  --correct_areas <bool> [false]
  --nocorrectareas <bool> [false]
  --noconserve <bool> [true]
  --nocheck <bool> [false]
  --sparse_constraints <bool> [false]
  --volumetric <bool> [false]
  --mono2 <bool> [false]
  --mono3 <bool> [false]
  --in_data <string> [""]
  --out_data <string> [""]
  --var <string> [""]
  --ncol_name <string> ["ncol"]
  --out_double <bool> [true]
  --preserve <string> [""]
  --preserveall <bool> [false]
  --fillvalue <double> [0.000000]
  --out_format <string> ["Netcdf4"] [Classic|Offset64Bits|Netcdf4|Netcdf4Classic]
------------------------------------------------------------
Initializing dimensions of map
..Input mesh
..Output mesh
Loading input mesh
..SCRIP Format File detected
..1167285 duplicate nodes detected
..Mesh size: Nodes [478835] Elements [235160]
Loading output mesh
..SCRIP Format File detected
..64798 duplicate nodes detected
..Mesh size: Nodes [21602] Elements [21600]
Loading overlap mesh
..Mesh size: Nodes [632018] Elements [386577]
..Calculating input mesh Face areas
....Input Mesh Geometric Area: 8.873230289717204e+00 (7.061092309006118e-01)
..Calculating output mesh Face areas
....Output Mesh Geometric Area: 1.256637061435917e+01 (1.000000000000000e+00)
..Overlap mesh forward correspondence found
Calculating overlap mesh Face areas
..Overlap Mesh Area: 8.873230289716785e+00 (7.061092309005785e-01)
Correcting source/target areas to overlap mesh areas
Mesh size: Edges [714274]
Calculating offline map (bilin)
..EXCEPTION (./src/GridElements.h, Line 509) FacePair already has a full set of Faces.
dmarsico1 commented 1 year ago

Regarding the first issue, it appears to have been an issue with the tolerance for the iterative method used to solve for the weights int "intbilin." If you go to line 933 of "NewtonQuadrilateral" FiniteVolumeTools.cpp, and change the stopping criterion to 1e-30, then it works fine.

Regarding the second issue, it is related to the ConstructEdgeMap() function. The methods "bilin," "intbilin," "bilingb," and "intbilingb" all rely on the dual of the source mesh and its edge map. There appears to be no issue constructing the dual of the ocean mesh, but the error "FacePair already has a full set of Faces" occurs with ConstructEdgeMap(), i.e. when the edge map is being constructed. I'm not sure why there would be a problem, however. Do you notice any problems with the dual of the ocean grid?

mt5555 commented 1 year ago

Thanks. The 1e-30 tol for atm->ocn maps fixed the issue.

For ocn->atm maps, when you wrote"Did you notice any problems with the dual of the ocean grid?" - what did you mean? My understanding is that this is computed internally in TR that I dont have access to? Note that these are MPAS ocean grids, which only mesh the ocean so they are not global.

dmarsico1 commented 1 year ago

The first two steps of the schemes that are currently failing are to construct the dual of the source mesh and its edgemap. Something isn't working in one of these steps. Is there any reason why a mesh that is not global would cause either the function "Dual" (GridElements.cpp line 3128) or "ConstructEdgeMap" to fail? I've only tested these schemes on global meshes.

There is nothing wrong with the underlying algorithms here; they will work as long as the dual of the source mesh and its edgemap are constructed properly.

dmarsico1 commented 1 year ago

There's a new pull request that should fix the problems for FV->FV bilin for MPAS oEC60to30v3. Also, it should work for the function LinearRemapGeneralizedBarycentric, but there doesn't seem to be a command line parameter to use this function. I haven't implemented the changes for intbilin or intbilinbg yet.

rljacob commented 1 year ago

Is this fixed by PR #112 ?

dmarsico1 commented 1 year ago

Yes, this is fixed by PR #112 (at least for bilin and bilingb. I'm working on the intbilin and intbilingb).