Commit 4eec5991 authored by Philippe Delandmeter's avatar Philippe Delandmeter
Browse files

slim3d Interface. Not all boundary conditions are ready yet

parent 59eeb233
......@@ -841,6 +841,75 @@ def get_distance(region, physical_tags):
dist = distance.distance(region._mesh, physical_tags, evaluationPoints=xyz)
return np.array(dist)
def extrude(mesh_file_name, bath_file_name, nb_layers=None, z_layers=None, layers_function=None, mesh_file_name_out='', factor_show=0, periodicity=None):
if (nb_layers and z_layers) or (nb_layers and layers_function) or (z_layers and layers_function):
dgpy.Msg.Fatal("Only one among nb_layers, z_layers and layers_function can be defined!")
if periodicity:
shiftOperation = periodicity[0]
cutTags = periodicity[1]
pasteTags = periodicity[2]
mapFilename = periodicity[3]
slim_private.periodicMap.periodicMap(mesh_file_name, '_tmp_mesh_file', 1, shiftOperation, cutTags, pasteTags, mapFilename)
slim_private.shutil.copyfile('_tmp_mesh_file', mesh_file_name)
groups = dgpy.dgGroupCollection(mesh_file_name)
bathFunc = slim_private._load_function(bath_file_name, groups)
bathDof = dgpy.dgDofContainer(groups, 1)
bathDof.interpolate(bathFunc)
globNodeId2Bath = {}
dataProx = dgpy.fullMatrixDouble()
for iGroup in range(groups.getNbElementGroups()) :
group = groups.getElementGroup(iGroup)
for iElement in range(group.getNbElements()) :
el = group.getElement(iElement)
bathDof.getElementProxy(iGroup,iElement, dataProx)
for iPt in range(group.getNbNodes()) :
vId = groups.mesh().gmshVertexId(el.vertex(iPt))
globNodeId2Bath[vId] = dataProx(iPt,0)
def layersForZ(h) :
for i in range(len(z_layers)) :
if z_layers[i] > h :
return i
return len(z_layers)
def getLayers(e, v) :
h = globNodeId2Bath[v[3]]
if nb_layers:
return [z * h /nb_layers for z in range(nb_layers + 1)]
elif z_layers:
return [z_layers[i] for i in range(layersForZ(h))]
elif layers_function:
return layers_function(v[0], v[1], h)
else:
dgpy.Msg.Fatal("Neither nb_layers nor z_layers nor layers_function was defined")
def getLayersShow(e, v) :
layers = getLayers(e, v)
for i in range(len(layers)):
layers[i] = factor_show * layers[i]
return layers
if mesh_file_name_out:
mesh3d = mesh_file_name_out
mesh3dShow = mesh_file_name_out[:-4] + '_show.msh'
else:
mesh3d = mesh_file_name[:-4] + '_3d.msh'
mesh3dShow = mesh_file_name[:-4] + '_3d_show.msh'
if periodicity:
slim_private.extrudePeriodic.extrudePeriodic(slim_private.mesh.mesh(mesh_file_name), getLayers, mapFilename, mapFilename).write(mesh3d)
else:
slim_private.extrude.extrude(slim_private.mesh.mesh(mesh_file_name), getLayers).write(mesh3d)
if factor_show > 0:
slim_private.extrude.extrude(slim_private.mesh.mesh(mesh_file_name), getLayersShow).write(mesh3dShow)
return
def fetch_ftp(file_name, data_dir=None):
"""This function downloads some data from the SLIM servers if available. The data is downloaded in the user global_data_directory. Data is only downloaded if it does not exist locally or is older than the server data.
......
......@@ -7,6 +7,11 @@ import sys
from . import dgftp
from os import path, makedirs, remove
import shutil
import dgpy.mesh as mesh
import dgpy.extrude as extrude
import dgpy.periodicMap as periodicMap
import dgpy.extrudePeriodic as extrudePeriodic
from dgpy.slim3d_setup import slim3d_setup
import time
class _velocity_extractor :
......@@ -189,8 +194,8 @@ def _is_string(s) :
if sys.version_info >= (3, 0):
return isinstance(s, str)
else :
return isinstance(s, basestring)
return isinstance(s, basestring)
def _check_if_partionned(mesh_file_name):
f = open(mesh_file_name,'r')
line = f.readline()
......@@ -209,6 +214,7 @@ def _cpu_time(seconds):
h, m = divmod(m, 60)
return "%d:%02d:%02d" % (h, m, s)
<<<<<<< 158678f887f029c3ceab77a02a388f6a13e9187c
def _new_default_linear_system(domain, equation) :
try :
sys = dgpy.linearSystemPETScBlockDouble()
......@@ -224,3 +230,31 @@ def _new_default_linear_system(domain, equation) :
except :
dgpy.Msg.Fatal("No linear system available")
return sys
=======
def _findTopAndBottomTags(mesh_file_name):
bottomTags = []
topTags = []
fin = open(mesh_file_name, 'r')
for i in range(5):
l = fin.readline()
num = int(l)
for i in range(num):
l = fin.readline()
(trash1, tag, trash2) = (l.split("\""))
if tag[0:3] == 'top':
topTags.append(tag)
if tag[0:6] == 'bottom':
bottomTags.append(tag)
fin.close()
return (topTags, bottomTags)
exit(-1)
class OpenBnd_3d:
def __init__(self, tag, uv, eta, salinity, temperature):
self.tag = tag
self.uv = uv
self.eta = eta
self.salinity = salinity
self.temperature = temperature
>>>>>>> slim3d Interface. Not all boundary conditions are ready yet
......@@ -42,6 +42,7 @@ namespace std {
%apply double &OUTPUT { double *SDev, double *TDev}
%extend slim3dParameters{
static void setG(double f){slim3dParameters::g = f;}
static double getG(){return slim3dParameters::g;}
static void setRho0(double f){slim3dParameters::rho0 = f;}
static void setVonKarman(double f){slim3dParameters::vonKarman = f;}
static void setNuv0(double f){slim3dParameters::nuv0 = f;}
......
......@@ -1035,10 +1035,10 @@ void slim3dSolverEquations::initialize()
eta2dEq->setLaxFriedrichs(lf);
uv2dEq->setLaxFriedrichs(lf);
uvHorDivEq->setBathymetryMinCG(bathCGMin2d);
uvHorDivEq->setUV(uv3d);
uvHorDivEq->setEta(eta2d);
uvHorDivEq->setBoundarySymmetry(_solver->physicalTags2d);
//uvHorDivEq->setBathymetryMinCG(bathCGMin2d);
//uvHorDivEq->setUV(uv3d);
//uvHorDivEq->setEta(eta2d);
//uvHorDivEq->setBoundarySymmetry(_solver->physicalTags2d);
if (_solver->getSolveUVGrad()){
uvGradEq = new dgFEVectorGradient(uv3d);
......@@ -1047,13 +1047,13 @@ void slim3dSolverEquations::initialize()
uvGradEq->setBoundarySymmetry(_solver->physicalTags2d);
}
depthIntegratorRho->setFunction(funcs->rhoFunc);
for (size_t i = 0; i < _solver->physicalTags2d.size(); i++ ) {
if ( _solver->physicalTags2d[i] != "top" )
depthIntegratorRho->setBoundary0Flux(_solver->physicalTags2d[i]);
else
depthIntegratorRho->setBoundarySymmetry(_solver->physicalTags2d[i], "");
}
//depthIntegratorRho->setFunction(funcs->rhoFunc);
//for (size_t i = 0; i < _solver->physicalTags2d.size(); i++ ) {
// if ( _solver->physicalTags2d[i] != "top" )
// depthIntegratorRho->setBoundary0Flux(_solver->physicalTags2d[i]);
// else
// depthIntegratorRho->setBoundarySymmetry(_solver->physicalTags2d[i], "");
//}
if (_solver->getSolveGMVel()){
GMEq = new dgEddyTransportFlux(funcs->eadySpeedFunc, funcs->GMIndepTermFunc, dofs->nbDof3d->getFunction());
......
......@@ -400,9 +400,10 @@ void smagorinskyLimited::operator()(functorCache &cm, fullMatrix<double> &val) c
}
}
okubo::okubo(double factor, const functor *area) {
okubo::okubo(double factor, double max, const functor *area) {
_factor = factor;
_area = area;
_max = max;
}
void okubo::operator()(functorCache &cm, fullMatrix<double> &val) const {
......@@ -414,7 +415,7 @@ void okubo::operator()(functorCache &cm, fullMatrix<double> &val) const {
/* typical mesh size = sqrt of area */
delta = sqrt(area(i,0));
/* Okubo formulation */
val(i, 0) = _factor * pow(delta, 1.15);
val(i, 0) = std::min(_max, _factor * pow(delta, 1.15));
}
}
......
......@@ -542,7 +542,7 @@ class smagorinskyLimited : public functor {
/// \brief Class that computes the horizontal turbulent diffusivity, based on Okubo (1971)
class okubo : public functor {
/// \brief factor parameter
double _factor;
double _factor, _max;
/// \brief area of each triangle
const functor *_area;
public:
......@@ -551,7 +551,7 @@ class okubo : public functor {
* \param factor Factor parameter.
* \param area Area (of each triangle) function
*/
okubo(double factor, const functor *area);
okubo(double factor, double max, const functor *area);
/**
* \brief Computing the horizontal turbulent diffusivity.
*/
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment