Commit b45bd1d2 authored by Derek Steinmoeller's avatar Derek Steinmoeller
Browse files

switch from forwards euler to 4th order low-storage explicit Runge-Kutta (LSERK4).

parent 7c023b0e
......@@ -137,6 +137,7 @@ int main(int argc, char ** argv) {
DTArray res_utheta(local_lbound,local_extent,local_storage);
DTArray res_ur(local_lbound,local_extent,local_storage);
DTArray res_eta(local_lbound,local_extent,local_storage);
// Allocate array for depth profile
DTArray H(local_lbound,local_extent,local_storage);
......@@ -287,63 +288,76 @@ int main(int argc, char ** argv) {
//update forcing:
Fx = 0.0; Fy = 0.0;
//step eta:
//eta_p = eta_m - dt*((eta_n u)_x + (eta_n v)_y)
temp = r2d*((H+eta)*ur);
mygrad.setup_array(&temp,FOURIER,NONE,CHEBY);
mygrad.get_dz(&rhs, false); //get r derivative, store in rhs
rhs = rinv*rhs; // (1/r)*d(r*F_r)/dr where F_r = (H+eta)u_r
temp = rinv*(H+eta)*utheta; // d((1/r)*(H+eta)*utheta)/dtheta
mygrad.setup_array(&temp,FOURIER,NONE,CHEBY);
mygrad.get_dx(&rhs, true); //get theta derivative. This time, we sum with old result (cumulative = true).
eta -= dt*rhs;
// filter free surface
if (FILTER_ON == true)
filter3(eta,XY_xform,FOURIER,NONE,CHEBY,F_CUTOFF,F_ORDER,F_STREN);
//impose explicit Neuman conditions on eta.
double neurhs1 = 0.0;
double neurhs2 = 0.0;
for (int j=local_lbound(1); j<(local_lbound(1)+local_extent(1)); j++) {
neurhs1 = sum(-Dr(0,Range(1,N_R-2))*eta(j,0,Range(1,N_R-2)));
neurhs2 = sum(-Dr(N_R-1,Range(1,N_R-2))*eta(j,0,Range(1,N_R-2)));
// perform 2x2 matrix product at each point
eta(j,0,0) = neumatinv(0,0)*neurhs1 + neumatinv(0,1)*neurhs2;
eta(j,0,N_R-1) = neumatinv(1,0)*neurhs1+ neumatinv(1,1)*neurhs2;
}
//form a1 and a2:
//nonlinear advection
mygrad.setup_array(&ur,FOURIER,NONE,CHEBY);
mygrad.get_dz(&temp,false); //d(ur)/dr
a1 = -ur*temp; // ur * d(ur)/dr
mygrad.get_dx(&temp,false); //d(ur)/d\theta
a1 -= rinv*utheta*temp; // (1/r)*utheta * d(ur)/dtheta
a1 += rinv*utheta*utheta; // Geometric term.
// a2
mygrad.setup_array(&utheta,FOURIER,NONE,CHEBY);
mygrad.get_dz(&temp,false); //d(utheta)/dr
a2 = -ur*temp; // ur * d(utheta)/dr
mygrad.get_dx(&temp,false); //d(utheta)/d\theta
a2 -= rinv*utheta*temp; // (1/r)(utheta) * d(utheta)/dtheta
a2 -= rinv*utheta*ur; // Geometric term.
// pressure gradient
mygrad.setup_array(&eta,FOURIER,NONE,CHEBY);
mygrad.get_dz(&temp,false); // d(eta)/dr
a1 -= G*temp;
mygrad.get_dx(&temp,false); // (1/r)*d(eta)_/dtheta
a2 -= G*rinv*temp;
// Coriolis and forcing
a1 += (ROT_F+ROT_B*r(kk)*sin(theta(ii)))*utheta + Fx;
a2 -= (ROT_F+ROT_B*r(kk)*sin(theta(ii)))*ur + Fy;
//done forming a1 & a2
// Initialize with zero RK residual.
res_utheta = 0.0; res_ur = 0.0; res_eta = 0.0;
for (int i=0; i < LSERK4::numStages; i++) {
//step eta:
//eta_p = eta_m - dt*((eta_n u)_x + (eta_n v)_y)
temp = -r2d*((H+eta)*ur);
mygrad.setup_array(&temp,FOURIER,NONE,CHEBY);
mygrad.get_dz(&rhs, false); //get r derivative, store in rhs
rhs = rinv*rhs; // (1/r)*d(r*F_r)/dr where F_r = (H+eta)u_r
temp = -rinv*(H+eta)*utheta; // d((1/r)*(H+eta)*utheta)/dtheta
mygrad.setup_array(&temp,FOURIER,NONE,CHEBY);
mygrad.get_dx(&rhs, true); //get theta derivative. This time, we sum with old result (cumulative = true).
res_eta = LSERK4::rk4a[i]*res_eta + dt*rhs;
eta += LSERK4::rk4b[i]*res_eta;
// filter free surface
if (FILTER_ON == true)
filter3(eta,XY_xform,FOURIER,NONE,CHEBY,F_CUTOFF,F_ORDER,F_STREN);
//impose explicit Neuman conditions on eta.
double neurhs1 = 0.0;
double neurhs2 = 0.0;
for (int j=local_lbound(1); j<(local_lbound(1)+local_extent(1)); j++) {
neurhs1 = sum(-Dr(0,Range(1,N_R-2))*eta(j,0,Range(1,N_R-2)));
neurhs2 = sum(-Dr(N_R-1,Range(1,N_R-2))*eta(j,0,Range(1,N_R-2)));
// perform 2x2 matrix product at each point
eta(j,0,0) = neumatinv(0,0)*neurhs1 + neumatinv(0,1)*neurhs2;
eta(j,0,N_R-1) = neumatinv(1,0)*neurhs1+ neumatinv(1,1)*neurhs2;
}
//form a1 and a2:
//nonlinear advection
mygrad.setup_array(&ur,FOURIER,NONE,CHEBY);
mygrad.get_dz(&temp,false); //d(ur)/dr
a1 = -ur*temp; // ur * d(ur)/dr
mygrad.get_dx(&temp,false); //d(ur)/d\theta
a1 -= rinv*utheta*temp; // (1/r)*utheta * d(ur)/dtheta
a1 += rinv*utheta*utheta; // Geometric term.
// a2
mygrad.setup_array(&utheta,FOURIER,NONE,CHEBY);
mygrad.get_dz(&temp,false); //d(utheta)/dr
a2 = -ur*temp; // ur * d(utheta)/dr
mygrad.get_dx(&temp,false); //d(utheta)/d\theta
a2 -= rinv*utheta*temp; // (1/r)(utheta) * d(utheta)/dtheta
a2 -= rinv*utheta*ur; // Geometric term.
// pressure gradient
mygrad.setup_array(&eta,FOURIER,NONE,CHEBY);
mygrad.get_dz(&temp,false); // d(eta)/dr
a1 -= G*temp;
mygrad.get_dx(&temp,false); // (1/r)*d(eta)_/dtheta
a2 -= G*rinv*temp;
// Coriolis and forcing
a1 += (ROT_F+ROT_B*r(kk)*sin(theta(ii)))*utheta + Fx;
a2 -= (ROT_F+ROT_B*r(kk)*sin(theta(ii)))*ur + Fy;
//done forming a1 & a2
// Perform Runge-Kutta incremental step..
res_ur = LSERK4::rk4a[i]*res_ur + dt*a1;
res_utheta = LSERK4::rk4a[i]*res_utheta + dt*a2;
ur += LSERK4::rk4b[i]*res_ur;
utheta += LSERK4::rk4b[i]*res_utheta;
} // RK-stages loop.
//compute - divergence of \vec{a}, store in rhs for helmholtz problem
......
Supports Markdown
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