Revision 94cd779e7d1352a014fc3ca4c98eeb0a439a2e09 authored by Eh Tan on 14 September 2007, 23:39:59 UTC, committed by Eh Tan on 14 September 2007, 23:39:59 UTC
1 parent fedce4a
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" {
#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) :
size_(boundary.size()),
toleranceOutflow_(E->control.tole_comp),
nwght(size_, 0)
{
computeWeightedNormal(boundary, E);
double total_area = computeTotalArea(comm, sink);
normalize(total_area);
}
AreaWeightedNormal::~AreaWeightedNormal()
{}
void AreaWeightedNormal::imposeConstraint(Velo& V,
const MPI_Comm& comm,
const Sink& sink) const
{
journal::info_t info("CitcomS-AreaWeightedNormal-outflow");
double outflow = computeOutflow(V, comm, sink);
info << journal::at(__HERE__)
<< "Net outflow: "
<< outflow << journal::endl;
if (std::abs(outflow) > toleranceOutflow_) {
reduceOutflow(V, outflow, sink);
outflow = computeOutflow(V, comm, sink);
info << journal::at(__HERE__)
<< "Net outflow after correction (SHOULD BE ZERO !): "
<< outflow << journal::endl;
}
}
// private functions
void AreaWeightedNormal::computeWeightedNormal(const Boundary& boundary,
const All_variables* E)
{
const int nodes_per_element = 8;
const int vpoints_per_side = 4;
const int side_normal[DIM][DIM] = {{SIDE_NORTH, -1, SIDE_SOUTH},
{SIDE_WEST, -1, SIDE_EAST},
{SIDE_BOTTOM, -1, SIDE_TOP}};
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);
if(bnode < 0) continue;
for(int j=0; j<DIM; j++) {
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];
}
} // 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);
return total_area;
}
void AreaWeightedNormal::normalize(double total_area)
{
for(int i=0; i<size_; i++)
for(int j=0; j<DIM; j++)
nwght[j][i] /= total_area;
}
double AreaWeightedNormal::computeOutflow(const Velo& V,
const MPI_Comm& comm,
const Sink& sink) const
{
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
Computing file changes ...