https://github.com/geodynamics/citcoms
Tip revision: f8f8ac9d2debc000a1b14bbe2468aa3a3b16d3bd authored by Rajesh Kommu on 18 April 2014, 22:33:36 UTC
Minor cleanup related to caps_per_proc removal
Minor cleanup related to caps_per_proc removal
Tip revision: f8f8ac9
AreaWeightedNormal.cc
// -*- C++ -*-
//
//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
//
//<LicenseText>
//
// CitcomS.py by Eh Tan, Eun-seo Choi, and Pururav Thoutireddy.
// Copyright (C) 2002-2005, California Institute of Technology.
//
// This program 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 of the License, or
// (at your option) any later version.
//
// This program 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 this program; if not, write to the Free Software
// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
//
//</LicenseText>
//
//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
//
#include "config.h"
#include <cmath>
#include <iostream>
#include "global_defs.h"
#include "journal/diagnostics.h"
#include "Boundary.h"
#include "Exchanger/UtilTemplate.h"
#include "Exchanger/Array2D.h"
#include "Exchanger/Sink.h"
#include "AreaWeightedNormal.h"
extern "C" {
// for definition of SIDE_NORTH etc.
#include "element_definitions.h"
}
//using Exchanger::Array2D;
using Exchanger::DIM;
using Exchanger::Sink;
AreaWeightedNormal::AreaWeightedNormal(const MPI_Comm& comm,
const Boundary& boundary,
const Sink& sink,
const All_variables* E) :
nwght(boundary.size(), 0)
{
journal::debug_t debug("CitcomS-Exchanger");
debug << journal::at(__HERE__) << journal::endl;
computeWeightedNormal(boundary, E);
double inv_total_area = 1 / computeTotalArea(comm, sink);
// normalize
for(int n=0; n<sink.size(); n++) {
int i = sink.meshNode(n);
for(int j=0; j<DIM; j++)
nwght[j][i] *= inv_total_area;
}
}
AreaWeightedNormal::~AreaWeightedNormal()
{}
void AreaWeightedNormal::imposeConstraint(Velo& V,
const MPI_Comm& comm,
const Sink& sink,
const All_variables* E) const
{
journal::debug_t debug("CitcomS-Exchanger");
debug << journal::at(__HERE__) << journal::endl;
double outflow = computeOutflow(V, comm, sink);
reduceOutflow(V, outflow, sink);
double new_outflow = computeOutflow(V, comm, sink);
if(E->parallel.me == 0) {
fprintf(stderr, "Net outflow amended from %e to %e\n",
outflow, new_outflow);
fprintf(E->fp, "Net outflow amended from %e to %e\n",
outflow, new_outflow);
}
}
// private functions
void AreaWeightedNormal::computeWeightedNormal(const Boundary& boundary,
const All_variables* E)
{
const int nodes_per_element = 8;
const int vpoints_per_side = 4;
// converting normal vector [-1, 0, 1] to side index,
const int side_normal[DIM][DIM] = {{SIDE_NORTH, -1, SIDE_SOUTH},
{SIDE_WEST, -1, SIDE_EAST},
{SIDE_BOTTOM, -1, SIDE_TOP}};
/* For each node belong to boundary elements, check all 6 faces,
* if the node is on the face and the face is part of the boundary,
* compute the surface area by summing the 2D determinants. */
for (int m=1; m<=E->sphere.caps_per_proc; m++)
for (int es=1; es<=E->boundary.nel; es++) {
int el = E->boundary.element[m][es];
for(int n=1; n<=nodes_per_element; n++) {
int node = E->ien[m][el].node[n];
int bnode = boundary.bnode(node);
// skipping nodes not on the boundary
if(bnode < 0) continue;
for(int j=0; j<DIM; j++) {
// normal is either -1, 1, or 0
int normal = boundary.normal(j, bnode);
if(normal) {
int side = side_normal[j][normal+1];
for(int k=1; k<=vpoints_per_side; k++)
nwght[j][bnode] += normal * E->boundary.det[m][side][k][es] / vpoints_per_side;
}
} // end of loop over dim
} // end of loop over element nodes
} // end of loop over elements
}
double AreaWeightedNormal::computeTotalArea(const MPI_Comm& comm,
const Sink& sink) const
{
double total_area = 0;
for(int n=0; n<sink.size(); n++) {
int i = sink.meshNode(n);
for(int j=0; j<DIM; j++)
total_area += std::abs(nwght[j][i]);
}
Exchanger::util::gatherSum(comm, total_area);
journal::info_t info("CitcomS-AreaWeightedNormal-area");
info << journal::at(__HERE__)
<< "Total surface area: "
<< total_area << journal::endl;
/** debug **
for(int j=0; j<DIM; j++) {
double in, out;
in = out = 0;
for(int n=0; n<sink.size(); n++) {
int i = sink.meshNode(n);
if(nwght[j][i] > 0) out += nwght[j][i];
else in += nwght[j][i];
}
Exchanger::util::gatherSum(comm, in);
Exchanger::util::gatherSum(comm, out);
journal::info_t info("CitcomS-AreaWeightedNormal-area");
info << journal::at(__HERE__)
<< "Partial surface area: " << j << " "
<< in << " " << out << journal::endl;
}
*/
return total_area;
}
double AreaWeightedNormal::computeOutflow(const Velo& V,
const MPI_Comm& comm,
const Sink& sink) const
{
/* integrate dot(V,n) over the surface */
double outflow = 0;
for(int n=0; n<sink.size(); n++) {
int i = sink.meshNode(n);
for(int j=0; j<DIM; j++)
outflow += V[j][n] * nwght[j][i];
}
Exchanger::util::gatherSum(comm, outflow);
return outflow;
}
int AreaWeightedNormal::sign(double number) const
{
return (number > 0) ? 1 : ((number < 0) ? -1 : 0);
}
void AreaWeightedNormal::reduceOutflow(Velo& V, double outflow,
const Sink& sink) const
{
for(int n=0; n<sink.size(); n++) {
int i = sink.meshNode(n);
for(int j=0; j<DIM; j++)
V[j][n] -= sign(nwght[j][i]) * outflow;
}
}
// version
// $Id$
// End of file