Commit 5397149f authored by Christopher Subich's avatar Christopher Subich

Merge branch 'master' into 'master'

Fix cfl calculation for mapped grids

Closes #6

See merge request !5
parents 9d7fa5dc a6e16637
......@@ -65,6 +65,7 @@ namespace NSIntegrator {
/* This will have to include ranges (arrays) for topography */
double Lx, Ly, Lz; // domain lengths
Array<double,1> Lzs; // domain depth at each horizontal position
double visco;
Control * usercode;
......@@ -141,6 +142,7 @@ namespace NSIntegrator {
tx(user->type_x()), ty(user->type_y()), tz(user->type_z()),
/* Lengths */
Lx(user->length_x()), Ly(user->length_y()), Lz(user->length_z()),
Lzs(split_range(szx)),
/* Viscosity */
visco(user->get_visco()),
/* User code */
......@@ -313,6 +315,9 @@ namespace NSIntegrator {
gradient->set_jac(secondDim,secondDim,bl_y/Ly);
// Calculate the vertical depths
blitz::Range all = blitz::Range::all();
Lzs = abs(zgrid(all,0,szz-1) - zgrid(all,0,0));
} else {
gradient->set_jac(firstDim,firstDim,bl_x/Lx);
gradient->set_jac(secondDim,secondDim,bl_y/Ly);
......
......@@ -802,7 +802,11 @@ namespace NSIntegrator {
grid dz is not a constant. So, use two terms of a Taylor
Series to give one */
double pin = M_PI/(szz-1);
timez = (Lz/2)*min(abs((pin*sin(kk*pin)+pin*pin/2*cos(kk*pin))/(1e-8 + abs(ws[0](ii,jj,kk)))));
if (mapped_problem) {
timez = min(abs(Lzs(ii)/2 * (pin*sin(kk*pin)+pin*pin/2*cos(kk*pin))/(1e-8 + abs(ws[0](ii,jj,kk)))));
} else {
timez = min(abs( Lz/2 * (pin*sin(kk*pin)+pin*pin/2*cos(kk*pin))/(1e-8 + abs(ws[0](ii,jj,kk)))));
}
} else {
timez = (Lz/szz) / (1e-6 + max(abs(ws[0])));
}
......
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