Commit 92608567 authored by Philippe Delandmeter's avatar Philippe Delandmeter Committed by Valentin Vallaeys
Browse files

slim3d: add testcases

new reStart + accurate parallel testcases
add columbia testcase
updated moving mesh & meshStar CG
slim3dNetCDFIO: update for z coordinates
parent f22f75c1
import slimPre
import numpy as np
pre_data_dir_base = 'data_'+slimPre.partition_nb()+'/'
pre_data_dir = pre_data_dir_base+slimPre.partition_id()+'/'
slimPre.make_directory(pre_data_dir)
mesh_file_name = slimPre.fetch_ftp('slim_data/columbia/meshes/columbia_v1.msh')
mesh_file_name = slimPre.partition_mesh(mesh_file_name, pre_data_dir_base+"mesh2d.msh")
tides_initial_time = '2006-05-1 07:45:00'
tides_final_time = '2006-05-14 00:00:00'
tides_time_step = 900
tides_physical_tags = ["sea_west", "sea_north", "sea_south"]
time_zone = -8.*3600
nb_layers = 10
print('Preprocessing Bathymetry')
bath_file_name = slimPre.fetch_ftp('slim_data/columbia/bath/bath_smooth_v1/bath_smooth_v1.idx', pre_data_dir_base)
bath_file = slimPre.fetch_ftp('slim_data/columbia/bath/bath_smooth_v1/bath_smooth_v1_COMP_0.msh', pre_data_dir_base)
zLayers = [0,1,3,5,10,20,50,100,300,800]
def layers_func(x, y, h):
n = nb_layers-1
f = nb_layers-1
det = (1-pow(f, (n+1.)/n))/(1-pow(f, 1./n))
layers = []
for z in range(n + 1):
a = 0
for i in range(z):
a+= pow(f, float(i)/n) * h /det
if (z < len(zLayers)):
a = min(a, zLayers[z])
layers.append(a)
layers.append(h)
return layers
print('Extruding mesh')
slimPre.extrude(mesh_file_name, bath_file_name, layers_function=layers_func, mesh_file_name_out=pre_data_dir_base+'mesh3d.msh', factor_show=500)
print('Loading 3D mesh')
mesh_file_name = pre_data_dir_base + 'mesh3d.msh'
mesh = slimPre.Mesh(mesh_file_name, mesh_proj="+proj=lcc +datum=NAD27 +lon_0=-120d30 +lat_1=46 +lat_2=44d20 +lat_0=43d40 +x_0=609601.2192024384 +y_0=0")
region_global = slimPre.Region(mesh)
lonlat_global = slimPre.Coordinate_system(region_global, data_proj="+proj=latlong +ellps=WGS84 +nadgrids=WO")
lonlat_coordinates = lonlat_global.coordinates
lonlat_coordinates_deg = lonlat_global.coordinates * 180. / np.pi
print('Preprocessing Coriolis')
corio = 2 * 7.292e-5 * np.sin(lonlat_coordinates[:,1])
slimPre.write_file(pre_data_dir+'coriolis.nc', region=region_global, time=None, data=[('coriolis', corio)])
#slimPre.netcdf_to_msh(mesh_file_name, pre_data_dir+'coriolis.nc', 'coriolis', 'coriolis')
print('Preprocessing Open Sea Boundary')
timeTide = slimPre.Time(initial_time=tides_initial_time, final_time=tides_final_time, time_step=tides_time_step)
region = slimPre.Region(mesh, physical_tags=tides_physical_tags)
(ssh,ssu,ssv) = slimPre.tpxo_tide(region, timeTide)
for i in range(len(timeTide._time)):
# from GMT to PST
timeTide._time[i] += time_zone
delta = (timeTide._time[i] - timeTide._time[0]) / (12 * 60. * 60.)
if delta > 1:
continue
if len(ssh) > 1:
ssh[i,:] *= max(0,min(delta,1))
ssu[i,:] *= max(0,min(delta,1))
ssv[i,:] *= max(0,min(delta,1))
slimPre.write_file(pre_data_dir+'tides.nc', region=region, time=timeTide, data=[('h', ssh),('u',ssu),('v',ssv)])
print('Preprocessing river discharge')
bath_file_name = slimPre.fetch_ftp('slim_data/columbia/forcings/beaver_ramp2.nc', pre_data_dir_base)
print('Preprocessing river salt')
slimPre.write_file(pre_data_dir+'river_salinity.nc', region=None, time=None, data=[('river_salinity', 0)])
print('Preprocessing open bnd salt')
slimPre.write_file(pre_data_dir+'open_salinity.nc', region=None, time=None, data=[('salt', 34)])
print('Preprocessing temp')
slimPre.write_file(pre_data_dir+'temp.nc', region=None, time=None, data=[('temp', 10)])
print('preprocessing done')
slimPre.exit(0)
import slim3d
import slimPre
import shutil
sim_Ti = '2006-05-1 00:00:00'
sim_Tf = '2006-05-14 00:30:00'
dt = 30.
sim_export = 900.
sim_exportFullRatio = -1#4*24 #each day
output_dir = 'output/'
pre_data_dir_base = 'data_'+slimPre.partition_nb()+'/'
pre_data_dir = pre_data_dir_base+slimPre.partition_id()+'/'
slimPre.make_directory(output_dir)
shutil.copyfile(pre_data_dir_base+'mesh3d.msh', output_dir+'mesh3d.msh')
shutil.copyfile(pre_data_dir_base+'mesh3d_show.msh', output_dir+'mesh3d_show.msh')
domain = slim3d.Domain(pre_data_dir_base+'mesh3d.msh')
equations = slim3d.Slim3d_equations(domain, salinity=True)
equations.set_coriolis((pre_data_dir+'coriolis.nc', 'coriolis'))
equations.set_horizontal_viscosity('smagorinsky', factor=0.5)
equations.set_implicit_vertical_diffusion(True)
equations.set_vertical_diffusivity('gotm')
equations.set_vertical_viscosity('gotm')
equations.set_limiter(True)
equations.set_bottom_friction(True, z0B=0.005, z0S=0.02)
equations.set_initial_temperature('netcdf', temperature=(pre_data_dir+'temp.nc', 'temp'))
equations.set_initial_salinity('netcdf', salinity=(pre_data_dir+'open_salinity.nc', 'salt'))
coast_tags = ["Land", "land", "island_0", "island_1", "island_2", "island_3",
"island_4", "island_5", "island_6", "island_7", "island_8", "island_9", "island_10"]
for coast in coast_tags:
equations.set_boundary_coast(coast)
equations.set_boundary_open('river',
flux=(pre_data_dir_base+'beaver_ramp2.nc', 'flux'),
salinity=(pre_data_dir+'river_salinity.nc', 'river_salinity'),
temperature=(pre_data_dir+'temp.nc', 'temp'))
open_tags = ["sea_west", "sea_north", "sea_south"]
for tag in open_tags:
equations.set_boundary_open(tag,
eta=(pre_data_dir+'tides.nc', 'h'),
u=(pre_data_dir+'tides.nc', 'u'),
v=(pre_data_dir+'tides.nc', 'v'),
temperature=(pre_data_dir+'temp.nc', 'temp'),
salinity=(pre_data_dir+'open_salinity.nc', 'salt'),
transport=True)
time_loop = slim3d.Loop(equations, time_step=dt,
export_time_step=sim_export,
ratio_full_export=sim_exportFullRatio,
initial_time=sim_Ti, final_time=sim_Tf,
output_directory=output_dir)
time_loop.export_elevation()
time_loop.export_salinity()
time_loop.export_uv()
time_loop.export_uv2d()
time_loop.loop()
\ No newline at end of file
......@@ -25,10 +25,6 @@ class Slim3d_equations :
slimSolver = self._slimSolver
slimSolver.setSolveS(salinity)
slimSolver.setSolveT(temperature)
slimSolver.setSolveUVImplicitVerticalDiffusion(False)
slimSolver.setSolveSImplicitVerticalDiffusion(False)
slimSolver.setSolveTImplicitVerticalDiffusion(False)
self._verticalDiffImplicit = False
slimSolver.setFlagUVLimiter(True)
if salinity:
slimSolver.setFlagSLimiter(True)
......@@ -63,21 +59,8 @@ class Slim3d_equations :
self._add_shear = None
self._verticalAdaptation = {}
def set_implicit_vertical_diffusion(self, flag):
slimSolver = self._slimSolver
self._verticalDiffImplicit = flag
if flag:
if self._vertical_viscosity:
slimSolver.setSolveUVImplicitVerticalDiffusion(True)
if slimSolver.getSolveS() and self._vertical_diffusivity:
slimSolver.setSolveSImplicitVerticalDiffusion(True)
if slimSolver.getSolveT() and self._vertical_diffusivity:
slimSolver.setSolveTImplicitVerticalDiffusion(True)
def set_implicit_vertical_advection(self, flag):
slimSolver = self._slimSolver
self._verticalAdvImplicit = flag
slimSolver.setSolveSImplicitVerticalAdvection(flag)
def set_implicit_vertical(self, flag):
self._slimSolver.setImplicitVertical(flag)
def set_linear_2d_equations(self, flag):
self._linear2d = flag
......@@ -136,8 +119,6 @@ class Slim3d_equations :
self._vertical_viscosity = mode
self._vertical_viscosity_constant_value = constant_value
self._vertical_viscosity_paca_values = paca_values
if self._verticalDiffImplicit:
self._slimSolver.setSolveUVImplicitVerticalDiffusion(True)
def set_vertical_diffusivity(self, mode, constant_value=1e-5, paca_values=[1e-5, 10]):
......@@ -146,10 +127,6 @@ class Slim3d_equations :
self._vertical_diffusivity = mode
self._vertical_diffusivity_constant_value = constant_value
self._vertical_diffusivity_paca_values = paca_values
if self._slimSolver.getSolveS() and self._vertical_diffusivity:
self._slimSolver.setSolveSImplicitVerticalDiffusion(True)
if self._slimSolver.getSolveT() and self._vertical_diffusivity:
self._slimSolver.setSolveTImplicitVerticalDiffusion(True)
def set_gotm_option_file(self,f):
self._gotm_option_file = f
......@@ -259,7 +236,7 @@ class Slim3d_equations :
class Loop:
"""Temporal solver"""
def __init__(self, equations, time_step=3600, export_time_step=3600, ratio_full_export=-1, initial_time='1970-01-01 00:00:00', final_time=0, output_directory="./output"):
def __init__(self, equations, time_step=3600, export_time_step=3600, ratio_full_export=-1, initial_time='1970-01-01 00:00:00', final_time=0, output_directory="./output", full_export_directory="./full_export"):
self._slimSolver = equations._slimSolver
self._equations = equations
self._dt = time_step
......@@ -275,6 +252,9 @@ class Loop:
else:
self.final_time = self.initial_time + final_time
self._odir = output_directory
self._full_odir = full_export_directory
self._restart_dir = None
self._restart_ind = -1
self._export_uv = False
self._export_w = False
self._export_z = True
......@@ -286,6 +266,7 @@ class Loop:
self._export_nuv = False
self._export_uv_nc = False
self._export_S_nc = False
self._export_T_nc = False
self._export_eta_nc = False
self._export_z_nc = False
self._export_rho = False
......@@ -293,8 +274,6 @@ class Loop:
self._timeIntegrator = None
self._export = []
self._ncExport = []
self._restart_dir = None
self._restart_ind = -1
self._timer = False
def export_uv(self, vect=True):
......
......@@ -128,8 +128,8 @@ def slim3d_setup(loop):
f.write(gotmFileStr)
f.close()
slimSolver.turbulenceSetupFile = 'gotmturb.nml'
if (slimSolver.getSolveUVImplicitVerticalDiffusion() or slimSolver.getSolveTurbulence()) and not slimSolver.getComputeBottomFriction():
dgpy.Msg.Fatal('If Vertical implicit diffusion or using gotm, bottom friction must be true')
if (slimSolver.getSolveTurbulence()) and not slimSolver.getComputeBottomFriction():
dgpy.Msg.Fatal('If using gotm, bottom friction must be true')
if eq._horizontal_viscosity == 'smagorinsky' :
slimSolver.setSolveUVGrad(True)
......@@ -144,7 +144,7 @@ def slim3d_setup(loop):
if eq._wind_stress[0] == 'speed':
eq._wind_speed_x = slim_private._load_function(eq._wind_stress[1], slimSolver.groups2d)
eq._wind_speed_y = slim_private._load_function(eq._wind_stress[2], slimSolver.groups2d)
f.windFunc = dgpy.functorNumpy(lambda cm : mergeUVNumpy(cm, eq.wind_speed_x, eq.wind_speed_y))
f.windFunc = dgpy.functorNumpy(lambda cm : mergeUVNumpy(cm, eq._wind_speed_x, eq._wind_speed_y))
f.windStressFunc = dgpy.windCubicStress(f.windFunc, eq._wind_stress[3])
# TODO Fatal if no wind speed and solving sediment
elif eq._wind_stress[0] == 'stress':
......@@ -273,13 +273,8 @@ def slim3d_setup(loop):
bottomTags = eq._domain._bottomTags
topTags = eq._domain._topTags
g = eq._gravity
# Open Boundaries
eta = f.eta2dFunc
uv = d.uvDof.getFunction()
uvAv2d = f.uvAv2dFunc
for openBnd in eq._boundary_open:
openBnd.setup()
......@@ -307,7 +302,10 @@ def slim3d_setup(loop):
horMomBndWall = horMomEq.newBoundaryWall()
horMomEq.addBoundaryCondition(eq._boundary_coast, horMomBndWall)
horMomEq.addBoundaryCondition('vertical_bottom', horMomBndWall)
horMomEq.addBoundaryCondition(bottomTags,horMomEq.newBoundaryBottom(f.z0BFunc, f.zBotDist2dPrecompFunc))
if slimSolver.getComputeBottomFriction() :
horMomEq.addBoundaryCondition(bottomTags,horMomEq.newBoundaryBottom(f.z0BFunc, f.zBotDist2dPrecompFunc))
else :
horMomEq.setBoundary0Flux(bottomTags)
wEq = e.wEq
wEq.setBoundarySymmetry(topTags)
......@@ -376,7 +374,9 @@ def slim3d_setup(loop):
timeIntegrator.setExportInterval(loop._export_dt)
#timeIntegrator.setExportAtEveryTimeStep(False)
#timeIntegrator.setFullExportAtEveryTimeStep(False)
timeIntegrator.setFullExportIntervalRatio(loop._ratio_full_export )
if loop._ratio_full_export > 0:
timeIntegrator.setFullExportIntervalRatio(loop._ratio_full_export )
timeIntegrator.initFullExporter(loop._full_odir)
timeIntegrator.setTimeStep(dt)
timeIntegrator.setInitTime(loop.initial_time)
timeIntegrator.setEndTime(loop.final_time)
......
......@@ -629,7 +629,7 @@ def residual_flow(region, flow, bathymetry, data_file_name=None):
def transportNumpy(cm) : # eta is neglected here
f = cm.get(bathFunc)*flow/section
return numpy.hstack([-f*nx,-f*ny])
return np.hstack([-f*nx,-f*ny])
f = dgpy.functorNumpy(transportNumpy)
vals = region._evaluateFunctor(f, 2)
if region._empty:
......
......@@ -1710,7 +1710,6 @@ static double _importIdx(dgDofContainer &dof, std::string filename, std::map<int
std::map<int, int> shiftedCMap;
for (auto it:componentMap)
shiftedCMap[it.first] = it.second - next_field;
double time;
int step;
_readMsh(dof, dirname + sname, time, step, shiftedCMap);
}
......
......@@ -6,6 +6,10 @@
#include <map>
#include "dgGroupOfElements.h"
#include "dgConfig.h"
#ifdef HAVE_MPI
#include "mpi.h"
#endif
/**Stores information about vertical column stucture.
* Each 3d element belongs to a certain column of stacked elements.
......@@ -109,14 +113,26 @@ class dgExtrusion {
size_t maxLev = 0;
for (size_t i = 0; i < getNbVerticals(false); i++)
maxLev = std::max(maxLev,getNbDGPointsInVertical(i));
return maxLev;
size_t nLevMax;
#ifdef HAVE_MPI
MPI_Allreduce(&maxLev, &nLevMax, 1, MPI_INT, MPI_MAX, MPI_COMM_WORLD);
#else
nLevMax = maxLev;
#endif
return nLevMax;
};
/** Return max number of CG nodes in any vertical edge **/
size_t getMaxNbCGPoints() const {
size_t maxLev = 0;
for (size_t i = 0; i < getNbVerticals(false); i++)
maxLev = std::max(maxLev,getNbCGPointsInVertical(i));
return maxLev;
size_t nLevMax;
#ifdef HAVE_MPI
MPI_Allreduce(&maxLev, &nLevMax, 1, MPI_INT, MPI_MAX, MPI_COMM_WORLD);
#else
nLevMax = maxLev;
#endif
return nLevMax;
};
/** Find 2d node for given 1D vertical.
* Given node iVertical returns (iGroup2d,iElem2d,iNode2d) corresponding to the 2d element group. */
......
......@@ -41,9 +41,6 @@ void dgConservationLawSW3dContinuity::riemann(functorCache &cache, fullMatrix<do
double g = slim3dParameters::g;
// determine which element is below
for(int i=0; i< val.size1(); i++) {
if (cache.interfaceGroup()->physicalTag() == "coastL"){
printf("wEq: eta %g -- %g, u %g -- %g h %g -- %g\n", etaIn(i,0), etaOut(i,0), UIn(i,0), UOut(i,0), bathL(i,0), bathR(i,0));
}
double nx = normals(i,0);
double ny = normals(i,1);
double nz = normals(i,2);
......
......@@ -97,11 +97,8 @@ void dgConservationLawSW3dTracer::interfaceTerm(functorCache &cm, fullMatrix<dou
}
void dgConservationLawSW3dTracer::diffusiveFlux(functorCache &cm, fullMatrix<double> &val) const {
const fullMatrix<double> &solGrad = cm.get(cm.solutionGradient());
const fullMatrix<double> &solGradE = cm.get(_explicitCGrad ? *_explicitCGrad : cm.solutionGradient());
const fullMatrix<double> &nuH = cm.get(*_kappaH);
const fullMatrix<double> &nuV = cm.get(*_kappaV);
int nQP = cm.nEvaluationPoint();
val.resize(nQP, 3, true);
if (_kappaH) {
const fullMatrix<double> &nuH = cm.get(*_kappaH);
const fullMatrix<double> &solGradE = cm.get(_explicitSolGrad ? *_explicitSolGrad : cm.solutionGradient());
......@@ -244,7 +241,7 @@ void dgConservationLawSW3dTracer::addIPTerm(functorCache &cm, fullMatrix<double>
solutionJumpPenalty += 0.5*ipf(i,0)*(cEL(i,0)-cER(i,0))*kappaH;
}
if (_kappaV) {
double kappaV = std::max((*kappaVL)(i,0),(*kappaVR)(i,0))*(ny*nz);
double kappaV = std::max((*kappaVL)(i,0),(*kappaVR)(i,0))*(nz*nz);
solutionJumpPenalty += 0.5*ipf(i,0)*(cL(i,0)-cR(i,0))*kappaV;
}
val(i,0) += -(meanNormalFlux + solutionJumpPenalty);
......
......@@ -54,12 +54,10 @@ void slim3dSolverDofs::initialize()
uvDof = new dgDofContainer(*groups3d, 2);
if ( _solver->_solveS || _solver->getSolveSImplicitVerticalDiffusion() ||
_solver->_solveTurbulence ) {
if ( _solver->_solveS || _solver->_solveTurbulence ) {
SDof = new dgDofContainer(*groups3d, 1);
}
if ( _solver->_solveT || _solver->getSolveTImplicitVerticalDiffusion() ||
_solver->_solveTurbulence ) {
if ( _solver->_solveT|| _solver->_solveTurbulence ) {
TDof = new dgDofContainer(*groups3d, 1);
}
bathDof2d = new dgDofContainer(*groups2d, 1);
......@@ -455,21 +453,9 @@ void slim3dSolverFunctions::initialize()
kappavFunc = functionSumNew(dofs->kappavDof->getFunction(),_kapv0Func);
}
}
else {
if ( !nuvFunc && _solver->getSolveUVImplicitVerticalDiffusion() ) {
Msg::Warning("nuvFunction not defined, using molecular diffusivity");
nuvFunc = _nuv0Func;
}
if ( !kappavFunc && (_solver->getSolveSImplicitVerticalDiffusion() ||
_solver->getSolveTImplicitVerticalDiffusion() ||
_solver->getSolveSedImplicitVerticalDiffusion() ) ) {
Msg::Warning("kappavFunction not defined, using molecular diffusivity");
kappavFunc = _kapv0Func;
}
}
if (_solver->getFlagJumpDiffusionTracerVert()){
if (_solver->_solveS){
if (_solver->_solveSDiff)
if (_solver->getImplicitVertical())
jumpDiffTracerVertS = new dgJumpDiffusionTracer(*_solver->groups3d, 0, dofs->wDof3d, 0., false, coeffCFLv, jumpPercentS, maxDiffS);
else
jumpDiffTracerVertS = new dgJumpDiffusionTracer(*_solver->groups3d, 0, dofs->wDof3d, approxDT, false, coeffCFLv, jumpPercentS, maxDiffS);
......@@ -481,7 +467,7 @@ void slim3dSolverFunctions::initialize()
kappavS = kappavJumpS;
}
if (_solver->_solveT){
if (_solver->_solveTDiff)
if (_solver->getImplicitVertical())
jumpDiffTracerVertT = new dgJumpDiffusionTracer(*_solver->groups3d, 0, dofs->wDof3d, 0., false, coeffCFLv, jumpPercentT, maxDiffT);
else
jumpDiffTracerVertT = new dgJumpDiffusionTracer(*_solver->groups3d, 0, dofs->wDof3d, approxDT, false, coeffCFLv, jumpPercentT, maxDiffT);
......@@ -834,28 +820,14 @@ void slim3dSolverEquations::initialize()
if ( _solver->_solveS ) {
SEq->setBathymetry(funcs->bathFunc2d);
if (funcs->kappahS)
SEq->setKappaH(funcs->kappahS);
if (!_solver->_solveSDiff && funcs->kappavS) {
Msg::Warning("Vertical diffusion equation of S disabled, assigning kappav to S equation");
SEq->setKappaV(funcs->kappavS);
}
}
if ( _solver->_solveSDiff ) {
//SDiffEq->setBoundary0Flux(_solver->physicalTags2d);
SEq->setKappaH(funcs->kappahS);
SEq->setKappaV(funcs->kappavS);
}
if ( _solver->_solveT ) {
TEq->setBathymetry(funcs->bathFunc2d);
if (funcs->kappahT)
TEq->setKappaH(funcs->kappahT);
if (!_solver->_solveTDiff && funcs->kappavT) {
Msg::Warning("Vertical diffusion equation of T disabled, assigning kappav to T equation");
TEq->setKappaV(funcs->kappavT);
}
}
if ( _solver->_solveTDiff ) {
//TDiffEq->setBoundary0Flux(_solver->physicalTags2d);
TEq->setKappaH(funcs->kappahT);
TEq->setKappaV(funcs->kappavT);
}
if ( _solver->_advectTKE ) {
......
......@@ -317,6 +317,7 @@ void slim3dNetCDFIO::exportUGRIDSingleComp(const dgDofContainer* dof, int iComp,
// construct data arrays
std::vector<double> dVec(nNode2d*nVert,-99.0);
std::vector<double> kBotVec(nNode2d,0.0);
size_t iPoint, iLevel, index2d;
size_t iGroup2d, iElem2d, iNode2d;
dgGroupCollection& groups = _columns->getGroups3d();
......@@ -330,13 +331,11 @@ void slim3dNetCDFIO::exportUGRIDSingleComp(const dgDofContainer* dof, int iComp,
_columns->mapVEdgeTo2dNode(iPoint,iGroup2d,iElem2d,iNode2d);
index2d = _ncWriter2d->get2dNodeLocalId(iGroup2d,iElem2d,iNode2d);
dVec[index2d*nVert+iLevel] = vals(iNode,iComp);
kBotVec[index2d] = _columns->getNbDGPointsInVertical(iPoint);
// dVec[iPoint*nVert+iLevel] = vals(iNode,iComp);
}
}
}
std::vector<double> kBotVec(nNode2d,0.0);
for (size_t iP = 0; iP < _columns->getNbVerticals(); iP++ )
kBotVec[iP] = _columns->getNbDGPointsInVertical(iP);
// write data array
MPI_Offset start[3], count[3];
start[0] = _ncWriter2d->_nodeOffset;
......
......@@ -7,11 +7,11 @@
slim3dSolver::slim3dSolver(const std::string meshFile3d, std::vector<std::string> bottomBndTags, std::vector<std::string> topBndTags, const std::string periodicMapFile) {
_solveW = true;
_solveUV = _solveUVDiff = true;
_solveUVAdvVert = false;
_solveS = _solveSDiff = _solveSAdvVert = false;
_solveT = _solveTDiff = _solveTAdvVert = false;
_solveSed = _solveSedDiff = _solveSedAdvVert = false;
_solveUV = true;
_implicitVertical = true;
_solveS = false;
_solveT = false;
_solveSed = false;
_numSedEq = 1;
_solveTurbulence = false;
_advectTKE = true;
......@@ -48,7 +48,6 @@ slim3dSolver::slim3dSolver(const std::string meshFile3d, std::vector<std::string
prismConverter = NULL;
movingMesh = NULL;
//zDiff1dEq = NULL;
uvLimiter = tkeLimiter = epsLimiter = sedLimiter= NULL;
jumpDiffTracerS = jumpDiffTracerT = jumpDiffTracerVertS = jumpDiffTracerVertT = NULL;
dgCG3d = dgCG2d = NULL;
copyField2d = NULL;
......@@ -121,10 +120,6 @@ slim3dSolver::~slim3dSolver()
if ( movingMesh ) delete movingMesh;
if ( _verticalAdaptMesh ) delete _verticalAdaptMesh;
// if ( zDiff1dEq ) delete zDiff1dEq;
if ( uvLimiter ) delete uvLimiter;
if ( tkeLimiter ) delete tkeLimiter;
if ( epsLimiter ) delete epsLimiter;
if ( sedLimiter ) delete sedLimiter;
if ( jumpDiffTracerS ) delete jumpDiffTracerS;
if ( jumpDiffTracerT ) delete jumpDiffTracerT;
if ( jumpDiffTracerVertS ) delete jumpDiffTracerVertS;
......@@ -161,21 +156,6 @@ slim3dSolverEquations* slim3dSolver::createEquations()
}
// if (_useAdaptiveVerticalGrid || _useAdaptiveVerticalGridWithDiffEq )
// zDiff1dEq = new slim3dZDiff1dEq(columnInfo, functions->getDt3d());
if (_useUVLimiter){
//uvLimiter = new dgSlopeLimiterVertex (_equations->horMomEq,groups3d, nodalVolume3d->get(), _useOptimalLimiter);
//uvLimiter = new dgSlopeLimiter4 (_equations->horMomEq,groups3d, _dofs->uvDof);
uvLimiter = new dgSlopeLimiterUV (_equations->horMomEq,groups3d);
//uvLimiter = new dgSlopeLimiterOriginal (_equations->horMomEq,groups3d);
}
if ( _solveSedDiff && !_solveSed && _useSedLimiter){
//sedLimiter = new dgSlopeLimiter4 (_equations->sedDiffEq,groups3d, _dofs->sedDof);
sedLimiter = new dgSlopeLimiterOriginal (_equations->sedEq[0],groups3d);
}
else if (_useSedLimiter && _solveSed){
//sedLimiter = new dgSlopeLimiter4 (_equations->sedEq,groups3d, _dofs->sedDof);
sedLimiter = new dgSlopeLimiterOriginal (_equations->sedEq[0],groups3d);
sedLimiterGround2d = new dgSlopeLimiterNoEq (groups2d);
}
if (getFlagJumpDiffusionTracer()){
jumpDiffTracerS = functions->jumpDiffTracerS;
jumpDiffTracerT = functions->jumpDiffTracerT;
......@@ -184,12 +164,6 @@ slim3dSolverEquations* slim3dSolver::createEquations()
jumpDiffTracerVertS = functions->jumpDiffTracerVertS;
jumpDiffTracerVertT = functions->jumpDiffTracerVertT;
}
if ( _advectTKE ) {
//tkeLimiter = new dgSlopeLimiter4 (_equations->tkeAdvEq,groups3d, _dofs->tkeDof);
tkeLimiter = new dgSlopeLimiterOriginal (_equations->tkeAdvEq,groups3d);
//epsLimiter = new dgSlopeLimiter4 (_equations->epsAdvEq,groups3d, _dofs->epsDof);
epsLimiter = new dgSlopeLimiterOriginal (_equations->epsAdvEq,groups3d);
}
_initialized = true;
return _equations;
}
......
......@@ -19,10 +19,11 @@ class slim3dExporter;
class slim3dSolver {
private:
bool _initialized;
bool _solveS, _solveSDiff, _solveSAdvVert;
bool _solveT, _solveTDiff, _solveTAdvVert;
bool _solveUV, _solveUVDiff, _solveUVAdvVert;
bool _solveSed, _solveSedDiff, _solveSedAdvVert;
bool _solveS;
bool _solveT;
bool _solveUV;
bool _solveSed;
bool _implicitVertical;
int _numSedEq;
bool _solveW;
bool _solveTurbulence, _advectTKE, _haveMovingMesh, _bottomFriction;
......@@ -72,7 +73,6 @@ public:
dgSW3dMovingMesh *movingMesh;
//slim3dZDiff1dEq *zDiff1dEq;
dgSW3dVerticalBottomRemover *verticalBottomRemover;
dgLimiter *uvLimiter, *tkeLimiter, *epsLimiter, *sedLimiter, *sedLimiterGround2d;
dgJumpDiffusionTracer *jumpDiffTracerS, *jumpDiffTracerT, *jumpDiffTracerVertS, *jumpDiffTracerVertT;
dgCGMeanFilter *dgCG3d, *dgCG2d;
// dgCGStructure *dgCG3d, *dgCG2d;
......@@ -122,29 +122,8 @@ public:
void setSolveUV(bool f) { _solveUV = f; };