diff --git a/swe_roe/ex1.c b/swe_roe/ex1.c index ef6d2c1..2ae6bde 100644 --- a/swe_roe/ex1.c +++ b/swe_roe/ex1.c @@ -14,16 +14,16 @@ struct _n_User { MPI_Comm comm; DM da; PetscInt Nt, Nx, Ny; - PetscReal dt, hx, hy; + PetscReal Lx, Ly, dt, dx, dy; PetscReal hu, hd; - PetscReal tiny_h; + PetscReal max_time, tiny_h; PetscInt dof, rank, size; Vec F, G, B; Vec subdomain; PetscInt xs, ys, xm, ym, xe, ye; PetscInt gxs, gxm, gys, gym, gxe, gye; - PetscBool debug, save, add_building; - PetscInt tstep; + PetscBool debug, add_building; + PetscInt save, tstep; }; extern PetscErrorCode RHSFunction(TS, PetscReal, Vec, Vec, void *); @@ -72,7 +72,7 @@ static PetscErrorCode SetInitialCondition(Vec X, User user) { PetscCall(VecZeroEntries(X)); for (PetscInt j = ys; j < ys + ym; j = j + 1) { for (PetscInt i = xs; i < xs + xm; i = i + 1) { - if (j < 95) { + if (j < 95 / user->dy) { x_ptr[j][i][0] = user->hu; } else { x_ptr[j][i][0] = user->hd; @@ -94,14 +94,11 @@ static PetscErrorCode SetInitialCondition(Vec X, User user) { PetscErrorCode Add_Buildings(User user) { PetscFunctionBeginUser; - DM da = user->da; - PetscReal hx = user->hx; - PetscReal hy = user->hy; - - PetscReal bu = 30 / hx; - PetscReal bd = 105 / hx; - PetscReal bl = 95 / hy; - PetscReal br = 105 / hy; + DM da = user->da; + PetscInt bu = 30 / user->dx; + PetscInt bd = 105 / user->dx; + PetscInt bl = 95 / user->dy; + PetscInt br = 105 / user->dy; PetscScalar ***b_ptr; PetscCall(DMDAVecGetArrayDOF(da, user->B, &b_ptr)); @@ -150,7 +147,10 @@ PetscErrorCode RHSFunction(TS ts, PetscReal t, Vec X, Vec F, void *ptr) { User user = (User)ptr; DM da = user->da; - PetscBool save = user->save; + PetscReal dx = user->dx; + PetscReal dy = user->dy; + PetscReal area = dx * dy; + PetscInt save = user->save; user->tstep = user->tstep + 1; @@ -188,7 +188,7 @@ PetscErrorCode RHSFunction(TS ts, PetscReal t, Vec X, Vec F, void *ptr) { for (PetscInt j = user->ys; j < user->ys + user->ym; j = j + 1) { for (PetscInt i = user->xs; i < user->xs + user->xm; i = i + 1) { for (PetscInt k = 0; k < user->dof; k = k + 1) { - f_ptr1[j][i][k] = -(f_ptr[j][i + 1][k] - f_ptr[j][i][k] + g_ptr[j + 1][i][k] - g_ptr[j][i][k]); + f_ptr1[j][i][k] = -(f_ptr[j][i + 1][k] * dx - f_ptr[j][i][k] * dx + g_ptr[j + 1][i][k] * dy - g_ptr[j][i][k] * dy) / area; } } } @@ -203,7 +203,7 @@ PetscErrorCode RHSFunction(TS ts, PetscReal t, Vec X, Vec F, void *ptr) { PetscCall(DMRestoreLocalVector(da, &localF)); PetscCall(DMRestoreLocalVector(da, &localG)); - if (save) { + if (save == 1) { char fname[PETSC_MAX_PATH_LEN]; sprintf(fname, "outputs/ex1_Nx_%d_Ny_%d_dt_%f_%d.dat", user->Nx, user->Ny, user->dt, user->tstep); @@ -550,8 +550,6 @@ PetscErrorCode solver(PetscReal hl, PetscReal hr, PetscReal ul, PetscReal ur, Pe } int main(int argc, char **argv) { - PetscErrorCode ierr; - /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - * * Initialize program and set problem parameters * - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ @@ -559,47 +557,47 @@ int main(int argc, char **argv) { User user; PetscCall(PetscNew(&user)); - user->dt = 0.04; - user->Nt = 180; - user->dof = 3; // h, uh, vh - user->comm = PETSC_COMM_WORLD; - user->debug = PETSC_FALSE; - user->debug = PETSC_FALSE; + user->dt = 0.04; + user->max_time = 7.2; + user->dof = 3; // h, uh, vh + user->comm = PETSC_COMM_WORLD; + user->debug = PETSC_FALSE; + user->save = 0; // save = 1: save outputs for each time step; save = 2: save outputs at last time step MPI_Comm_size(user->comm, &user->size); MPI_Comm_rank(user->comm, &user->rank); PetscPrintf(user->comm, "Running on %d processors! \n", user->size); /*Default number of cells (box), cellsize in x-direction, and y-direction */ - user->Nx = 200; - user->Ny = 200; - user->hx = 1.0; - user->hy = 1.0; + user->Lx = 200; + user->Ly = 200; + user->dx = 1.0; + user->dy = 1.0; user->hu = 10.0; // water depth for the upstream of dam [m] user->hd = 5.0; // water depth for the downstream of dam [m] user->tiny_h = 1e-7; - ierr = PetscOptionsBegin(user->comm, NULL, "2D Mesh Options", ""); - PetscCall(ierr); + PetscOptionsBegin(user->comm, NULL, "2D Mesh Options", ""); { - PetscCall(PetscOptionsInt("-Nx", "Number of cells in X", "", user->Nx, &user->Nx, NULL)); - PetscCall(PetscOptionsInt("-Ny", "Number of cells in Y", "", user->Ny, &user->Ny, NULL)); - PetscCall(PetscOptionsInt("-Nt", "Number of time steps", "", user->Nt, &user->Nt, NULL)); - PetscCall(PetscOptionsReal("-hx", "dx", "", user->hx, &user->hx, NULL)); - PetscCall(PetscOptionsReal("-hy", "dy", "", user->hy, &user->hy, NULL)); + PetscCall(PetscOptionsReal("-t", "simulation time", "", user->max_time, &user->max_time, NULL)); + PetscCall(PetscOptionsReal("-Lx", "Length in X", "", user->Lx, &user->Lx, NULL)); + PetscCall(PetscOptionsReal("-Ly", "Length in Y", "", user->Ly, &user->Ly, NULL)); + PetscCall(PetscOptionsReal("-dx", "dx", "", user->dx, &user->dx, NULL)); + PetscCall(PetscOptionsReal("-dy", "dy", "", user->dy, &user->dy, NULL)); PetscCall(PetscOptionsReal("-hu", "hu", "", user->hu, &user->hu, NULL)); PetscCall(PetscOptionsReal("-hd", "hd", "", user->hd, &user->hd, NULL)); PetscCall(PetscOptionsReal("-dt", "dt", "", user->dt, &user->dt, NULL)); PetscCall(PetscOptionsBool("-b", "Add buildings", "", user->add_building, &user->add_building, NULL)); PetscCall(PetscOptionsBool("-debug", "debug", "", user->debug, &user->debug, NULL)); - PetscCall(PetscOptionsBool("-save", "save outputs", "", user->save, &user->save, NULL)); + PetscCall(PetscOptionsInt("-save", "save outputs", "", user->save, &user->save, NULL)); } - ierr = PetscOptionsEnd(); - PetscCall(ierr); + PetscOptionsEnd(); assert(user->hu >= 0.); assert(user->hd >= 0.); - PetscReal max_time = user->Nt * user->dt; - PetscPrintf(user->comm, "Max simulation time is %f\n", max_time); + user->Nx = user->Lx / user->dx; + user->Ny = user->Ly / user->dy; + user->Nt = user->max_time / user->dt; + PetscPrintf(user->comm, "Max simulation time is %f\n", user->max_time); /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - * * Initialize DMDA @@ -608,7 +606,7 @@ int main(int argc, char **argv) { user->dof, 1, NULL, NULL, &user->da)); PetscCall(DMSetFromOptions(user->da)); PetscCall(DMSetUp(user->da)); - PetscCall(DMDASetUniformCoordinates(user->da, 0.0, user->Nx * user->hx, 0.0, user->Ny * user->hy, 0.0, 0.0)); + PetscCall(DMDASetUniformCoordinates(user->da, 0.0, user->Lx, 0.0, user->Ly, 0.0, 0.0)); /* { @@ -654,10 +652,12 @@ int main(int argc, char **argv) { char fname[PETSC_MAX_PATH_LEN]; sprintf(fname, "outputs/ex1_Nx_%d_Ny_%d_dt_%f_IC.dat", user->Nx, user->Ny, user->dt); - PetscViewer viewer; - PetscCall(PetscViewerBinaryOpen(user->comm, fname, FILE_MODE_WRITE, &viewer)); - PetscCall(VecView(X, viewer)); - PetscCall(PetscViewerDestroy(&viewer)); + if (user->save > 0) { + PetscViewer viewer; + PetscCall(PetscViewerBinaryOpen(user->comm, fname, FILE_MODE_WRITE, &viewer)); + PetscCall(VecView(X, viewer)); + PetscCall(PetscViewerDestroy(&viewer)); + } } /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - * @@ -681,7 +681,7 @@ int main(int argc, char **argv) { PetscCall(TSSetType(ts, TSEULER)); PetscCall(TSSetDM(ts, user->da)); PetscCall(TSSetRHSFunction(ts, R, RHSFunction, user)); - PetscCall(TSSetMaxTime(ts, max_time)); + PetscCall(TSSetMaxTime(ts, user->max_time)); PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER)); PetscCall(TSSetSolution(ts, X)); PetscCall(TSSetTimeStep(ts, user->dt)); @@ -696,6 +696,16 @@ int main(int argc, char **argv) { * - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ PetscCall(TSSolve(ts, X)); + if (user->save == 2) { + char fname[PETSC_MAX_PATH_LEN]; + sprintf(fname, "outputs/ex1_Nx_%d_Ny_%d_dt_%f_%d.dat", user->Nx, user->Ny, user->dt, user->Nt); + + PetscViewer viewer; + PetscCall(PetscViewerBinaryOpen(user->comm, fname, FILE_MODE_WRITE, &viewer)); + PetscCall(VecView(X, viewer)); + PetscCall(PetscViewerDestroy(&viewer)); + } + PetscCall(VecDestroy(&X)); PetscCall(VecDestroy(&R)); PetscCall(VecDestroy(&user->F)); diff --git a/swe_roe/ex1_code_verification.m b/swe_roe/ex1_code_verification.m new file mode 100644 index 0000000..93695bc --- /dev/null +++ b/swe_roe/ex1_code_verification.m @@ -0,0 +1,123 @@ +clear;close all;clc + +% Run several simulations before runing the following scripts +% mpiexec -n 40 ./ex1 -hd 5.0 -save 2 -b -dx 10 -dy 10 -dt 0.0005 +% mpiexec -n 40 ./ex1 -hd 5.0 -save 2 -b -dx 2 -dy 2 -dt 0.0005 +% mpiexec -n 40 ./ex1 -hd 5.0 -save 2 -b -dx 1 -dy 1 -dt 0.0005 +% mpiexec -n 40 ./ex1 -hd 5.0 -save 2 -b -dx 0.5 -dy 0.5 -dt 0.0005 +% mpiexec -n 40 ./ex1 -hd 5.0 -save 2 -b -dx 0.2 -dy 0.2 -dt 0.0005 +% mpiexec -n 40 ./ex1 -hd 5.0 -save 2 -b -dx 0.1 -dy 0.1 -dt 0.0005 +% mpiexec -n 40 ./ex1 -hd 5.0 -save 2 -b -dx 0.05 -dy 0.05 -dt 0.0005 +% mpiexec -n 40 ./ex1 -hd 5.0 -save 2 -b -dx 0.02 -dy 0.02 -dt 0.0005<-- assuming this is exact solution +% mpiexec -n 40 ./ex1 -hd 5.0 -save 2 -b -dx 0.01 -dy 0.01 -dt 0.0005<-- unstable, smaller time step is needed, but too expensive + +% User need to provide petsc directory +addpath('/Users/xudo627/developments/petsc/share/petsc/matlab/'); +Lx = 200; % [m] +Ly = 200; % [m] +dof = 3; + +filename = 'outputs/ex1_Nx_10000_Ny_10000_dt_0.000500_14400.dat'; +data = PetscBinaryRead(filename); +strs = strsplit(filename,'_'); +for i = 1 : length(strs) + if strcmp(strs{i},'Nx') + Nx = str2num(strs{i+1}); + elseif strcmp(strs{i},'Ny') + Ny = str2num(strs{i+1}); + elseif strcmp(strs{i},'dt') + dt = str2num(strs{i+1}); + end +end + +figure; +data = reshape(data, [dof, length(data)/dof]); +hexact = reshape(data(1,:),[Nx Ny]); +uexact = reshape(data(1,:),[Nx Ny]); +vexact = reshape(data(1,:),[Nx Ny]); + +imagesc(hexact); colorbar; caxis([0 10]); colormap(jet); + +Nxexact = 10000; + +dx = []; +Nxs = [20 40 100 200 400 1000 2000]; +err_L1 = NaN(length(Nxs),3); +err_L2 = NaN(length(Nxs),3); +err_Max = NaN(length(Nxs),3); + +k = 1; +for Nx = Nxs + disp(['Nx = ' num2str(Nx)]); + Ny = Nx; + dx = [dx; Lx/Nx]; + + bu = 29 / dx(k); + bd = 104 / dx(k); + bl = 94 / dx(k); + br = 106 / dx(k); + ii = 1 : Nx; jj = 1 : Ny; + + filename = ['outputs/ex1_Nx_' num2str(Nx) '_Ny_' num2str(Ny) '_dt_0.000500_14400.dat']; + data = PetscBinaryRead(filename); + data = reshape(data, [dof, length(data)/dof]); + h = reshape(data(1,:),[Nx Ny]); + u = reshape(data(1,:),[Nx Ny]); + v = reshape(data(1,:),[Nx Ny]); + + h(ii <= bu,jj < br & jj > bl) = NaN; + h(ii > bd ,jj < br & jj > bl) = NaN; + u(ii <= bu,jj < br & jj > bl) = NaN; + u(ii > bd ,jj < br & jj > bl) = NaN; + v(ii <= bu,jj < br & jj > bl) = NaN; + v(ii > bd ,jj < br & jj > bl) = NaN; + + scale_factor = Nxexact / Nx; + hexact_tmp = convert_res(hexact,1,scale_factor) ./ scale_factor^2; + uhexact_tmp = convert_res(uexact.*hexact,1,scale_factor) ./ scale_factor^2; + vhexact_tmp = convert_res(vexact.*hexact,1,scale_factor) ./ scale_factor^2; + + for i = 1 : 3 + if i == 1 + err = abs(h - hexact_tmp); + elseif i == 2 + err = abs(u.*h - uhexact_tmp); + elseif i == 3 + err = abs(v.*h - vhexact_tmp); + end + err_L1(k,i) = nansum(abs(err(:)))/(Nx*Ny); + err_L2(k,i) = sqrt(nansum(err(:).^2)/(Nx*Ny)); + err_Max(k,i) = max(abs(err(:))); + end + + k = k + 1; +end + +figure; set(gcf,'Position',[10 10 900 300]); +varnames = {'h','uh','vh'}; +for i = 1 : 3 + subplot(1,3,i); + loglog(dx, err_L1(:,i) ,'bo-','LineWidth',2); hold on; grid on; + loglog(dx, err_L2(:,i) ,'g*-','LineWidth',2); + loglog(dx, err_Max(:,i),'rx-','LineWidth',2); + ydum = ylim; + loglog([dx(end) dx(1)], ydum(1).*[1, dx(1)/dx(end)],'k--','LineWidth',2); + title(varnames{i},'FontSize',15,'FontWeight','bold'); + if i == 1 + leg = legend('1-norm','2-norm','max-norm','first order'); + end +end + + +function Z = convert_res(A,u,d) + %u = 5; %upsampling factor + %d = 12; %downsampling factor + t = d/u; %reduction factor + [m,n]=size(A); + L1=speye(m); L2=speye(round(m/t))/u; + R1=speye(n); R2=speye(round(n/t))/u; + L=repelem(L2,1,d) * repelem(L1,u,1); + R=(repelem(R2,1,d) * repelem(R1,u,1)).'; + Z= L*(A*R); +end + diff --git a/swe_roe/show_ex1.m b/swe_roe/ex1_show.m similarity index 83% rename from swe_roe/show_ex1.m rename to swe_roe/ex1_show.m index 2654481..904070a 100644 --- a/swe_roe/show_ex1.m +++ b/swe_roe/ex1_show.m @@ -3,6 +3,8 @@ % User need to provide petsc directory addpath('/Users/xudo627/developments/petsc/share/petsc/matlab/'); +Lx = 200; % [m] +Ly = 200; % [m] Nt = length(dir('outputs/ex1*.dat')) - 1; dof = 3; @@ -22,8 +24,11 @@ dt = str2num(strs{i+1}); end end -x = 1 : Nx; -y = 1 : Ny; + +dx = Lx / Nx; +dy = Ly / Ny; +x = dx : dx : Lx; +y = dy : dy : Ly; [x,y] = meshgrid(x,y); % Show Initial Condition @@ -36,8 +41,8 @@ % Animation of DAM break imAlpha = ones(size(h0)); -imAlpha(1:30,96:105) = 0; -imAlpha(106:200,96:105) = 0; +imAlpha(1:30/dx,95/dy+1:105/dy) = 0; +imAlpha(105/dy+1:200/dy,95/dy+1:105/dy) = 0; for i = 1 : Nt file = dir(['outputs/ex1_*_' num2str(i) '.dat']); data = PetscBinaryRead(fullfile(file(1).folder,file(1).name)); @@ -50,12 +55,12 @@ title(['t = ' num2str(i*dt) 's'],'FontSize',15); end -h(1:31,95:106) = NaN; -h(105:200,95:106) = NaN; -u(1:31,95:106) = 0; -u(105:200,95:106) = 0; -v(1:31,95:106) = 0; -v(105:200,95:106) = 0; +h(1:30/dx+1,95/dy:105/dy+1) = NaN; +h(105/dx:200/dy,95/dy:105/dy+1) = NaN; +u(1:30/dx+1,95/dy:105/dy+1) = 0; +u(105/dx:200/dx,95/dy:105/dy+1) = 0; +v(1:30/dx+1,95/dy:105/dy+1) = 0; +v(105/dx:200/dx,95/dy:105/dy+1) = 0; % Water depth and velocity at last time step figure; set(gcf,'Position',[10 10 1200 500]);