Skip to content

Commit

Permalink
Fix typos in FaceDivFree interior interpolation. (#4048)
Browse files Browse the repository at this point in the history 
The problem was masked in test in amrex/Tests/DivFreePatch because all
three components were the same.

## Summary

## Additional background

## Checklist

The proposed changes:
- [ x] fix a bug or incorrect behavior in AMReX
- [ ] add new capabilities to AMReX
- [x] changes answers in the test suite to more than roundoff level
- [ ] are likely to significantly affect the results of downstream AMReX
users
- [ ] include documentation in the code and/or rst files, if appropriate
  • Loading branch information
@cgilet
cgilet authored Jul 30, 2024
1 parent 20e6f2e commit a8c5255
Show file tree
Hide file tree
Showing 3 changed files with 40 additions and 21 deletions.
16 changes: 8 additions & 8 deletions Src/AmrCore/AMReX_Interp_3D_C.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -294,13 +294,13 @@ facediv_int (int ci, int cj, int ck, int nf,
+ dy*(2*dx*dx+dy*dy)/(8*dz*xspys)*(w000+w012-w002-w010)
+ dy3/(8*dz*xspys)*(w100+w112-w102-w110);

fine[1](fi+1, fj+1, fk , nf) = Real(0.5)*(v001+v021)
fine[1](fi+1, fj+1, fk , nf) = Real(0.5)*(v100+v120)
+ dy*(2*dz*dz+dy*dy)/(8*dx*yspzs)*(u000+u210-u010-u200)
+ dy3/(8*dx*yspzs)*(u001+u211-u011-u201)
+ dy*(2*dx*dx+dy*dy)/(8*dz*xspys)*(w100+w112-w102-w110)
+ dy3/(8*dz*xspys)*(w000+w012-w002-w010);

fine[1](fi , fj+1, fk+1, nf) = Real(0.5)*(v100+v120)
fine[1](fi , fj+1, fk+1, nf) = Real(0.5)*(v001+v021)
+ dy*(2*dz*dz+dy*dy)/(8*dx*yspzs)*(u001+u211-u011-u201)
+ dy3/(8*dx*yspzs)*(u000+u210-u010-u200)
+ dy*(2*dx*dx+dy*dy)/(8*dz*xspys)*(w000+w012-w002-w010)
Expand All @@ -319,17 +319,17 @@ facediv_int (int ci, int cj, int ck, int nf,
+ dz3/(8*dy*zspxs)*(v100+v121-v101-v120);

fine[2](fi , fj+1, fk+1, nf) = Real(0.5)*(w010+w012)
+ dz*(2*dy*dy+dz*dz)/(8*dx*yspzs)*(u000+u201-u001-u200)
+ dz3/(8*dx*yspzs)*(u010+u211-u011-u210)
+ dz*(2*dx*dx+dz*dz)/(8*dy*zspxs)*(v100+v121-v101-v120)
+ dz3/(8*dy*zspxs)*(v000+v021-v001-v020);

fine[2](fi+1, fj , fk+1, nf) = Real(0.5)*(w100+w102)
+ dz*(2*dy*dy+dz*dz)/(8*dx*yspzs)*(u010+u211-u011-u210)
+ dz3/(8*dx*yspzs)*(u000+u201-u001-u200)
+ dz*(2*dx*dx+dz*dz)/(8*dy*zspxs)*(v000+v021-v001-v020)
+ dz3/(8*dy*zspxs)*(v100+v121-v101-v120);

fine[2](fi+1, fj , fk+1, nf) = Real(0.5)*(w100+w102)
+ dz*(2*dy*dy+dz*dz)/(8*dx*yspzs)*(u000+u201-u001-u200)
+ dz3/(8*dx*yspzs)*(u010+u211-u011-u210)
+ dz*(2*dx*dx+dz*dz)/(8*dy*zspxs)*(v100+v121-v101-v120)
+ dz3/(8*dy*zspxs)*(v000+v021-v001-v020);

fine[2](fi+1, fj+1, fk+1, nf) = Real(0.5)*(w110+w112)
+ dz*(2*dy*dy+dz*dz)/(8*dx*yspzs)*(u010+u211-u011-u210)
+ dz3/(8*dx*yspzs)*(u000+u201-u001-u200)
Expand Down
3 changes: 3 additions & 0 deletions Src/AmrCore/AMReX_Interpolater.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -565,6 +565,9 @@ public:
* of the divergence of the underlying crse cell. All fine cells overlying
* a given coarse cell will have the same divergence, even when the coarse
* grid divergence is spatially varying.
* Based on Vanella et. al. (doi:10.1016/j.jcp.2010.05.003, section 3.2),
* but solves the interior closure problem using least squares with an
* initial guess equal to the average of fine face values across the cell.
*/
class FaceDivFree
:
Expand Down
42 changes: 29 additions & 13 deletions Tests/DivFreePatch/main.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -94,11 +94,14 @@ void CoarsenToFine(MultiFab& div_refined_coarse,


Real MFdiff(const MultiFab& lhs, const MultiFab& rhs,
int strt_comp, int num_comp, int nghost, const std::string name = "")
int strt_comp, int num_comp, int nghost, const std::string name = "",
bool relative = false)
{
MultiFab temp(lhs.boxArray(), lhs.DistributionMap(), lhs.nComp(), nghost);
Copy(temp, lhs, strt_comp, strt_comp, num_comp, nghost);
temp.minus(rhs, strt_comp, num_comp, nghost);
if (relative) {
temp.divide(rhs, strt_comp, num_comp, nghost); }

if (name != "")
{ amrex::VisMF::Write(temp, std::string("pltfiles/" + name)); }
Expand Down Expand Up @@ -301,7 +304,7 @@ void main_main ()
// Setup initial value on the coarse faces.
for (int i=0; i<AMREX_SPACEDIM; ++i)
{
setupMF(c_mf_faces[i], 1);
setupMF(c_mf_faces[i], i);
}

amrex::UtilCreateDirectoryDestructive("pltfiles");
Expand Down Expand Up @@ -403,16 +406,21 @@ void main_main ()
amrex::VisMF::Write(f_mf_faces_wg[1], std::string("pltfiles/fwgy"));,
amrex::VisMF::Write(f_mf_faces_wg[2], std::string("pltfiles/fwgz")); );

amrex::Print() << " Max InterpFromCoarse divergence error: "
<< MFdiff(div_fine, div_refined_coarse, 0, 1, nghost_f, "diff") << '\n';
amrex::Print() << " Max InterpFromCoarse divergence error: absolute relative\n "
<< " "
<<MFdiff(div_fine, div_refined_coarse, 0, 1, nghost_f, "diff")
<< " "
<<MFdiff(div_fine, div_refined_coarse, 0, 1, nghost_f, "", true)
<< '\n';

// ***************************************************************

// Change coarse values, save current fine values for checking
// the final result.
amrex::Print()<<" Cyclically permute the velocity components. Should not change div\n";
for (int i=0; i<AMREX_SPACEDIM; ++i)
{
setupMF(c_mf_faces[i], 2, BoxArray(amrex::coarsen(f_geom.Domain(), ratio).convert(c_mf_faces[i].ixType())));
setupMF(c_mf_faces[i], (i+1)%3, BoxArray(amrex::coarsen(f_geom.Domain(), ratio).convert(c_mf_faces[i].ixType())));
Copy(f_mf_copy[i], f_mf_faces[i], 0, 0, 1, 0);
}

Expand All @@ -424,8 +432,13 @@ void main_main ()
CoarsenToFine(div_refined_coarse, div_coarse, c_geom, f_geom_all, ratio);
amrex::VisMF::Write(div_coarse, std::string("pltfiles/coarsetofineB"));

amrex::Print() << " Checking new adjustment hasn't changed solution on fine region: "
<< MFdiff(div_fine, div_refined_coarse, 0, 1, nghost_f, "diffFP") << '\n';
amrex::Print() << " Change to divergence on fine region: absolute relative\n "
<< " "
<<MFdiff(div_fine, div_refined_coarse, 0, 1, nghost_f)
<< " "
<<MFdiff(div_fine, div_refined_coarse, 0, 1, nghost_f, "", true)
<< '\n';



// ***************************************************************
Expand Down Expand Up @@ -485,20 +498,23 @@ void main_main ()
// ================================================

// Checking fine valid cells are identical.
Real max_diff = 0;
for (int i=0; i<AMREX_SPACEDIM; ++i)
{
Real max_i = std::abs( MFdiff(f_mf_copy[i], f_mf_faces[i], 0, 1, 0) );
max_diff = (max_diff > max_i) ? max_diff : max_i;
amrex::Print() << " Fine valid region maximum change, comp "<<i<<": "
<<MFdiff(f_mf_copy[i], f_mf_faces[i], 0, 1, 0)
<< '\n';
}
amrex::Print() << " Fine values maximum change: " << max_diff << '\n';

// Check fine divergence = coarse divergence in ghost cells.
calcDiv(f_mf_faces, div_fine, f_geom.CellSizeArray());
amrex::VisMF::Write(div_fine, std::string("pltfiles/fineFP"));

amrex::Print() << " Max FillPatchTwoLevels divergence error: "
<< MFdiff(div_fine, div_refined_coarse, 0, 1, nghost_f, "diffFP") << '\n';
amrex::Print() << " Max FillPatchTwoLevels divergence error: absolute relative\n "
<< " "
<<MFdiff(div_fine, div_refined_coarse, 0, 1, nghost_f, "diffFP")
<< " "
<<MFdiff(div_fine, div_refined_coarse, 0, 1, nghost_f, "", true)
<< '\n';

for (int i=0; i<AMREX_SPACEDIM; ++i)
{
Expand Down

0 comments on commit a8c5255

Please sign in to comment.