[cig-commits] [commit] tomographic: initial version (afafa83)
cig_noreply at geodynamics.org
cig_noreply at geodynamics.org
Wed Jun 11 11:58:59 PDT 2014
Repository : https://github.com/geodynamics/aspect
On branch : tomographic
Link : https://github.com/geodynamics/aspect/compare/0000000000000000000000000000000000000000...475fdb71218169cbf9959396a1c1c7ad3deaa43d
>---------------------------------------------------------------
commit afafa83dfdfbfc1551aefbe55d82b28f0bfda15f
Author: Timo Heister <timo.heister at gmail.com>
Date: Sun Jun 8 10:44:39 2014 -0400
initial version
>---------------------------------------------------------------
afafa83dfdfbfc1551aefbe55d82b28f0bfda15f
{benchmark/solkz => tomo}/CMakeLists.txt | 4 +-
tomo/convert.cc | 29 +++
tomo/notes.txt | 1 +
cookbooks/convection-box-3d.prm => tomo/tomo3d.prm | 25 +-
tomo/tomo_amr.cc | 280 +++++++++++++++++++++
5 files changed, 331 insertions(+), 8 deletions(-)
diff --git a/benchmark/solkz/CMakeLists.txt b/tomo/CMakeLists.txt
similarity index 77%
copy from benchmark/solkz/CMakeLists.txt
copy to tomo/CMakeLists.txt
index 8b69d1c..68c714f 100644
--- a/benchmark/solkz/CMakeLists.txt
+++ b/tomo/CMakeLists.txt
@@ -3,8 +3,8 @@ CMAKE_MINIMUM_REQUIRED(VERSION 2.8.8)
FIND_PACKAGE(Aspect REQUIRED HINTS ${ASPECT_DIR} ../ ../../ $ENV{ASPECT_DIR})
DEAL_II_INITIALIZE_CACHED_VARIABLES()
-SET(TARGET "solkz")
+SET(TARGET "tomo")
PROJECT(${TARGET})
-ADD_LIBRARY(${TARGET} SHARED solkz.cc)
+ADD_LIBRARY(${TARGET} SHARED tomo_amr.cc)
ASPECT_SETUP_PLUGIN(${TARGET})
diff --git a/tomo/convert.cc b/tomo/convert.cc
new file mode 100644
index 0000000..91f6019
--- /dev/null
+++ b/tomo/convert.cc
@@ -0,0 +1,29 @@
+#include <iostream>
+#include <fstream>
+
+int main()
+{
+ char c;
+ std::ifstream f_in("mt_gambier_512.ubc");
+ std::ofstream f_out("mt_gambier_512.data");
+
+ for (unsigned int x=0;x<512;++x)
+ for (unsigned int y=0;y<512;++y)
+ for (unsigned int z=0;z<512;++z)
+ {
+ char c;
+ f_in >> c;
+ // std::cout << (int)c << " ";
+
+ if (c==1)
+ c=0xff;
+
+ if (z<300 && x<300)
+ c=0;
+
+
+ f_out << c;
+ }
+
+
+}
diff --git a/tomo/notes.txt b/tomo/notes.txt
new file mode 100644
index 0000000..c21053b
--- /dev/null
+++ b/tomo/notes.txt
@@ -0,0 +1 @@
+http://people.physics.anu.edu.au/~aps110/network_comparison/
\ No newline at end of file
diff --git a/cookbooks/convection-box-3d.prm b/tomo/tomo3d.prm
similarity index 92%
copy from cookbooks/convection-box-3d.prm
copy to tomo/tomo3d.prm
index e53d4b8..1a04fa2 100644
--- a/cookbooks/convection-box-3d.prm
+++ b/tomo/tomo3d.prm
@@ -1,3 +1,6 @@
+
+set Additional shared libraries = ./libtomo.so
+
# A description of convection in a 3d box. See the manual for more information.
# At the top, we define the number of space dimensions we would like to
@@ -8,7 +11,7 @@ set Dimension = 3
# time system we want to work in and what the end time is. We
# also designate an output directory.
set Use years in output instead of seconds = false
-set End time = 1.0
+set End time = 0.0
set Output directory = output
# Then there are variables that describe the tolerance of
@@ -16,7 +19,7 @@ set Output directory = output
# be normalized. Here, we choose a zero average pressure
# at the surface of the domain (for the current geometry, the
# surface is defined as the top boundary).
-set Linear solver tolerance = 1e-7
+set Linear solver tolerance = 1e-10
set Temperature solver tolerance = 1e-10
set Pressure normalization = surface
@@ -146,12 +149,13 @@ end
# on. We here use a globally refined mesh without
# adaptive mesh refinement.
subsection Mesh refinement
- set Initial global refinement = 3
- set Initial adaptive refinement = 3
- set Time steps between mesh refinement = 15
+ set Initial global refinement = 10
+ set Initial adaptive refinement = 0
+ set Time steps between mesh refinement = 0
set Additional refinement times = 0.003
+set Strategy = TomoAMR
end
@@ -170,10 +174,19 @@ subsection Postprocess
set List of postprocessors = velocity statistics, temperature statistics, heat flux statistics, visualization
subsection Visualization
- set Time between graphical output = 0.0001
+ set Time between graphical output = 0
end
end
+subsection Compositional fields
+ set Number of fields = 1
+end
+
+subsection Compositional initial conditions
+ set Model name = tomo
+end
+
+
subsection Checkpointing
# The number of timesteps between performing checkpoints. If 0 and time
diff --git a/tomo/tomo_amr.cc b/tomo/tomo_amr.cc
new file mode 100644
index 0000000..ad25d81
--- /dev/null
+++ b/tomo/tomo_amr.cc
@@ -0,0 +1,280 @@
+/*
+ Copyright (C) 2011 - 2014 by the authors of the ASPECT code.
+
+ This file is part of ASPECT.
+
+ ASPECT is free software; you can redistribute it and/or modify
+ it under the terms of the GNU General Public License as published by
+ the Free Software Foundation; either version 2, or (at your option)
+ any later version.
+
+ ASPECT is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ GNU General Public License for more details.
+
+ You should have received a copy of the GNU General Public License
+ along with ASPECT; see the file doc/COPYING. If not see
+ <http://www.gnu.org/licenses/>.
+*/
+
+
+
+#include <aspect/mesh_refinement/interface.h>
+#include <aspect/compositional_initial_conditions/interface.h>
+#include <aspect/simulator_access.h>
+
+
+
+#include <deal.II/base/quadrature_lib.h>
+#include <deal.II/fe/fe_values.h>
+#include <math.h>
+
+
+using namespace dealii;
+
+class Data
+ {
+ public:
+
+
+ Data()
+ {
+ sizes[0]=512;
+ sizes[1]=512;
+ sizes[2]=512;
+ data.resize(512*512*512);
+
+ std::ifstream in("mt_gambier_512.data",std::ios::in | std::ios::binary);
+ in.read((char*)&data[0],512*512*512);
+ }
+
+
+ void compute_min_max(Point<2> &corner1, Point<2> &corner2, double &minval, double &maxval, double &mean, double &err) const
+ {
+
+ }
+
+
+ void compute_min_max(Point<3> &corner1, Point<3> &corner2, double &minval, double &maxval, double &mean, double &err) const
+ {
+ unsigned int x1 = std::min((unsigned int)(sizes[0]*corner1[0]), sizes[0]-1);
+ unsigned int x2 = std::min((unsigned int)(sizes[0]*corner2[0]), sizes[0]-1);
+ unsigned int y1 = std::min((unsigned int)(sizes[1]*corner1[1]), sizes[1]-1);
+ unsigned int y2 = std::min((unsigned int)(sizes[1]*corner2[1]), sizes[1]-1);
+ unsigned int z1 = std::min((unsigned int)(sizes[2]*corner1[2]), sizes[2]-1);
+ unsigned int z2 = std::min((unsigned int)(sizes[2]*corner2[2]), sizes[2]-1);
+
+ mean = 0.0;
+ unsigned int count = 0;
+ for (unsigned int x=x1;x<=x2;++x)
+ for (unsigned int y=y1;y<=y2;++y)
+ for (unsigned int z=z1;z<=z2;++z)
+ {
+ mean += (double)(data[x+sizes[0]*y + sizes[0]*sizes[1]*z])/255.0;
+ count += 1;
+ }
+
+ mean /= count;
+
+
+ err = 0;
+
+ minval = 1.0;
+ maxval = 0.0;
+ for (unsigned int x=x1;x<=x2;++x)
+ for (unsigned int y=y1;y<=y2;++y)
+ for (unsigned int z=z1;z<=z2;++z)
+ {
+ double val = (double)(data[x+sizes[0]*y + sizes[0]*sizes[1]*z])/255.0;
+ minval = std::min(minval, val);
+ maxval = std::max(maxval, val);
+ err += (val-mean)*(val-mean);
+ }
+
+ //err /= count;
+ //std::cout << minval << " " << maxval << " " << mean << " " << err << std::endl;
+
+ mean = err;
+
+ }
+
+ unsigned int sizes[3];
+ std::vector<unsigned char> data;
+
+
+ };
+
+Data data;
+
+
+
+namespace aspect
+{
+ namespace MeshRefinement
+ {
+
+
+
+
+
+ template <int dim>
+ class TomoAMR : public Interface<dim>,
+ public SimulatorAccess<dim>
+ {
+ public:
+ /**
+ * After cells have been marked for coarsening/refinement, apply
+ * additional criteria independent of the error estimate.
+ *
+ */
+ virtual
+ void
+ tag_additional_cells () const;
+
+ /**
+ * Declare the parameters this class takes through input files.
+ */
+ static
+ void
+ declare_parameters (ParameterHandler &prm);
+
+ /**
+ * Read the parameters this class declares from the parameter
+ * file.
+ */
+ virtual
+ void
+ parse_parameters (ParameterHandler &prm);
+
+ private:
+ };
+
+
+ template <int dim>
+ void
+ TomoAMR<dim>::tag_additional_cells () const
+ {
+ for (typename Triangulation<dim>::active_cell_iterator
+ cell = this->get_triangulation().begin_active();
+ cell != this->get_triangulation().end(); ++cell)
+ {
+ if (cell->is_locally_owned())
+ {
+ Point<dim> p1 = cell->center();
+ Point<dim> p2 = cell->center();
+ p1[0] -= cell->extent_in_direction(0)*0.5;
+ p2[0] += cell->extent_in_direction(0)*0.5;
+ p1[1] -= cell->extent_in_direction(1)*0.5;
+ p2[1] += cell->extent_in_direction(1)*0.5;
+ p1[2] -= cell->extent_in_direction(2)*0.5;
+ p2[2] += cell->extent_in_direction(2)*0.5;
+
+ double minval, maxval, mean, err;
+ data.compute_min_max(p1, p2, minval, maxval, mean, err);
+ if (maxval-minval < 0.1 || (err) < 5000.0)
+ cell->clear_refine_flag ();
+ }
+ }
+ }
+
+ template <int dim>
+ void
+ TomoAMR<dim>::
+ declare_parameters (ParameterHandler &prm)
+ {
+ prm.enter_subsection("Mesh refinement");
+ {
+
+ prm.enter_subsection("TomoAMR");
+ {
+
+ }
+ prm.leave_subsection();
+ }
+ prm.leave_subsection();
+ }
+
+ template <int dim>
+ void
+ TomoAMR<dim>::parse_parameters (ParameterHandler &prm)
+ {
+ prm.enter_subsection("Mesh refinement");
+ {
+ prm.enter_subsection("TomoAMR");
+
+
+ prm.leave_subsection();
+ }
+ prm.leave_subsection();
+ }
+
+ }
+
+
+
+
+ template <int dim>
+ class TomoInitial : public CompositionalInitialConditions::Interface<dim>,
+ public SimulatorAccess<dim>
+ {
+ public:
+ /**
+ * Constructor.
+ */
+
+ /**
+ * Return the initial composition as a function of position and number
+ * of compositional field.
+ */
+ virtual
+ double initial_composition (const Point<dim> &position, const unsigned int n_comp) const
+ {
+ double minval, maxval, mean, err;
+ Point<dim> p = position;
+ data.compute_min_max(p, p, minval, maxval, mean, err);
+ return (minval+maxval)*0.5;
+
+ }
+
+ /**
+ * Declare the parameters this class takes through input files. The
+ * default implementation of this function does not describe any
+ * parameters. Consequently, derived classes do not have to overload
+ * this function if they do not take any runtime parameters.
+ */
+ static
+ void
+ declare_parameters (ParameterHandler &prm)
+ {}
+
+ /**
+ * Read the parameters this class declares from the parameter file.
+ * The default implementation of this function does not read any
+ * parameters. Consequently, derived classes do not have to overload
+ * this function if they do not take any runtime parameters.
+ */
+ virtual
+ void
+ parse_parameters (ParameterHandler &prm)
+ {}
+
+ private:
+
+ };
+
+}
+
+// explicit instantiations
+namespace aspect
+{
+ namespace MeshRefinement
+ {
+ ASPECT_REGISTER_MESH_REFINEMENT_CRITERION(TomoAMR,
+ "TomoAMR",
+ "TODO")
+ }
+ ASPECT_REGISTER_COMPOSITIONAL_INITIAL_CONDITIONS(TomoInitial,
+ "tomo",
+ "")
+}
More information about the CIG-COMMITS
mailing list