[mpas-developers] Vector Reconstruction
Xylar S. Asay-Davis
xylar at lanl.gov
Wed Mar 24 16:44:31 MDT 2010
Hi Michael,
Great, I'm glad you've figured out more exactly what the problem is. I
would suggest that we modify the reconstruction code as you suggest below.
There's no reason that I can see to restrict the edgesOnCell to be
ordered the same way as the cellsOnCell, (though I guess I would have
expected this to happen by default in grid generation).
Please feel free to add the reconstruction code to the repository once
you're happy with it.
-Xylar
> Hi, Todd.
>
> You're right: removing -D_MPI should work to run the code in
> serial. I'd agree that moving this to a definition at the top of
> the Makefile would be a good way to simplify the work of compiling
> the code for serial execution.
>
> Regarding the differences in behavior of the vector reconstruction
> when changing the grid.nc file, I think I may have found the
> issue. The vector reconstruction code makes the assumption that
> edgesOnCell(i,iCell) is the edge separating iCell from
> cellsOnCell(i,iCell), i=1,nEdgesOnCell(iCell), and this assumption
> doesn't appear to hold for the grid.nc file that I've been testing
> with (based on some debugging code I added to the vector
> reconstruction module). One solution would be to rewrite the code
>
> do j=1,nEdgesOnCell(iCell)
> cell2 = cellsOnCell(j,iCell)
> iEdge = edgesOnCell(j,iCell)
> X2(1) = xCell(cell2)
> X2(2) = yCell(cell2)
> X2(3) = zCell(cell2)
> if (iCell == cellsOnEdge(1,iEdge)) then
> normals(:,j) = X2(:) - X1(:)
> else
> normals(:,j) = X1(:) - X2(:)
> endif
>
> call unit_vector_in_R3(normals(:,j))
>
> xLoc(:,j,iCell) = 0.5*(X2(:) + X1(:))
> enddo
>
> to
>
> do j=1,nEdgesOnCell(iCell)
> iEdge = edgesOnCell(j,iCell)
> if (iCell == cellsOnEdge(1,iEdge)) then
> cell2 = cellsOnEdge(2,iEdge)
> X2(1) = xCell(cell2)
> X2(2) = yCell(cell2)
> X2(3) = zCell(cell2)
> normals(:,j) = X2(:) - X1(:)
> else
> cell2 = cellsOnEdge(1,iEdge)
> X2(1) = xCell(cell2)
> X2(2) = yCell(cell2)
> X2(3) = zCell(cell2)
> normals(:,j) = X1(:) - X2(:)
> endif
>
> call unit_vector_in_R3(normals(:,j))
>
> xLoc(:,j,iCell) = 0.5*(X2(:) + X1(:))
> enddo
>
>
> to remove the assumption; doing this, I get reasonable values for
> uReconstructZ for SWTC2. Another solution might be to just require
> that grid.nc files provide edgesOnCell and cellsOnCell arrays that
> begin at the same angle and list edges/cells in counter-clockwise
> order, satisfying the assumption currently made in the vector
> reconstruction code; it sounds like our current repository
> versions of the grid generators might already satisfy this, but it
> would be something that we'd need to specify in case grids were to
> be generated through alternative codes. Given that the fix in the
> vector reconstruction code is relatively small, I'd lean in favor
> of imposing fewer requirements on grid.nc files at this point, but
> we might find that, down the road a bit, having the luxury of
> making this assumption would make our lives much easier. What do
> you all think?
>
> Cheers!
> Michael
>
>
> On Wed, Mar 24, 2010 at 04:00:43PM -0600, Todd Ringler wrote:
>>
>> Michael,
>>
>> When I want to run in serial mode I just remove the "-D_MPI" compile
>> flag under the trunk code. I am not sure if this is suppose to work,
>> but it does (at least for the times I have tried it).
>>
>> I recompiled the trunk sw for serial model and I get errors in z
>> velocity of about 8e-3.
>>
>> (If we think that the trunk code should work in serial mode when D_MPI
>> is omitted, maybe at some point we can move that option to the top of
>> the Makefile to allow folks to commit it out.)
>>
>> Cheers,
>> Todd
>>
>>
>> On Mar 24, 2010, at 3:33 PM, Michael Duda wrote:
>>
>> >Hi, Todd.
>> >
>> >This is definitely very odd. The executable is named swmodel
>> >rather than sw_model.exe since the code is from the swmodel branch
>> >(formerly in the trunk), which doesn't support multiple cores
>> >under the same directory structure. In our current trunk, the
>> >executable is always named $(CORE)_model.exe, whereas we had
>> >simply called it swmodel in the previous trunk/swmodel code.
>> >
>> >The grid is the quasi-uniform 10242 mesh that I've had posted at
>> >http://www.mmm.ucar.edu/people/duda/files/mpas, and was created
>> >based on the generating points you had sent sometime late last
>> >year. I don't think the problem is connected with the mesh, since
>> >the 'u' and 'v' fields (i.e., the normal and tangential velocity
>> >fields) in the output.nc file look correct when I plot them with
>> >ncl/cells.ncl. The height field from the simulation remains
>> >stable, so I don't think there's a problem with the underlying
>> >simulation.
>> >
>> >If you or Xylar could place a tar file of the code, either
>> >trunk/mpas or branches/swmodel, that you've used in your tests
>> >with the new module_vector_reconstruction.F, I could peruse that
>> >and see whether I can spot any key differences in how the module
>> >is used. Just let me know if there's any testing that I can do to
>> >help out here.
>> >
>> >Cheers,
>> >Michael
>> >
>> >
>> >On Wed, Mar 24, 2010 at 03:24:58PM -0600, Todd Ringler wrote:
>> >>
>> >>Hi Michael,
>> >>
>> >>I also reproduced your results.
>> >>
>> >>How can this be? We must be missing something obvious here.
>> >>
>> >>I notice that the executable is named "swmodel" and sw_model.exe --
>> >>was the reconstruction tested on top of the trunk? Also, is this the
>> >>quasi uniform 10242 mesh? built with the current version of the grid
>> >>generator?
>> >>
>> >>Cheers,
>> >>Todd
>> >>
>> >>On Mar 24, 2010, at 3:00 PM, Xylar S. Asay-Davis wrote:
>> >>
>> >>>Hi Michael,
>> >>>
>> >>>I'm able to reproduce your results. It doesn't look good! I'll see
>> >>>what
>> >>>I can figure out.
>> >>>
>> >>>-Xylar
>> >>>
>> >>>
>> >>>>Hi, Todd.
>> >>>>
>> >>>>I wonder whether there's something odd with how I've added the new
>> >>>>module_vector_reconstruction.F code and Registry entry. I had been
>> >>>>testing with branches/swmodel for simplicity, and the changes are
>> >>>>really trivial. If you wouldn't mind, could you give a try with
>> >>>>the tar file I've been testing with? I've placed this at
>> >>>>http://www.mmm.ucar.edu/people/duda/files/swmodel_vectors.tar.gz.
>> >>>>The namelist is set up to run just 5 time steps, and the input
>> >>>>files are all included, so it should just be a matter of compiling
>> >>>>and running. Also included in the tar file is the output file that
>> >>>>I generated for SWTC2; doing an ncdump or using ncview, you can
>> >>>>see the large values in uReconstructZ.
>> >>>>
>> >>>>If you're able to get correct results with my tar file, I'd not
>> >>>>object to committing the new module_vector_reconstruction.F; I can
>> >>>>propose to add this in a new shared directory 'operators', to be
>> >>>>created in trunk/mpas/src with appropriate changes to the
>> >>>>Makefiles. However, if you can reproduce the large values of
>> >>>>uReconstructZ with my tar file, I can look for any differences if
>> >>>>you could supply the code that you were able to successfully run
>> >>>>with. Does that sound like a reasonable way forward?
>> >>>>
>> >>>>Cheers,
>> >>>>Michael
>> >>>>
>> >>>>
>> >>>>On Wed, Mar 24, 2010 at 12:06:19PM -0600, Todd Ringler wrote:
>> >>>>>
>> >>>>>Hi Michael,
>> >>>>>
>> >>>>>I can not reproduce your finding. Below is the rbf reconstruction
>> >>>>>(vectors) and z-component of velocity (colors with colorbar, units
>> >>>>>m/
>> >>>>>s). This is SWTC#2 on the quasi-uniform 2562 grid. I am happy to
>> >>>>>followup over the phone. Have to run, back at 2 pm.
>> >>>>>
>> >>>>>Cheers,
>> >>>>>Todd
>> >>>>>
>> >>>>
>> >>>>
>> >>>>
>> >>>>>
>> >>>>>On Mar 23, 2010, at 11:56 AM, Michael Duda wrote:
>> >>>>>
>> >>>>>>Hi, Xylar.
>> >>>>>>
>> >>>>>>Thanks for the quick reply. In SWTC2, the maximum zonal velocity
>> >>>>>>is around 38 m/s at the equator, so 10 m/s for uReconstructZ does
>> >>>>>>seem out of place. Just to make sure I understand how to
>> >>>>>>interpret
>> >>>>>>uReconstruct{X,Y,Z}, these are the components of the velocity
>> >>>>>>vector in Cartesian space, with the z-axis positive through the
>> >>>>>>north pole of the sphere, right? I had been testing on a fairly
>> >>>>>>coarse grid (10242 cells), so that could be a factor; I'll do a
>> >>>>>>little more testing at higher resolution.
>> >>>>>>
>> >>>>>>Cheers,
>> >>>>>>Michael
>> >>>>>>
>> >>>>>>
>> >>>>>>On Tue, Mar 23, 2010 at 11:41:18AM -0600, Xylar S. Asay-Davis
>> >>>>>>wrote:
>> >>>>>>>Hi Michael,
>> >>>>>>>
>> >>>>>>>I don't think your attachment made it through.
>> >>>>>>>
>> >>>>>>>I'll take a look at shallow water test case 2 to see what
>> >>>>>>>might be
>> >>>>>>>the
>> >>>>>>>problem there. The value of uReconstructZ won't machine zero
>> >>>>>>>for a
>> >>>>>>>purely
>> >>>>>>>zonal flow because the reconstruction in its current form will
>> >>>>>>>only
>> >>>>>>>guarantee exact reconstruct of a constant flow in the tangent
>> >>>>>>>plane
>> >>>>>>>of the
>> >>>>>>>cell, where as a zonal flow has components slightly out of the
>> >>>>>>>tangent
>> >>>>>>>plane. In my tests, I was seeing uReconstructZ smaller than
>> >>>>>>>uReconstructX
>> >>>>>>>and uReconstructY by a factor of around 1e-6.
>> >>>>>>>
>> >>>>>>>-Xylar
>> >>>>>>>
>> >>>>>>>
>> >>>>>>>>Hi, Xylar.
>> >>>>>>>>
>> >>>>>>>>The module looks good to me -- very clean code. The only
>> >>>>>>>>changes
>> >>>>>>>>I'd suggest are to the comments at the beginning of the
>> >>>>>>>>init_reconstruct() routine, and a few changes to indentation,
>> >>>>>>>>as
>> >>>>>>>>in the attached module_vector_reconstruction.F.
>> >>>>>>>>
>> >>>>>>>>This is most likely a problem in how I'm performing the
>> >>>>>>>>computation given that both you and Todd have checked the
>> >>>>>>>>reconstructed vectors, but I'm having trouble getting a zonal
>> >>>>>>>>and
>> >>>>>>>>meridional wind from the reconstructed vector in (x,y,z). Even
>> >>>>>>>>plotting the magnitude of the reconstructed vector, sqrt(x^2 +
>> >>>>>>>>y^2
>> >>>>>>>>+ z^2), I'm getting some unexpected values in the reconstructed
>> >>>>>>>>wind speed field for shallow water test case 2 (with alpha=0).
>> >>>>>>>>Given that we should have purely zonal flow for this test case,
>> >>>>>>>>I'm suprised to see values in uReconstructZ of order 10. Do
>> >>>>>>>>either
>> >>>>>>>>you or Todd have any thoughts or suggestions?
>> >>>>>>>>
>> >>>>>>>>Cheers,
>> >>>>>>>>Michael
>> >>>>>>>>
>> >>>>>>>>
>> >>>>>>>>On Mon, Mar 22, 2010 at 03:30:39PM -0600, Xylar S. Asay-Davis
>> >>>>>>>>wrote:
>> >>>>>>>>>Attached is my proposed vector reconstruction code. I have
>> >>>>>>>>>run
>> >>>>>>>>>several
>> >>>>>>>>>sanity checks and it seems to work fine. However, I haven't
>> >>>>>>>>>been
>> >>>>>>>>>able
>> >>>>>>>>>to
>> >>>>>>>>>get a full ocean core run with this code. It is not clear
>> >>>>>>>>>to me
>> >>>>>>>>>if this
>> >>>>>>>>>has to do with this code or with recent changes that Mark has
>> >>>>>>>>>made to
>> >>>>>>>>>the
>> >>>>>>>>>ocean core relating to 3D density. I'll keep working with
>> >>>>>>>>>Mark
>> >>>>>>>>>to try
>> >>>>>>>>>to
>> >>>>>>>>>figure that out, but in the mean time it would be great to
>> >>>>>>>>>have
>> >>>>>>>>>Todd
>> >>>>>>>>>and/or Michael test out this code to see if it will work for
>> >>>>>>>>>you.
>> >>>>>>>>>
>> >>>>>>>>>A registry entry for the coeffs_reconst is required as follows
>> >>>>>>>>>(already
>> >>>>>>>>>in
>> >>>>>>>>>the ocean core):
>> >>>>>>>>>var real coeffs_reconstruct ( R3 maxEdges nCells ) -
>> >>>>>>>>>coeffs_reconstruct
>> >>>>>>>>>- -
>> >>>>>>>>>
>> >>>>>>>>>Registry entries for matrix_reconstruct, normal and
>> >>>>>>>>>rbf_value
>> >>>>>>>>>are no
>> >>>>>>>>>longer required by this code.
>> >>>>>>>>>
>> >>>>>>>>>-Xylar
>> >>>>>>>>
>> >>>>>>>>>_______________________________________________
>> >>>>>>>>>mpas-developers mailing list
>> >>>>>>>>>mpas-developers at mailman.ucar.edu
>> >>>>>>>>>http://mailman.ucar.edu/mailman/listinfo/mpas-developers
>> >>>>>>>>
>> >>>>>>>>_______________________________________________
>> >>>>>>>>mpas-developers mailing list
>> >>>>>>>>mpas-developers at mailman.ucar.edu
>> >>>>>>>>http://mailman.ucar.edu/mailman/listinfo/mpas-developers
>> >>>>>>>>
>> >>>>>>>
>> >>>>>>>
>> >>>>>>>***********************
>> >>>>>>>Xylar S. Asay-Davis
>> >>>>>>>E-mail: xylar at lanl.gov
>> >>>>>>>Phone: (505) 606-0025
>> >>>>>>>Fax: (505) 665-2659
>> >>>>>>>CNLS, MS B258
>> >>>>>>>Los Alamos National Laboratory
>> >>>>>>>Los Alamos, NM 87545
>> >>>>>>>***********************
>> >>>>>>>
>> >>>>>>_______________________________________________
>> >>>>>>mpas-developers mailing list
>> >>>>>>mpas-developers at mailman.ucar.edu
>> >>>>>>http://mailman.ucar.edu/mailman/listinfo/mpas-developers
>> >>>>>
>> >>>>
>> >>>>_______________________________________________
>> >>>>mpas-developers mailing list
>> >>>>mpas-developers at mailman.ucar.edu
>> >>>>http://mailman.ucar.edu/mailman/listinfo/mpas-developers
>> >>>>
>> >>>
>> >>>
>> >>>***********************
>> >>>Xylar S. Asay-Davis
>> >>>E-mail: xylar at lanl.gov
>> >>>Phone: (505) 606-0025
>> >>>Fax: (505) 665-2659
>> >>>CNLS, MS B258
>> >>>Los Alamos National Laboratory
>> >>>Los Alamos, NM 87545
>> >>>***********************
>> >>>
>> >>
>>
> _______________________________________________
> mpas-developers mailing list
> mpas-developers at mailman.ucar.edu
> http://mailman.ucar.edu/mailman/listinfo/mpas-developers
>
***********************
Xylar S. Asay-Davis
E-mail: xylar at lanl.gov
Phone: (505) 606-0025
Fax: (505) 665-2659
CNLS, MS B258
Los Alamos National Laboratory
Los Alamos, NM 87545
***********************
More information about the mpas-developers
mailing list