[cig-commits] commit: Add in gradient detection.
Mercurial
hg at geodynamics.org
Mon Apr 11 12:14:39 PDT 2011
changeset: 148:cfcbe2a4f7e8
user: Walter Landry <wlandry at caltech.edu>
date: Mon Apr 11 01:12:26 2011 -0700
files: FACStokes.h FACStokes/FACStokes.C FACStokes/applyGradientDetector.C wscript
description:
Add in gradient detection.
diff -r a43af5d15c5a -r cfcbe2a4f7e8 FACStokes.h
--- a/FACStokes.h Mon Apr 11 01:10:32 2011 -0700
+++ b/FACStokes.h Mon Apr 11 01:12:26 2011 -0700
@@ -82,6 +82,18 @@ namespace SAMRAI {
int finest_level);
//@}
+
+ virtual void
+ applyGradientDetector(const tbox::Pointer<hier::BasePatchHierarchy> hierarchy,
+ const int level_number,
+ const double error_data_time,
+ const int tag_index,
+ const bool initial_time,
+ const bool uses_richardson_extrapolation);
+
+ void computeAdaptionEstimate(pdat::CellData<double>& estimate_data,
+ const pdat::CellData<double>& soln_cell_data)
+ const;
//@{ @name appu::VisDerivedDataStrategy virtuals
@@ -169,6 +181,8 @@ namespace SAMRAI {
*
* These are initialized in the constructor and never change.
*/
+
+ double d_adaptation_threshold;
public:
int p_id, cell_viscosity_id, edge_viscosity_id, dp_id, p_exact_id,
p_rhs_id, v_id, v_rhs_id;
diff -r a43af5d15c5a -r cfcbe2a4f7e8 FACStokes/FACStokes.C
--- a/FACStokes/FACStokes.C Mon Apr 11 01:10:32 2011 -0700
+++ b/FACStokes/FACStokes.C Mon Apr 11 01:12:26 2011 -0700
@@ -139,6 +139,10 @@ namespace SAMRAI {
1 for coarsening
operator */);
+ d_adaptation_threshold=database->getDoubleWithDefault("adaption_threshold",
+ // 0.5);
+ 1e1);
+
/*
* Specify an implementation of solv::RobinBcCoefStrategy for the
* solver to use. We use the implementation
diff -r a43af5d15c5a -r cfcbe2a4f7e8 FACStokes/applyGradientDetector.C
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/FACStokes/applyGradientDetector.C Mon Apr 11 01:12:26 2011 -0700
@@ -0,0 +1,169 @@
+#include "FACStokes.h"
+#include "SAMRAI/geom/CartesianGridGeometry.h"
+#include "SAMRAI/pdat/MDA_Access.h"
+#include "SAMRAI/pdat/ArrayDataAccess.h"
+
+void SAMRAI::FACStokes::applyGradientDetector
+(const tbox::Pointer<hier::BasePatchHierarchy> hierarchy_,
+ const int ln,
+ const double ,
+ const int tag_index,
+ const bool ,
+ const bool )
+{
+ const tbox::Pointer<hier::PatchHierarchy> hierarchy__ = hierarchy_;
+ hier::PatchHierarchy& hierarchy = *hierarchy__;
+ tbox::Pointer<geom::CartesianGridGeometry>
+ grid_geometry_ = hierarchy.getGridGeometry();
+ hier::PatchLevel& level =
+ (hier::PatchLevel &) * hierarchy.getPatchLevel(ln);
+ hier::PatchLevel::Iterator pi;
+ int ntag = 0, ntotal = 0;
+ double maxestimate = 0;
+ for (pi.initialize(level); pi; pi++) {
+ hier::Patch& patch = **pi;
+ tbox::Pointer<hier::PatchData>
+ tag_data = patch.getPatchData(tag_index);
+ ntotal += patch.getBox().numberCells().getProduct();
+ if (tag_data.isNull()) {
+ TBOX_ERROR(
+ "Data index " << tag_index << " does not exist for patch.\n");
+ }
+ tbox::Pointer<pdat::CellData<int> > tag_cell_data_ = tag_data;
+ if (tag_cell_data_.isNull()) {
+ TBOX_ERROR("Data index " << tag_index << " is not cell int data.\n");
+ }
+ tbox::Pointer<hier::PatchData>
+ soln_data = patch.getPatchData(p_id);
+ if (soln_data.isNull()) {
+ TBOX_ERROR("Data index " << p_id
+ << " does not exist for patch.\n");
+ }
+ tbox::Pointer<pdat::CellData<double> > soln_cell_data_ = soln_data;
+ if (soln_cell_data_.isNull()) {
+ TBOX_ERROR("Data index " << p_id
+ << " is not cell int data.\n");
+ }
+ pdat::CellData<double>& soln_cell_data = *soln_cell_data_;
+ pdat::CellData<int>& tag_cell_data = *tag_cell_data_;
+ pdat::CellData<double> estimate_data(patch.getBox(),
+ 1,
+ hier::IntVector(d_dim, 0));
+ computeAdaptionEstimate(estimate_data,
+ soln_cell_data);
+ tag_cell_data.fill(0);
+ hier::Box::Iterator i;
+ for (i.initialize(patch.getBox()); i; i++) {
+ const pdat::CellIndex cell_index(*i);
+ if (maxestimate < estimate_data(cell_index)) maxestimate =
+ estimate_data(cell_index);
+
+ tbox::plog << "estimate "
+ << cell_index[0] << " "
+ << cell_index[1] << " "
+ << d_adaptation_threshold << " "
+ << estimate_data(cell_index) << " "
+ << std::boolalpha
+ << (estimate_data(cell_index) > d_adaptation_threshold) << " "
+ << "\n";
+ if (estimate_data(cell_index) > d_adaptation_threshold) {
+ tag_cell_data(cell_index) = 1;
+ ++ntag;
+ }
+ }
+ }
+ tbox::plog << "Adaption threshold is " << d_adaptation_threshold << "\n";
+ tbox::plog << "Number of cells tagged on level " << ln << " is "
+ << ntag << "/" << ntotal << "\n";
+ tbox::plog << "Max estimate is " << maxestimate << "\n";
+}
+
+
+void SAMRAI::FACStokes::computeAdaptionEstimate
+(pdat::CellData<double>& estimate_data,
+ const pdat::CellData<double>& soln_cell_data) const
+{
+ const int* lower = &estimate_data.getBox().lower()[0];
+ const int* upper = &estimate_data.getBox().upper()[0];
+ if (d_dim == tbox::Dimension(2)) {
+ MDA_AccessConst<double, 2, MDA_OrderColMajor<2> > co =
+ pdat::ArrayDataAccess::access<2, double>(soln_cell_data.getArrayData());
+ MDA_Access<double, 2, MDA_OrderColMajor<2> > es =
+ pdat::ArrayDataAccess::access<2, double>(estimate_data.getArrayData());
+ int i, j;
+ double estimate, est0, est1, est2, est3, est4, est5;
+ for (j = lower[1]; j <= upper[1]; ++j) {
+ for (i = lower[0]; i <= upper[0]; ++i) {
+ est0 =
+ tbox::MathUtilities<double>::Abs(co(i + 1, j)
+ + co(i - 1,j) - 2 * co(i, j));
+ est1 =
+ tbox::MathUtilities<double>::Abs(co(i, j + 1)
+ + co(i,j - 1) - 2 * co(i, j));
+ est2 = 0.5
+ * tbox::MathUtilities<double>::Abs(co(i + 1, j+ 1)
+ + co(i - 1, j - 1) - 2 * co(i, j));
+ est3 = 0.5
+ * tbox::MathUtilities<double>::Abs(co(i + 1, j- 1)
+ + co(i - 1, j + 1) - 2 * co(i, j));
+ est4 = tbox::MathUtilities<double>::Max(est0, est1);
+ est5 = tbox::MathUtilities<double>::Max(est2, est3);
+ estimate = tbox::MathUtilities<double>::Max(est4, est5);
+ es(i, j) = estimate;
+ }
+ }
+ }
+ if (d_dim == tbox::Dimension(3)) {
+ MDA_AccessConst<double, 3, MDA_OrderColMajor<3> > co =
+ pdat::ArrayDataAccess::access<3, double>(soln_cell_data.getArrayData());
+ MDA_Access<double, 3, MDA_OrderColMajor<3> > es =
+ pdat::ArrayDataAccess::access<3, double>(estimate_data.getArrayData());
+ // math::PatchCellDataOpsReal<double> cops;
+ // cops.printData( soln_cell_data_, soln_cell_data_->getGhostBox(), tbox::plog );
+ int i, j, k;
+ double estimate, est0, est1, est2, est3, est4, est5, est6, est7, est8,
+ esta, estb, estc, estd, este, estf, estg;
+ for (k = lower[2]; k <= upper[2]; ++k) {
+ for (j = lower[1]; j <= upper[1]; ++j) {
+ for (i = lower[0]; i <= upper[0]; ++i) {
+ est0 =
+ tbox::MathUtilities<double>::Abs(co(i + 1, j, k)
+ + co(i - 1,j,k) - 2 * co(i, j, k));
+ est1 =
+ tbox::MathUtilities<double>::Abs(co(i, j + 1, k)
+ + co(i,j - 1,k) - 2 * co(i, j, k));
+ est2 =
+ tbox::MathUtilities<double>::Abs(co(i, j, k + 1)
+ + co(i,j,k - 1) - 2 * co(i, j, k));
+ est3 = 0.5 * tbox::MathUtilities<double>::Abs(co(i,j + 1,k + 1)
+ + co(i, j - 1, k - 1)
+ - 2 * co(i, j, k));
+ est4 = 0.5 * tbox::MathUtilities<double>::Abs(co(i,j + 1,k - 1)
+ + co(i, j - 1, k + 1)
+ - 2 * co(i, j, k));
+ est5 = 0.5 * tbox::MathUtilities<double>::Abs(co(i + 1,j,k + 1)
+ + co(i - 1, j, k - 1)
+ - 2 * co(i, j, k));
+ est6 = 0.5 * tbox::MathUtilities<double>::Abs(co(i + 1,j,k - 1)
+ + co(i - 1, j, k + 1)
+ - 2 * co(i, j, k));
+ est7 = 0.5 * tbox::MathUtilities<double>::Abs(co(i + 1,j + 1,k)
+ + co(i - 1, j - 1, k)
+ - 2 * co(i, j, k));
+ est8 = 0.5 * tbox::MathUtilities<double>::Abs(co(i + 1,j - 1,k)
+ + co(i - 1, j + 1, k)
+ - 2 * co(i, j, k));
+ esta = tbox::MathUtilities<double>::Max(est0, est1);
+ estb = tbox::MathUtilities<double>::Max(est2, est3);
+ estc = tbox::MathUtilities<double>::Max(est4, est5);
+ estd = tbox::MathUtilities<double>::Max(est6, est7);
+ este = tbox::MathUtilities<double>::Max(esta, estb);
+ estf = tbox::MathUtilities<double>::Max(estc, estd);
+ estg = tbox::MathUtilities<double>::Max(este, estf);
+ estimate = tbox::MathUtilities<double>::Max(estg, est8);
+ es(i, j, k) = estimate;
+ }
+ }
+ }
+ }
+}
diff -r a43af5d15c5a -r cfcbe2a4f7e8 wscript
--- a/wscript Mon Apr 11 01:10:32 2011 -0700
+++ b/wscript Mon Apr 11 01:12:26 2011 -0700
@@ -13,6 +13,7 @@ def build(bld):
source = ['main.C',
'FACStokes/FACStokes.C',
'FACStokes/fix_viscosity.C',
+ 'FACStokes/applyGradientDetector.C',
'FACStokes/initializeLevelData.C',
'FACStokes/packDerivedDataIntoDoubleBuffer.C',
'FACStokes/resetHierarchyConfiguration.C',
More information about the CIG-COMMITS
mailing list